Credit Risk & Risk Management Strategy¶

1. Load packages and dataset¶

Credit cards have become an essential part of the modern financial system, allowing customers to make purchases and payments in advance. However, this convenience is also associated with credit default risk, as some customers may spend beyond their repayment capacity and fail to pay back their debt on time.

This research is based on a dataset of 25,000 customers and aims to identify the trends and patterns behind credit card defaults. In addition, the project develops a predictive tool that enables users to estimate the probability of customer default in the following month with only a few clicks. The findings of this research also provide practical strategies for reducing default risk and improving credit risk management.

In [62]:
# Core Data and Numerical Tools
import numpy as np # scientific computing with arrays
import pandas as pd # data structures and analysis tools
from scipy.stats import randint, uniform # statistical distributions for hyperparameter sampling
import random

# Visualization Tools
import matplotlib.pyplot as plt # general purpose data visualization
import seaborn as sns # statistical data visualization

# Preprocessing and Utilities
from sklearn.model_selection import train_test_split # split data into train and validation subsets
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder # standardize features and encode categorical features
from sklearn.compose import ColumnTransformer # column transformation
from sklearn.pipeline import Pipeline # data processing and modeling steps
from sklearn.feature_selection import VarianceThreshold # remove features with low variance
from collections import Counter

# Models
from sklearn.tree import plot_tree, DecisionTreeClassifier # Decision Tree
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier # Gradient Boost and Random Forest
from xgboost import XGBClassifier # XGBoost
from xgboost import plot_importance # visualization of feature importance scores
from sklearn.model_selection import GridSearchCV, StratifiedKFold, RandomizedSearchCV # tuning and cross-validation tools
import numpy as np
# Evaluation Metrics and Display
from sklearn.metrics import mean_squared_error, roc_curve, auc, r2_score, accuracy_score, roc_auc_score, precision_score, recall_score, f1_score, make_scorer # model performance evaluation
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, RocCurveDisplay, brier_score_loss # output analysis and display
from sklearn.calibration import calibration_curve # probability calibration
from sklearn.linear_model import LogisticRegression # probability calibration
# Save model for UI design
import joblib
In [11]:
df = pd.read_csv(r"C:\Users\nguye\Downloads\train_dataset_creditscoring.csv")

2. Data exploration¶

In [12]:
df.shape
Out[12]:
(25247, 27)
In [13]:
df.describe()
Out[13]:
Customer_ID marriage sex education LIMIT_BAL age pay_0 pay_2 pay_3 pay_4 ... Bill_amt6 pay_amt1 pay_amt2 pay_amt3 pay_amt4 pay_amt5 pay_amt6 AVG_Bill_amt PAY_TO_BILL_ratio next_month_default
count 25247.000000 25247.000000 25247.000000 25247.000000 25247.000000 25121.000000 25247.000000 25247.000000 25247.000000 25247.000000 ... 25247.000000 25247.000000 2.524700e+04 25247.000000 25247.000000 25247.000000 25247.000000 25247.000000 25247.000000 25247.000000
mean 17640.000000 1.551907 0.604111 1.852101 168342.060443 35.438199 -0.042857 -0.159544 -0.190359 -0.241415 ... 38806.221029 5718.624966 6.047352e+03 5288.910651 4865.960834 4906.766828 5270.499287 44859.647485 0.362962 0.190399
std 7288.325459 0.522629 0.489050 0.797379 129892.784807 9.174998 1.099315 1.173990 1.172636 1.146753 ... 59182.792531 16806.842125 2.400962e+04 17851.879609 15979.116544 15860.726852 17960.816915 62819.226119 5.047206 0.392624
min 5017.000000 0.000000 0.000000 0.000000 10000.000000 21.000000 -2.000000 -2.000000 -2.000000 -2.000000 ... 0.000000 0.000000 0.000000e+00 0.000000 0.000000 0.000000 0.000000 -56043.170000 -546.930000 0.000000
25% 11328.500000 1.000000 0.000000 1.000000 50000.000000 28.000000 -1.000000 -1.000000 -1.000000 -1.000000 ... 1241.710000 999.985000 9.219100e+02 399.990000 300.150000 262.365000 130.070000 4858.670000 0.040000 0.000000
50% 17640.000000 2.000000 1.000000 2.000000 140000.000000 34.000000 0.000000 0.000000 0.000000 0.000000 ... 17102.580000 2145.020000 2.026830e+03 1844.300000 1500.100000 1513.790000 1500.040000 21102.830000 0.090000 0.000000
75% 23951.500000 2.000000 1.000000 2.000000 240000.000000 41.000000 0.000000 0.000000 0.000000 0.000000 ... 49245.195000 5031.150000 5.000190e+03 4600.640000 4014.990000 4099.890000 4018.780000 57136.580000 0.590000 0.000000
max 30263.000000 3.000000 1.000000 6.000000 1000000.000000 79.000000 8.000000 8.000000 8.000000 7.000000 ... 961663.620000 873551.980000 1.684259e+06 896040.150000 621000.080000 426529.180000 528666.150000 877313.830000 205.380000 1.000000

8 rows × 27 columns

In [93]:
# Histograms
df_plot = df.drop(columns=['Customer_ID']) # Drop ID
axes = df_plot.hist(bins = 30, figsize = (20,15),
             xlabelsize=14,   # x-axis label size
            ylabelsize=14)   # y-axis label size

for ax in axes.ravel():
    ax.set_title(ax.get_title(), fontsize=14)   # subplot title
    ax.tick_params(axis='both', labelsize=14)   # tick labels

# Increase distance between subplots (default is around 0.2)
plt.subplots_adjust(hspace=0.5, wspace=0.3)

plt.show()
In [15]:
df.dtypes
Out[15]:
Customer_ID             int64
marriage                int64
sex                     int64
education               int64
LIMIT_BAL               int64
age                   float64
pay_0                   int64
pay_2                   int64
pay_3                   int64
pay_4                   int64
pay_5                   int64
pay_6                   int64
Bill_amt1             float64
Bill_amt2             float64
Bill_amt3             float64
Bill_amt4             float64
Bill_amt5             float64
Bill_amt6             float64
pay_amt1              float64
pay_amt2              float64
pay_amt3              float64
pay_amt4              float64
pay_amt5              float64
pay_amt6              float64
AVG_Bill_amt          float64
PAY_TO_BILL_ratio     float64
next_month_default      int64
dtype: object

3. Data wrangling¶

3.1. Checking missing data & filling missing data¶

In [16]:
# Checking missing data
df.isnull().sum()
Out[16]:
Customer_ID             0
marriage                0
sex                     0
education               0
LIMIT_BAL               0
age                   126
pay_0                   0
pay_2                   0
pay_3                   0
pay_4                   0
pay_5                   0
pay_6                   0
Bill_amt1               0
Bill_amt2               0
Bill_amt3               0
Bill_amt4               0
Bill_amt5               0
Bill_amt6               0
pay_amt1                0
pay_amt2                0
pay_amt3                0
pay_amt4                0
pay_amt5                0
pay_amt6                0
AVG_Bill_amt            0
PAY_TO_BILL_ratio       0
next_month_default      0
dtype: int64
In [17]:
# Filling missing ages with the median age of the dataset
median_age = df[["age"]].median()
print(median_age)

age    34.0
dtype: float64
In [18]:
df_train = df.fillna(median_age)
In [19]:
# Check again
df_train.isnull().sum()
Out[19]:
Customer_ID           0
marriage              0
sex                   0
education             0
LIMIT_BAL             0
age                   0
pay_0                 0
pay_2                 0
pay_3                 0
pay_4                 0
pay_5                 0
pay_6                 0
Bill_amt1             0
Bill_amt2             0
Bill_amt3             0
Bill_amt4             0
Bill_amt5             0
Bill_amt6             0
pay_amt1              0
pay_amt2              0
pay_amt3              0
pay_amt4              0
pay_amt5              0
pay_amt6              0
AVG_Bill_amt          0
PAY_TO_BILL_ratio     0
next_month_default    0
dtype: int64

3.2. Data wrangling¶

In [20]:
# Define wrangling steps
def standardization(df_train):
    # Late_payment_count
    df_train["late_payment_count"] = (
        (df_train[[f"pay_{i}" for i in [0, 2, 3, 4, 5, 6]]] >= 1)
        .any(axis=1)
        .astype(int)
    )

    #df_train["minimum_payment_count"] = (
    #    (df_train[[f"pay_{i}" for i in [0, 2, 3, 4, 5, 6]]] == 0)
    #    .any(axis=1)
    #    .astype(int)
    #)

    # Bill trend
    for i in range(1, 6):

        df_train[f"bill_trend{i}"] = np.where(
            df_train[f"Bill_amt{i}"] <= 0,
            1,  # 100%
            (
                (df_train[f"Bill_amt{i+1}"] - df_train[f"Bill_amt{i}"])
                / df_train[f"Bill_amt{i}"]
            )
        )

    df_train["average_bill_trend"] = df_train[
        [f"bill_trend{i}" for i in range(1, 6)]
    ].mean(axis=1)

    # Average payment
    df_train["average_payment"] = df_train[
        [f"pay_amt{i}" for i in range(1, 6)]
    ].mean(axis=1)

    df_train["avg_monthly_pay_to_bill"] = np.where(
        df_train[f"AVG_Bill_amt"] <= 0, 1,
        df_train["average_payment"] / (df_train["AVG_Bill_amt"])
    )

    # Pay to bill ratio
    for i in range(2, 7):

        df_train[f"pay_to_bill_month{i}"] = np.where(
            df_train[f"Bill_amt{i-1}"] <= 0,
            1,
            df_train[f"pay_amt{i}"] / df_train[f"Bill_amt{i-1}"]
        )

    # Trend of pay-to-bill
    for i in range(2, 6):

        df_train[f"pay_to_bill_trend{i}"] = np.where(
            df_train[f"pay_to_bill_month{i}"] == 0,
            1,
            (
                (
                    df_train[f"pay_to_bill_month{i+1}"]
                    - df_train[f"pay_to_bill_month{i}"]
                )
                / df_train[f"pay_to_bill_month{i}"]
            )
        )

    df_train["average_pay_to_bill_trend"] = df_train[
        [f"pay_to_bill_trend{i}" for i in range(2, 6)]
    ].mean(axis=1)


    df_train['credit_utilization'] = (
        df_train['AVG_Bill_amt'] / (df_train['LIMIT_BAL'] + 1)
    )

    return df_train
In [21]:
df_train_clean = standardization(df_train)
df_train_clean.shape
Out[21]:
(25247, 47)

3.3. Features defining¶

In [22]:
# Defines features
X = df_train_clean.drop(columns = ["next_month_default", 'Customer_ID', #'pay_2', 'pay_3', 'pay_4', 'pay_5', 'pay_6',
        'Bill_amt1', 'Bill_amt2', 'Bill_amt3',
        'Bill_amt4', 'Bill_amt5', 'Bill_amt6',
        'pay_amt1', 'pay_amt2', 'pay_amt3',
        'pay_amt4', 'pay_amt5', 'pay_amt6',
         'PAY_TO_BILL_ratio', 'bill_trend1', 'bill_trend2', 'bill_trend3', 'bill_trend4',
       'bill_trend5', 'pay_to_bill_month2',
       'pay_to_bill_month3', 'pay_to_bill_month4', 'pay_to_bill_month5',
        'pay_to_bill_trend2', 'pay_to_bill_trend3',
       'pay_to_bill_trend4', 'pay_to_bill_month6', 'pay_to_bill_trend5','bill_trend1', 'bill_trend2', 'bill_trend3',
       'bill_trend4', 'bill_trend5',])
X.columns
Out[22]:
Index(['marriage', 'sex', 'education', 'LIMIT_BAL', 'age', 'pay_0', 'pay_2',
       'pay_3', 'pay_4', 'pay_5', 'pay_6', 'AVG_Bill_amt',
       'late_payment_count', 'average_bill_trend', 'average_payment',
       'avg_monthly_pay_to_bill', 'average_pay_to_bill_trend',
       'credit_utilization'],
      dtype='str')
In [23]:
# Check correlation to avoid multicollinearity
X.corr()
Out[23]:
marriage sex education LIMIT_BAL age pay_0 pay_2 pay_3 pay_4 pay_5 pay_6 AVG_Bill_amt late_payment_count average_bill_trend average_payment avg_monthly_pay_to_bill average_pay_to_bill_trend credit_utilization
marriage 1.000000 -0.031876 -0.149530 -0.102085 -0.413662 0.020694 0.023440 0.031506 0.031848 0.031855 0.031112 -0.023068 0.002214 -0.006318 -0.008207 -0.006554 0.006177 0.048741
sex -0.031876 1.000000 0.019173 0.023435 -0.089775 -0.051158 -0.066689 -0.060550 -0.060808 -0.052867 -0.043636 -0.023827 -0.021602 0.005914 -0.004793 0.002711 0.009040 -0.065878
education -0.149530 0.019173 1.000000 -0.220720 0.174701 0.105422 0.126021 0.113370 0.108537 0.098257 0.081798 0.007969 0.023707 -0.003717 -0.059183 -0.045693 0.002560 0.171771
LIMIT_BAL -0.102085 0.023435 -0.220720 1.000000 0.142283 -0.268792 -0.294822 -0.283323 -0.264973 -0.246225 -0.230497 0.302158 -0.208700 0.011631 0.324793 0.061479 -0.002501 -0.385545
age -0.413662 -0.089775 0.174701 0.142283 1.000000 -0.041480 -0.053802 -0.053781 -0.046872 -0.052584 -0.047208 0.055521 -0.019785 0.007701 0.039095 0.012034 0.009523 -0.040146
pay_0 0.020694 -0.051158 0.105422 -0.268792 -0.041480 1.000000 0.667648 0.569559 0.535307 0.507752 0.470861 0.192420 0.606112 -0.024968 -0.106446 -0.112753 -0.013359 0.415533
pay_2 0.023440 -0.066689 0.126021 -0.294822 -0.053802 0.667648 1.000000 0.766510 0.662837 0.624032 0.575934 0.239527 0.405748 -0.036673 -0.087358 -0.164790 -0.013115 0.502862
pay_3 0.031506 -0.060550 0.113370 -0.283323 -0.053781 0.569559 0.766510 1.000000 0.775747 0.686512 0.630682 0.238982 0.420419 -0.040655 -0.064298 -0.153145 -0.017331 0.491577
pay_4 0.031848 -0.060808 0.108537 -0.264973 -0.046872 0.535307 0.662837 0.775747 1.000000 0.818685 0.711637 0.249515 0.375933 -0.026633 -0.043284 -0.145403 -0.012109 0.494344
pay_5 0.031855 -0.052867 0.098257 -0.246225 -0.052584 0.507752 0.624032 0.686512 0.818685 1.000000 0.812498 0.262093 0.343576 -0.023802 -0.023273 -0.143899 -0.013044 0.492094
pay_6 0.031112 -0.043636 0.081798 -0.230497 -0.047208 0.470861 0.575934 0.630682 0.711637 0.812498 1.000000 0.268905 0.351793 -0.015388 -0.004260 -0.138060 -0.011486 0.486606
AVG_Bill_amt -0.023068 -0.023827 0.007969 0.302158 0.055521 0.192420 0.239527 0.238982 0.249515 0.262093 0.268905 1.000000 -0.062796 -0.009848 0.331930 -0.111962 -0.002040 0.547846
late_payment_count 0.002214 -0.021602 0.023707 -0.208700 -0.019785 0.606112 0.405748 0.420419 0.375933 0.343576 0.351793 -0.062796 1.000000 -0.007980 -0.127197 -0.016517 -0.009667 0.128258
average_bill_trend -0.006318 0.005914 -0.003717 0.011631 0.007701 -0.024968 -0.036673 -0.040655 -0.026633 -0.023802 -0.015388 -0.009848 -0.007980 1.000000 0.041313 0.008161 0.024201 -0.020496
average_payment -0.008207 -0.004793 -0.059183 0.324793 0.039095 -0.106446 -0.087358 -0.064298 -0.043284 -0.023273 -0.004260 0.331930 -0.127197 0.041313 1.000000 0.043223 0.017736 0.039664
avg_monthly_pay_to_bill -0.006554 0.002711 -0.045693 0.061479 0.012034 -0.112753 -0.164790 -0.153145 -0.145403 -0.143899 -0.138060 -0.111962 -0.016517 0.008161 0.043223 1.000000 0.001811 -0.172906
average_pay_to_bill_trend 0.006177 0.009040 0.002560 -0.002501 0.009523 -0.013359 -0.013115 -0.017331 -0.012109 -0.013044 -0.011486 -0.002040 -0.009667 0.024201 0.017736 0.001811 1.000000 -0.005907
credit_utilization 0.048741 -0.065878 0.171771 -0.385545 -0.040146 0.415533 0.502862 0.491577 0.494344 0.492094 0.486606 0.547846 0.128258 -0.020496 0.039664 -0.172906 -0.005907 1.000000
In [24]:
# Define target variable
y = df_train_clean["next_month_default"]

4. Data visualization¶

In [25]:
fig = plt.figure(figsize=(10, 7))

plt.pie(
    df_train_clean["next_month_default"].value_counts(),
    labels=["Not Defaulted", "Defaulted"],
    autopct='%1.1f%%',      # show percentage
    startangle=90,          # rotate for nicer view
    pctdistance=0.75,       # position of percentage text
    labeldistance=1.1       # position of labels
)

plt.title("Distribution of Target Variable")
plt.legend(["Not Defaulted", "Defaulted"])
plt.show()

Out of the 25,000 customers in the dataset, 19% defaulted within one month, indicating a class imbalance that must be taken into account during the modelling process.

In [26]:
# Distribution of age (box plot)

for x in ["marriage", "sex", "education", "late_payment_count"]:
    boxplot = sns.boxplot(x="next_month_default", y=x, data=df_train_clean)
    plt.show()

No significant differences were detected between the default and non-default groups with regard to education, marital status, and sex. However, customers who are unable to repay in the following month tend to have a history of late payments.

In [27]:
for x in ["avg_monthly_pay_to_bill", "average_bill_trend", "credit_utilization"]:
    boxplot = sns.boxplot(x="next_month_default", y=x, data=df_train_clean)
    plt.yscale("log")  # Set x-axis limits to focus on the box plot
    plt.show()

According to the box plots, customers who default within one month tend to make lower payments relative to their bill amounts and also tend to utilize a larger proportion of their credit limits.

In [28]:
for x in ["age", "LIMIT_BAL", "AVG_Bill_amt"]:
    plt.hist(df_train_clean[df_train_clean["next_month_default"] == 1][x], bins=30, color = "orange", alpha=0.5)
    plt.hist(df_train_clean[df_train_clean["next_month_default"] == 0][x], bins=30, color = "blue", alpha=0.2)
    plt.xlabel(x)
    plt.ylabel("Frequency")
    plt.title(f"Distribution of {x} for Defaulted Customers")
    plt.show()

However, customers who are likely to default and those who are not show no substantial differences in age, credit limit, or average bill amount.

5. Modelling¶

Preprocessing¶

In this section, we develop two models: XGBoost and Random Forest to predict customer default in the following month. After choosing the best performing model, the model will be calibrated and used to develope the credit scoring tool. The F1-score is used as the primary evaluation metric because it balances precision and recall for the imbalanced target. Model development and threshold selection use the training/validation data, while the final test set is kept separate for performance assessment. Because the selected tree-based model uses class weights, its raw output probabilities are calibrated separately before they are interpreted as Probability of Default (PD).

In [29]:
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
In [30]:
# Split the data into training, calibration, validation, and final test sets.
# The calibration set is kept separate from model fitting so that probability calibration
# does not use the same observations used to train the selected tree-based model.
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, stratify=y, test_size=0.40, random_state=42
)

# 10% calibration, 10% validation, 20% final test
X_cal, X_holdout, y_cal, y_holdout = train_test_split(
    X_temp, y_temp, stratify=y_temp, test_size=0.75, random_state=42
)
X_val, X_test, y_val, y_test = train_test_split(
    X_holdout, y_holdout, stratify=y_holdout, test_size=0.67, random_state=42
)

print(f"Training set: {len(X_train):,}")
print(f"Calibration set: {len(X_cal):,}")
print(f"Validation set: {len(X_val):,}")
print(f"Test set: {len(X_test):,}")

Training set: 15,148
Calibration set: 2,524
Validation set: 2,499
Test set: 5,076
In [31]:
features = X.columns

First, a preprocessing pipeline is defined for the models. For the tree-based models, class imbalance is addressed and variables showing little or no difference between the default and non-default groups are removed.

In [32]:
# Define preprocessing step
trees_preprocessor = ColumnTransformer(
    [
        ("zv_filter", VarianceThreshold(1e-4), features)
    ],
    remainder='passthrough' # Pass through numerical features unscaled
)
In [33]:
# Balance the training data without requiring imbalanced-learn.
# Only the training set is oversampled; calibration, validation, and test sets
# retain their original default rate.
print("Before balancing:", Counter(y_train))

rng = np.random.default_rng(42)
train_df = pd.concat([X_train, y_train.rename("target")], axis=1)
majority = train_df[train_df["target"] == 0]
minority = train_df[train_df["target"] == 1]

minority_oversampled = minority.sample(
    n=len(majority),
    replace=True,
    random_state=42
)

balanced_df = pd.concat([majority, minority_oversampled]).sample(
    frac=1,
    random_state=42
)

X_balanced = balanced_df.drop(columns="target")
y_balanced = balanced_df["target"]

print("After balancing:", Counter(y_balanced))

Before balancing: Counter({0: 12264, 1: 2884})
After balancing: Counter({1: 12264, 0: 12264})

XGBoost¶

In [34]:
# 1. XGBoost classifier
#cale_pos_weight = (len(y_balanced) - sum(y_balanced))/sum(y_balanced) # handle class imbalance

xgb_model = XGBClassifier(
    objective="binary:logistic",
    #scale_pos_weight=scale_pos_weight,
    eval_metric="aucpr",
    use_label_encoder=False
)
In [35]:
# 2. Pipeline
pipeline = Pipeline(steps=
 [
    ("preprocess", trees_preprocessor),
    ("model", xgb_model)
])
In [36]:
# 3. Hyperparameter
params = {
    "model__n_estimators": [300, 500, 700],
    #"model__scale_pos_weight": [scale_pos_weight],
    "model__max_depth": [3, 5, 7],
    "model__min_child_weight": [3, 5, 7],
    "model__gamma": [0, 0.1, 0.3, 0.5],
    "model__colsample_bytree": [0.6, 0.8, 1.0],
    "model__learning_rate": [0.01, 0.05, 0.1]
}
In [37]:
# 4. Cross-validation
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=9)
#scorer = make_scorer(roc_auc_score, needs_proba=True)
In [38]:
# 5. Hyperparameter tuning
search = RandomizedSearchCV(
    estimator=pipeline,
    param_distributions=params,
    n_iter=20,  # like size=20 in R
    scoring='f1',
    cv=cv,
    verbose=2,
    n_jobs=-1,
    random_state=42
)
In [39]:
# 6. Fit model
search.fit(X_balanced, y_balanced)

Fitting 10 folds for each of 20 candidates, totalling 200 fits

c:\Users\nguye\AppData\Local\Programs\Python\Python311\Lib\site-packages\xgboost\training.py:200: UserWarning: [20:37:31] WARNING: C:\actions-runner\_work\xgboost\xgboost\src\learner.cc:782: 
Parameters: { "use_label_encoder" } are not used.

  bst.update(dtrain, iteration=i, fobj=obj)
Out[39]:
RandomizedSearchCV(cv=StratifiedKFold(n_splits=10, random_state=9, shuffle=True),
                   estimator=Pipeline(steps=[('preprocess',
                                              ColumnTransformer(remainder='passthrough',
                                                                transformers=[('zv_filter',
                                                                               VarianceThreshold(threshold=0.0001),
                                                                               Index(['marriage', 'sex', 'education', 'LIMIT_BAL', 'age', 'pay_0', 'pay_2',
       'pay_3', 'pay_4', 'pay_5', 'pay_6', 'AVG_Bill_amt',
       'late_pay...
                                                            multi_strategy=None,
                                                            n_estimators=None,
                                                            n_jobs=None,
                                                            num_parallel_tree=None, ...))]),
                   n_iter=20, n_jobs=-1,
                   param_distributions={'model__colsample_bytree': [0.6, 0.8,
                                                                    1.0],
                                        'model__gamma': [0, 0.1, 0.3, 0.5],
                                        'model__learning_rate': [0.01, 0.05,
                                                                 0.1],
                                        'model__max_depth': [3, 5, 7],
                                        'model__min_child_weight': [3, 5, 7],
                                        'model__n_estimators': [300, 500, 700]},
                   random_state=42, scoring='f1', verbose=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Parameters
estimator estimator: estimator object

An object of that type is instantiated for each grid point.
This is assumed to implement the scikit-learn estimator interface.
Either estimator needs to provide a ``score`` function,
or ``scoring`` must be passed. 		Pipeline(step...=None, ...))])
param_distributions param_distributions: dict or list of dicts

Dictionary with parameters names (`str`) as keys and distributions
or lists of parameters to try. Distributions must provide a ``rvs``
method for sampling (such as those from scipy.stats.distributions).
If a list is given, it is sampled uniformly.
If a list of dicts is given, first a dict is sampled uniformly, and
then a parameter is sampled using that dict as above. 		{'model__colsample_bytree': [0.6, 0.8, ...], 'model__gamma': [0, 0.1, ...], 'model__learning_rate': [0.01, 0.05, ...], 'model__max_depth': [3, 5, ...], ...}
n_iter n_iter: int, default=10

Number of parameter settings that are sampled. n_iter trades
off runtime vs quality of the solution. 		20
scoring scoring: str, callable, list, tuple or dict, default=None

Strategy to evaluate the performance of the cross-validated model on
the test set.

If `scoring` represents a single score, one can use:

- a single string (see :ref:`scoring_string_names`);
- a callable (see :ref:`scoring_callable`) that returns a single value;
- `None`, the `estimator`'s
:ref:`default evaluation criterion ` is used.

If `scoring` represents multiple scores, one can use:

- a list or tuple of unique strings;
- a callable returning a dictionary where the keys are the metric
names and the values are the metric scores;
- a dictionary with metric names as keys and callables as values.

See :ref:`multimetric_grid_search` for an example.

If None, the estimator's score method is used. 		'f1'
n_jobs n_jobs: int, default=None

Number of jobs to run in parallel.
``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
``-1`` means using all processors. See :term:`Glossary `
for more details.

.. versionchanged:: v0.20
`n_jobs` default changed from 1 to None 		-1
refit refit: bool, str, or callable, default=True

Refit an estimator using the best found parameters on the whole
dataset.

For multiple metric evaluation, this needs to be a `str` denoting the
scorer that would be used to find the best parameters for refitting
the estimator at the end.

Where there are considerations other than maximum score in
choosing a best estimator, ``refit`` can be set to a function which
returns the selected ``best_index_`` given the ``cv_results_``. In that
case, the ``best_estimator_`` and ``best_params_`` will be set
according to the returned ``best_index_`` while the ``best_score_``
attribute will not be available.

The refitted estimator is made available at the ``best_estimator_``
attribute and permits using ``predict`` directly on this
``RandomizedSearchCV`` instance.

Also for multiple metric evaluation, the attributes ``best_index_``,
``best_score_`` and ``best_params_`` will only be available if
``refit`` is set and all of them will be determined w.r.t this specific
scorer.

See ``scoring`` parameter to know more about multiple metric
evaluation.

See :ref:`this example
`
for an example of how to use ``refit=callable`` to balance model
complexity and cross-validated score.

.. versionchanged:: 0.20
Support for callable added. 		True
cv cv: int, cross-validation generator or an iterable, default=None

Determines the cross-validation splitting strategy.
Possible inputs for cv are:

- None, to use the default 5-fold cross validation,
- integer, to specify the number of folds in a `(Stratified)KFold`,
- :term:`CV splitter`,
- An iterable yielding (train, test) splits as arrays of indices.

For integer/None inputs, if the estimator is a classifier and ``y`` is
either binary or multiclass, :class:`StratifiedKFold` is used. In all
other cases, :class:`KFold` is used. These splitters are instantiated
with `shuffle=False` so the splits will be the same across calls.

Refer :ref:`User Guide ` for the various
cross-validation strategies that can be used here.

.. versionchanged:: 0.22
``cv`` default value if None changed from 3-fold to 5-fold. 		StratifiedKFo... shuffle=True)
verbose verbose: int

Controls the verbosity: the higher, the more messages.

- >1 : the computation time for each fold and parameter candidate is
displayed;
- >2 : the score is also displayed;
- >3 : the fold and candidate parameter indexes are also displayed
together with the starting time of the computation. 		2
pre_dispatch pre_dispatch: int, or str, default='2*n_jobs'

Controls the number of jobs that get dispatched during parallel
execution. Reducing this number can be useful to avoid an
explosion of memory consumption when more jobs get dispatched
than CPUs can process. This parameter can be:

- None, in which case all the jobs are immediately created and spawned. Use
this for lightweight and fast-running jobs, to avoid delays due to on-demand
spawning of the jobs
- An int, giving the exact number of total jobs that are spawned
- A str, giving an expression as a function of n_jobs, as in '2*n_jobs' 		'2*n_jobs'
random_state random_state: int, RandomState instance or None, default=None

Pseudo random number generator state used for random uniform sampling
from lists of possible values instead of scipy.stats distributions.
Pass an int for reproducible output across multiple
function calls.
See :term:`Glossary `. 		42
error_score error_score: 'raise' or numeric, default=np.nan

Value to assign to the score if an error occurs in estimator fitting.
If set to 'raise', the error is raised. If a numeric value is given,
FitFailedWarning is raised. This parameter does not affect the refit
step, which will always raise the error. 		nan
return_train_score return_train_score: bool, default=False

If ``False``, the ``cv_results_`` attribute will not include training
scores.
Computing training scores is used to get insights on how different
parameter settings impact the overfitting/underfitting trade-off.
However computing the scores on the training set can be computationally
expensive and is not strictly required to select the parameters that
yield the best generalization performance.

.. versionadded:: 0.19

.. versionchanged:: 0.21
Default value was changed from ``True`` to ``False`` 		False
Parameters
transformers transformers: list of tuples

List of (name, transformer, columns) tuples specifying the
transformer objects to be applied to subsets of the data.

name : str
Like in Pipeline and FeatureUnion, this allows the transformer and
its parameters to be set using ``set_params`` and searched in grid
search.
transformer : {'drop', 'passthrough'} or estimator
Estimator must support :term:`fit` and :term:`transform`.
Special-cased strings 'drop' and 'passthrough' are accepted as
well, to indicate to drop the columns or to pass them through
untransformed, respectively.
columns : str, array-like of str, int, array-like of int, array-like of bool, slice or callable
Indexes the data on its second axis. Integers are interpreted as
positional columns, while strings can reference DataFrame columns
by name. A scalar string or int should be used where
``transformer`` expects X to be a 1d array-like (vector),
otherwise a 2d array will be passed to the transformer.
A callable is passed the input data `X` and can return any of the
above. To select multiple columns by name or dtype, you can use
:obj:`make_column_selector`. 		[('zv_filter', ...)]
remainder remainder: {'drop', 'passthrough'} or estimator, default='drop'

By default, only the specified columns in `transformers` are
transformed and combined in the output, and the non-specified
columns are dropped. (default of ``'drop'``).
By specifying ``remainder='passthrough'``, all remaining columns that
were not specified in `transformers`, but present in the data passed
to `fit` will be automatically passed through. This subset of columns
is concatenated with the output of the transformers. For dataframes,
extra columns not seen during `fit` will be excluded from the output
of `transform`.
By setting ``remainder`` to be an estimator, the remaining
non-specified columns will use the ``remainder`` estimator. The
estimator must support :term:`fit` and :term:`transform`.
Note that using this feature requires that the DataFrame columns
input at :term:`fit` and :term:`transform` have identical order. 		'passthrough'
sparse_threshold sparse_threshold: float, default=0.3

If the output of the different transformers contains sparse matrices,
these will be stacked as a sparse matrix if the overall density is
lower than this value. Use ``sparse_threshold=0`` to always return
dense. When the transformed output consists of all dense data, the
stacked result will be dense, and this keyword will be ignored. 		0.3
n_jobs n_jobs: int, default=None

Number of jobs to run in parallel.
``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
``-1`` means using all processors. See :term:`Glossary `
for more details. 		None
transformer_weights transformer_weights: dict, default=None

Multiplicative weights for features per transformer. The output of the
transformer is multiplied by these weights. Keys are transformer names,
values the weights. 		None
verbose verbose: bool, default=False

If True, the time elapsed while fitting each transformer will be
printed as it is completed. 		False
verbose_feature_names_out verbose_feature_names_out: bool, str or Callable[[str, str], str], default=True

- If True, :meth:`ColumnTransformer.get_feature_names_out` will prefix
all feature names with the name of the transformer that generated that
feature. It is equivalent to setting
`verbose_feature_names_out="{transformer_name}__{feature_name}"`.
- If False, :meth:`ColumnTransformer.get_feature_names_out` will not
prefix any feature names and will error if feature names are not
unique.
- If ``Callable[[str, str], str]``,
:meth:`ColumnTransformer.get_feature_names_out` will rename all the features
using the name of the transformer. The first argument of the callable is the
transformer name and the second argument is the feature name. The returned
string will be the new feature name.
- If ``str``, it must be a string ready for formatting. The given string will
be formatted using two field names: ``transformer_name`` and ``feature_name``.
e.g. ``"{feature_name}__{transformer_name}"``. See :meth:`str.format` method
from the standard library for more info.

.. versionadded:: 1.0

.. versionchanged:: 1.6
`verbose_feature_names_out` can be a callable or a string to be formatted. 		True
force_int_remainder_cols force_int_remainder_cols: bool, default=False

This parameter has no effect.

.. note::
If you do not access the list of columns for the remainder columns
in the `transformers_` fitted attribute, you do not need to set
this parameter.

.. versionadded:: 1.5

.. versionchanged:: 1.7
The default value for `force_int_remainder_cols` will change from
`True` to `False` in version 1.7.

.. deprecated:: 1.7
`force_int_remainder_cols` is deprecated and will be removed in 1.9. 		'deprecated' 
Index(['marriage', 'sex', 'education', 'LIMIT_BAL', 'age', 'pay_0', 'pay_2',
       'pay_3', 'pay_4', 'pay_5', 'pay_6', 'AVG_Bill_amt',
       'late_payment_count', 'average_bill_trend', 'average_payment',
       'avg_monthly_pay_to_bill', 'average_pay_to_bill_trend',
       'credit_utilization'],
      dtype='str')
Parameters
threshold threshold: float, default=0

Features with a training-set variance lower than this threshold will
be removed. The default is to keep all features with non-zero variance,
i.e. remove the features that have the same value in all samples. 		0.0001 
[]
passthrough
Parameters
objective objective: typing.Union[str, xgboost.sklearn._SklObjWProto, typing.Callable[[typing.Any, typing.Any], typing.Tuple[numpy.ndarray, numpy.ndarray]], NoneType]

Specify the learning task and the corresponding learning objective or a custom
objective function to be used.

For custom objective, see :doc:`/tutorials/custom_metric_obj` and
:ref:`custom-obj-metric` for more information, along with the end note for
function signatures. 		'binary:logistic'
base_score base_score: typing.Union[float, typing.List[float], NoneType]

The initial prediction score of all instances, global bias. 		None
booster None
callbacks callbacks: typing.Optional[typing.List[xgboost.callback.TrainingCallback]]

List of callback functions that are applied at end of each iteration.
It is possible to use predefined callbacks by using
:ref:`Callback API `.

.. note::

States in callback are not preserved during training, which means callback
objects can not be reused for multiple training sessions without
reinitialization or deepcopy.

.. code-block:: python

for params in parameters_grid:
# be sure to (re)initialize the callbacks before each run
callbacks = [xgb.callback.LearningRateScheduler(custom_rates)]
reg = xgboost.XGBRegressor(**params, callbacks=callbacks)
reg.fit(X, y) 		None
colsample_bylevel colsample_bylevel: typing.Optional[float]

Subsample ratio of columns for each level. 		None
colsample_bynode colsample_bynode: typing.Optional[float]

Subsample ratio of columns for each split. 		None
colsample_bytree colsample_bytree: typing.Optional[float]

Subsample ratio of columns when constructing each tree. 		0.8
device device: typing.Optional[str]

.. versionadded:: 2.0.0

Device ordinal, available options are `cpu`, `cuda`, and `gpu`. 		None
early_stopping_rounds early_stopping_rounds: typing.Optional[int]

.. versionadded:: 1.6.0

- Activates early stopping. Validation metric needs to improve at least once in
every **early_stopping_rounds** round(s) to continue training. Requires at
least one item in **eval_set** in :py:meth:`fit`.

- If early stopping occurs, the model will have two additional attributes:
:py:attr:`best_score` and :py:attr:`best_iteration`. These are used by the
:py:meth:`predict` and :py:meth:`apply` methods to determine the optimal
number of trees during inference. If users want to access the full model
(including trees built after early stopping), they can specify the
`iteration_range` in these inference methods. In addition, other utilities
like model plotting can also use the entire model.

- If you prefer to discard the trees after `best_iteration`, consider using the
callback function :py:class:`xgboost.callback.EarlyStopping`.

- If there's more than one item in **eval_set**, the last entry will be used for
early stopping. If there's more than one metric in **eval_metric**, the last
metric will be used for early stopping. 		None
enable_categorical enable_categorical: bool

See the same parameter of :py:class:`DMatrix` for details. 		False
eval_metric eval_metric: typing.Union[str, typing.List[typing.Union[str, typing.Callable]], typing.Callable, NoneType]

.. versionadded:: 1.6.0

Metric used for monitoring the training result and early stopping. It can be a
string or list of strings as names of predefined metric in XGBoost (See
:doc:`/parameter`), one of the metrics in :py:mod:`sklearn.metrics`, or any
other user defined metric that looks like `sklearn.metrics`.

If custom objective is also provided, then custom metric should implement the
corresponding reverse link function.

Unlike the `scoring` parameter commonly used in scikit-learn, when a callable
object is provided, it's assumed to be a cost function and by default XGBoost
will minimize the result during early stopping.

For advanced usage on Early stopping like directly choosing to maximize instead
of minimize, see :py:obj:`xgboost.callback.EarlyStopping`.

See :doc:`/tutorials/custom_metric_obj` and :ref:`custom-obj-metric` for more
information.

.. code-block:: python

from sklearn.datasets import load_diabetes
from sklearn.metrics import mean_absolute_error
X, y = load_diabetes(return_X_y=True)
reg = xgb.XGBRegressor(
tree_method="hist",
eval_metric=mean_absolute_error,
)
reg.fit(X, y, eval_set=[(X, y)]) 		'aucpr'
feature_types feature_types: typing.Optional[typing.Sequence[str]]

.. versionadded:: 1.7.0

Used for specifying feature types without constructing a dataframe. See
the :py:class:`DMatrix` for details. 		None
feature_weights feature_weights: Optional[ArrayLike]

Weight for each feature, defines the probability of each feature being selected
when colsample is being used. All values must be greater than 0, otherwise a
`ValueError` is thrown. 		None
gamma gamma: typing.Optional[float]

(min_split_loss) Minimum loss reduction required to make a further partition on
a leaf node of the tree. 		0.1
grow_policy grow_policy: typing.Optional[str]

Tree growing policy.

- depthwise: Favors splitting at nodes closest to the node,
- lossguide: Favors splitting at nodes with highest loss change. 		None
importance_type None
interaction_constraints interaction_constraints: typing.Union[str, typing.List[typing.Tuple[str]], NoneType]

Constraints for interaction representing permitted interactions. The
constraints must be specified in the form of a nested list, e.g. ``[[0, 1], [2,
3, 4]]``, where each inner list is a group of indices of features that are
allowed to interact with each other. See :doc:`tutorial
` for more information 		None
learning_rate learning_rate: typing.Optional[float]

Boosting learning rate (xgb's "eta") 		0.1
max_bin max_bin: typing.Optional[int]

If using histogram-based algorithm, maximum number of bins per feature 		None
max_cat_threshold max_cat_threshold: typing.Optional[int]

.. versionadded:: 1.7.0

.. note:: This parameter is experimental

Maximum number of categories considered for each split. Used only by
partition-based splits for preventing over-fitting. Also, `enable_categorical`
needs to be set to have categorical feature support. See :doc:`Categorical Data
` and :ref:`cat-param` for details. 		None
max_cat_to_onehot max_cat_to_onehot: Optional[int]

.. versionadded:: 1.6.0

.. note:: This parameter is experimental

A threshold for deciding whether XGBoost should use one-hot encoding based split
for categorical data. When number of categories is lesser than the threshold
then one-hot encoding is chosen, otherwise the categories will be partitioned
into children nodes. Also, `enable_categorical` needs to be set to have
categorical feature support. See :doc:`Categorical Data
` and :ref:`cat-param` for details. 		None
max_delta_step max_delta_step: typing.Optional[float]

Maximum delta step we allow each tree's weight estimation to be. 		None
max_depth max_depth: typing.Optional[int]

Maximum tree depth for base learners. 		7
max_leaves max_leaves: typing.Optional[int]

Maximum number of leaves; 0 indicates no limit. 		None
min_child_weight min_child_weight: typing.Optional[float]

Minimum sum of instance weight(hessian) needed in a child. 		3
missing missing: float

Value in the data which needs to be present as a missing value. Default to
:py:data:`numpy.nan`. 		nan
monotone_constraints monotone_constraints: typing.Union[typing.Dict[str, int], str, NoneType]

Constraint of variable monotonicity. See :doc:`tutorial `
for more information. 		None
multi_strategy multi_strategy: typing.Optional[str]

.. versionadded:: 2.0.0

.. note:: This parameter is working-in-progress.

The strategy used for training multi-target models, including multi-target
regression and multi-class classification. See :doc:`/tutorials/multioutput` for
more information.

- ``one_output_per_tree``: One model for each target.
- ``multi_output_tree``: Use multi-target trees. 		None
n_estimators n_estimators: Optional[int]

Number of boosting rounds. 		500
n_jobs n_jobs: typing.Optional[int]

Number of parallel threads used to run xgboost. When used with other
Scikit-Learn algorithms like grid search, you may choose which algorithm to
parallelize and balance the threads. Creating thread contention will
significantly slow down both algorithms. 		None
num_parallel_tree None
random_state random_state: typing.Union[numpy.random.mtrand.RandomState, numpy.random._generator.Generator, int, NoneType]

Random number seed.

.. note::

Using gblinear booster with shotgun updater is nondeterministic as
it uses Hogwild algorithm. 		None
reg_alpha reg_alpha: typing.Optional[float]

L1 regularization term on weights (xgb's alpha). 		None
reg_lambda reg_lambda: typing.Optional[float]

L2 regularization term on weights (xgb's lambda). 		None
sampling_method sampling_method: typing.Optional[str]

Sampling method. Used only by the GPU version of ``hist`` tree method.

- ``uniform``: Select random training instances uniformly.
- ``gradient_based``: Select random training instances with higher probability
when the gradient and hessian are larger. (cf. CatBoost) 		None
scale_pos_weight scale_pos_weight: typing.Optional[float]

Balancing of positive and negative weights. 		None
subsample subsample: typing.Optional[float]

Subsample ratio of the training instance. 		None
tree_method tree_method: typing.Optional[str]

Specify which tree method to use. Default to auto. If this parameter is set to
default, XGBoost will choose the most conservative option available. It's
recommended to study this option from the parameters document :doc:`tree method
` 		None
validate_parameters validate_parameters: typing.Optional[bool]

Give warnings for unknown parameter. 		None
verbosity verbosity: typing.Optional[int]

The degree of verbosity. Valid values are 0 (silent) - 3 (debug). 		None
use_label_encoder False
In [40]:
# 7. Best parameters
print("Best parameters:", search.best_params_)
print("Best f1:", search.best_score_)

Best parameters: {'model__n_estimators': 500, 'model__min_child_weight': 3, 'model__max_depth': 7, 'model__learning_rate': 0.1, 'model__gamma': 0.1, 'model__colsample_bytree': 0.8}
Best f1: 0.9229417566536997
In [41]:
# 8. Final model
xgb_model = search.best_estimator_
In [42]:
# Predict on validation set

# Predict class labels
y_pred_xgb = xgb_model.predict(X_val)

# Predict probabilities (for ROC-AUC)
y_val_proba_xgb = xgb_model.predict_proba(X_val)[:, 1]  # probability for class 1
y_train_proba_xgb = xgb_model.predict_proba(X_balanced)[:, 1]
In [43]:
# Evaluate model

# Accuracy
accuracy = accuracy_score(y_true = y_val, y_pred = y_pred_xgb)
print("Accuracy:", accuracy)

# Classification report (precision, recall, F1-score)
print("Classification Report:\n", classification_report(y_true = y_val, y_pred = y_pred_xgb))

# ROC-AUC
auc_score_train = roc_auc_score(y_balanced, y_train_proba_xgb)
auc_score_val = roc_auc_score(y_val, y_val_proba_xgb)
print(f"AUC training: {auc_score_train:.4f}")
print(f"AUC validation: {auc_score_val:.4f}")

# Confusion matrix
cm = confusion_matrix(y_true = y_val, y_pred = y_pred_xgb, labels = xgb_model.classes_)
disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels= xgb_model.classes_)
disp.plot()
plt.show()

Accuracy: 0.7995198079231692
Classification Report:
               precision    recall  f1-score   support

           0       0.88      0.88      0.88      2023
           1       0.47      0.47      0.47       476

    accuracy                           0.80      2499
   macro avg       0.67      0.67      0.67      2499
weighted avg       0.80      0.80      0.80      2499

AUC training: 0.9990
AUC validation: 0.7658
In [44]:
# Feature importance
# 1. Extract fitted XGBClassifier
xgb_model_fitted = xgb_model.named_steps['model']

# 3. Feature importance
importance_dict = xgb_model_fitted.get_booster().get_score(importance_type='gain')

# 4. Map 'f0', 'f1', ... to real feature names
importance_df = pd.DataFrame({
    'Feature': features,
    'Importance': list(importance_dict.values())
}).sort_values(by='Importance', ascending=False)

print(importance_df)

# 5. Plot
plt.figure(figsize=(10,6))
sns.barplot(x='Importance', y='Feature', data=importance_df)
plt.title("XGBoost Feature Importance")
plt.xlabel("Importance (gain)")
plt.ylabel("Feature")
plt.tight_layout()
plt.show()

                      Feature  Importance
12         late_payment_count  100.131767
5                       pay_0   20.615599
6                       pay_2    6.115468
9                       pay_5    5.334560
17         credit_utilization    4.817808
8                       pay_4    4.780683
14            average_payment    4.742585
15    avg_monthly_pay_to_bill    4.687218
7                       pay_3    4.647171
3                   LIMIT_BAL    4.585955
2                   education    4.458379
13         average_bill_trend    4.420187
16  average_pay_to_bill_trend    4.359066
4                         age    4.285543
11               AVG_Bill_amt    4.244455
1                         sex    3.667030
0                    marriage    3.489526
10                      pay_6    3.223861

Random forest¶

In [45]:
rf = RandomForestClassifier(criterion='gini',
                            random_state=42,
                            class_weight='balanced')
In [46]:
# 2. Pipeline
pipeline = Pipeline(steps=[
    ("preprocess", trees_preprocessor), # Use the new specific preprocessor
    ("model", rf)
])
In [47]:
# 3. Hyperparameter
param_dist = {
    'model__n_estimators': randint(100, 300),
    'model__max_depth': randint(3, 15),
    'model__min_samples_split': randint(2, 10),
    'model__min_samples_leaf': randint(1, 5),
    'model__max_features': randint(1,6),
}
In [48]:
# 4. Cross-validation
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=9)
In [49]:
# 5. Hyperparameter tuning
random_search = RandomizedSearchCV(
    estimator=pipeline,
    param_distributions=param_dist,
    n_iter=20,
    scoring='average_precision',
    cv=cv,
    verbose=2,
    n_jobs=-1,
    random_state=42
)
In [50]:
# 6. Fit model
random_search.fit(X_balanced, y_balanced)

Fitting 10 folds for each of 20 candidates, totalling 200 fits
Out[50]:
RandomizedSearchCV(cv=StratifiedKFold(n_splits=10, random_state=9, shuffle=True),
                   estimator=Pipeline(steps=[('preprocess',
                                              ColumnTransformer(remainder='passthrough',
                                                                transformers=[('zv_filter',
                                                                               VarianceThreshold(threshold=0.0001),
                                                                               Index(['marriage', 'sex', 'education', 'LIMIT_BAL', 'age', 'pay_0', 'pay_2',
       'pay_3', 'pay_4', 'pay_5', 'pay_6', 'AVG_Bill_amt',
       'late_pay...
                                        'model__min_samples_leaf': <scipy.stats._distn_infrastructure.rv_discrete_frozen object at 0x000001F9838FCA90>,
                                        'model__min_samples_split': <scipy.stats._distn_infrastructure.rv_discrete_frozen object at 0x000001F9C8FFFCD0>,
                                        'model__n_estimators': <scipy.stats._distn_infrastructure.rv_discrete_frozen object at 0x000001F983F5CAD0>},
                   random_state=42, scoring='average_precision', verbose=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Parameters
estimator estimator: estimator object

An object of that type is instantiated for each grid point.
This is assumed to implement the scikit-learn estimator interface.
Either estimator needs to provide a ``score`` function,
or ``scoring`` must be passed. 		Pipeline(step...m_state=42))])
param_distributions param_distributions: dict or list of dicts

Dictionary with parameters names (`str`) as keys and distributions
or lists of parameters to try. Distributions must provide a ``rvs``
method for sampling (such as those from scipy.stats.distributions).
If a list is given, it is sampled uniformly.
If a list of dicts is given, first a dict is sampled uniformly, and
then a parameter is sampled using that dict as above. 		{'model__max_depth': <scipy.stats....001F9C93FAB90>, 'model__max_features': <scipy.stats....001F983D38FD0>, 'model__min_samples_leaf': <scipy.stats....001F9838FCA90>, 'model__min_samples_split': <scipy.stats....001F9C8FFFCD0>, ...}
n_iter n_iter: int, default=10

Number of parameter settings that are sampled. n_iter trades
off runtime vs quality of the solution. 		20
scoring scoring: str, callable, list, tuple or dict, default=None

Strategy to evaluate the performance of the cross-validated model on
the test set.

If `scoring` represents a single score, one can use:

- a single string (see :ref:`scoring_string_names`);
- a callable (see :ref:`scoring_callable`) that returns a single value;
- `None`, the `estimator`'s
:ref:`default evaluation criterion ` is used.

If `scoring` represents multiple scores, one can use:

- a list or tuple of unique strings;
- a callable returning a dictionary where the keys are the metric
names and the values are the metric scores;
- a dictionary with metric names as keys and callables as values.

See :ref:`multimetric_grid_search` for an example.

If None, the estimator's score method is used. 		'average_precision'
n_jobs n_jobs: int, default=None

Number of jobs to run in parallel.
``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
``-1`` means using all processors. See :term:`Glossary `
for more details.

.. versionchanged:: v0.20
`n_jobs` default changed from 1 to None 		-1
refit refit: bool, str, or callable, default=True

Refit an estimator using the best found parameters on the whole
dataset.

For multiple metric evaluation, this needs to be a `str` denoting the
scorer that would be used to find the best parameters for refitting
the estimator at the end.

Where there are considerations other than maximum score in
choosing a best estimator, ``refit`` can be set to a function which
returns the selected ``best_index_`` given the ``cv_results_``. In that
case, the ``best_estimator_`` and ``best_params_`` will be set
according to the returned ``best_index_`` while the ``best_score_``
attribute will not be available.

The refitted estimator is made available at the ``best_estimator_``
attribute and permits using ``predict`` directly on this
``RandomizedSearchCV`` instance.

Also for multiple metric evaluation, the attributes ``best_index_``,
``best_score_`` and ``best_params_`` will only be available if
``refit`` is set and all of them will be determined w.r.t this specific
scorer.

See ``scoring`` parameter to know more about multiple metric
evaluation.

See :ref:`this example
`
for an example of how to use ``refit=callable`` to balance model
complexity and cross-validated score.

.. versionchanged:: 0.20
Support for callable added. 		True
cv cv: int, cross-validation generator or an iterable, default=None

Determines the cross-validation splitting strategy.
Possible inputs for cv are:

- None, to use the default 5-fold cross validation,
- integer, to specify the number of folds in a `(Stratified)KFold`,
- :term:`CV splitter`,
- An iterable yielding (train, test) splits as arrays of indices.

For integer/None inputs, if the estimator is a classifier and ``y`` is
either binary or multiclass, :class:`StratifiedKFold` is used. In all
other cases, :class:`KFold` is used. These splitters are instantiated
with `shuffle=False` so the splits will be the same across calls.

Refer :ref:`User Guide ` for the various
cross-validation strategies that can be used here.

.. versionchanged:: 0.22
``cv`` default value if None changed from 3-fold to 5-fold. 		StratifiedKFo... shuffle=True)
verbose verbose: int

Controls the verbosity: the higher, the more messages.

- >1 : the computation time for each fold and parameter candidate is
displayed;
- >2 : the score is also displayed;
- >3 : the fold and candidate parameter indexes are also displayed
together with the starting time of the computation. 		2
pre_dispatch pre_dispatch: int, or str, default='2*n_jobs'

Controls the number of jobs that get dispatched during parallel
execution. Reducing this number can be useful to avoid an
explosion of memory consumption when more jobs get dispatched
than CPUs can process. This parameter can be:

- None, in which case all the jobs are immediately created and spawned. Use
this for lightweight and fast-running jobs, to avoid delays due to on-demand
spawning of the jobs
- An int, giving the exact number of total jobs that are spawned
- A str, giving an expression as a function of n_jobs, as in '2*n_jobs' 		'2*n_jobs'
random_state random_state: int, RandomState instance or None, default=None

Pseudo random number generator state used for random uniform sampling
from lists of possible values instead of scipy.stats distributions.
Pass an int for reproducible output across multiple
function calls.
See :term:`Glossary `. 		42
error_score error_score: 'raise' or numeric, default=np.nan

Value to assign to the score if an error occurs in estimator fitting.
If set to 'raise', the error is raised. If a numeric value is given,
FitFailedWarning is raised. This parameter does not affect the refit
step, which will always raise the error. 		nan
return_train_score return_train_score: bool, default=False

If ``False``, the ``cv_results_`` attribute will not include training
scores.
Computing training scores is used to get insights on how different
parameter settings impact the overfitting/underfitting trade-off.
However computing the scores on the training set can be computationally
expensive and is not strictly required to select the parameters that
yield the best generalization performance.

.. versionadded:: 0.19

.. versionchanged:: 0.21
Default value was changed from ``True`` to ``False`` 		False
Parameters
transformers transformers: list of tuples

List of (name, transformer, columns) tuples specifying the
transformer objects to be applied to subsets of the data.

name : str
Like in Pipeline and FeatureUnion, this allows the transformer and
its parameters to be set using ``set_params`` and searched in grid
search.
transformer : {'drop', 'passthrough'} or estimator
Estimator must support :term:`fit` and :term:`transform`.
Special-cased strings 'drop' and 'passthrough' are accepted as
well, to indicate to drop the columns or to pass them through
untransformed, respectively.
columns : str, array-like of str, int, array-like of int, array-like of bool, slice or callable
Indexes the data on its second axis. Integers are interpreted as
positional columns, while strings can reference DataFrame columns
by name. A scalar string or int should be used where
``transformer`` expects X to be a 1d array-like (vector),
otherwise a 2d array will be passed to the transformer.
A callable is passed the input data `X` and can return any of the
above. To select multiple columns by name or dtype, you can use
:obj:`make_column_selector`. 		[('zv_filter', ...)]
remainder remainder: {'drop', 'passthrough'} or estimator, default='drop'

By default, only the specified columns in `transformers` are
transformed and combined in the output, and the non-specified
columns are dropped. (default of ``'drop'``).
By specifying ``remainder='passthrough'``, all remaining columns that
were not specified in `transformers`, but present in the data passed
to `fit` will be automatically passed through. This subset of columns
is concatenated with the output of the transformers. For dataframes,
extra columns not seen during `fit` will be excluded from the output
of `transform`.
By setting ``remainder`` to be an estimator, the remaining
non-specified columns will use the ``remainder`` estimator. The
estimator must support :term:`fit` and :term:`transform`.
Note that using this feature requires that the DataFrame columns
input at :term:`fit` and :term:`transform` have identical order. 		'passthrough'
sparse_threshold sparse_threshold: float, default=0.3

If the output of the different transformers contains sparse matrices,
these will be stacked as a sparse matrix if the overall density is
lower than this value. Use ``sparse_threshold=0`` to always return
dense. When the transformed output consists of all dense data, the
stacked result will be dense, and this keyword will be ignored. 		0.3
n_jobs n_jobs: int, default=None

Number of jobs to run in parallel.
``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
``-1`` means using all processors. See :term:`Glossary `
for more details. 		None
transformer_weights transformer_weights: dict, default=None

Multiplicative weights for features per transformer. The output of the
transformer is multiplied by these weights. Keys are transformer names,
values the weights. 		None
verbose verbose: bool, default=False

If True, the time elapsed while fitting each transformer will be
printed as it is completed. 		False
verbose_feature_names_out verbose_feature_names_out: bool, str or Callable[[str, str], str], default=True

- If True, :meth:`ColumnTransformer.get_feature_names_out` will prefix
all feature names with the name of the transformer that generated that
feature. It is equivalent to setting
`verbose_feature_names_out="{transformer_name}__{feature_name}"`.
- If False, :meth:`ColumnTransformer.get_feature_names_out` will not
prefix any feature names and will error if feature names are not
unique.
- If ``Callable[[str, str], str]``,
:meth:`ColumnTransformer.get_feature_names_out` will rename all the features
using the name of the transformer. The first argument of the callable is the
transformer name and the second argument is the feature name. The returned
string will be the new feature name.
- If ``str``, it must be a string ready for formatting. The given string will
be formatted using two field names: ``transformer_name`` and ``feature_name``.
e.g. ``"{feature_name}__{transformer_name}"``. See :meth:`str.format` method
from the standard library for more info.

.. versionadded:: 1.0

.. versionchanged:: 1.6
`verbose_feature_names_out` can be a callable or a string to be formatted. 		True
force_int_remainder_cols force_int_remainder_cols: bool, default=False

This parameter has no effect.

.. note::
If you do not access the list of columns for the remainder columns
in the `transformers_` fitted attribute, you do not need to set
this parameter.

.. versionadded:: 1.5

.. versionchanged:: 1.7
The default value for `force_int_remainder_cols` will change from
`True` to `False` in version 1.7.

.. deprecated:: 1.7
`force_int_remainder_cols` is deprecated and will be removed in 1.9. 		'deprecated' 
Index(['marriage', 'sex', 'education', 'LIMIT_BAL', 'age', 'pay_0', 'pay_2',
       'pay_3', 'pay_4', 'pay_5', 'pay_6', 'AVG_Bill_amt',
       'late_payment_count', 'average_bill_trend', 'average_payment',
       'avg_monthly_pay_to_bill', 'average_pay_to_bill_trend',
       'credit_utilization'],
      dtype='str')
Parameters
threshold threshold: float, default=0

Features with a training-set variance lower than this threshold will
be removed. The default is to keep all features with non-zero variance,
i.e. remove the features that have the same value in all samples. 		0.0001 
[]
passthrough
Parameters
n_estimators n_estimators: int, default=100

The number of trees in the forest.

.. versionchanged:: 0.22
The default value of ``n_estimators`` changed from 10 to 100
in 0.22. 		157
criterion criterion: {"gini", "entropy", "log_loss"}, default="gini"

The function to measure the quality of a split. Supported criteria are
"gini" for the Gini impurity and "log_loss" and "entropy" both for the
Shannon information gain, see :ref:`tree_mathematical_formulation`.
Note: This parameter is tree-specific. 		'gini'
max_depth max_depth: int, default=None

The maximum depth of the tree. If None, then nodes are expanded until
all leaves are pure or until all leaves contain less than
min_samples_split samples. 		14
min_samples_split min_samples_split: int or float, default=2

The minimum number of samples required to split an internal node:

- If int, then consider `min_samples_split` as the minimum number.
- If float, then `min_samples_split` is a fraction and
`ceil(min_samples_split * n_samples)` are the minimum
number of samples for each split.

.. versionchanged:: 0.18
Added float values for fractions. 		5
min_samples_leaf min_samples_leaf: int or float, default=1

The minimum number of samples required to be at a leaf node.
A split point at any depth will only be considered if it leaves at
least ``min_samples_leaf`` training samples in each of the left and
right branches. This may have the effect of smoothing the model,
especially in regression.

- If int, then consider `min_samples_leaf` as the minimum number.
- If float, then `min_samples_leaf` is a fraction and
`ceil(min_samples_leaf * n_samples)` are the minimum
number of samples for each node.

.. versionchanged:: 0.18
Added float values for fractions. 		1
min_weight_fraction_leaf min_weight_fraction_leaf: float, default=0.0

The minimum weighted fraction of the sum total of weights (of all
the input samples) required to be at a leaf node. Samples have
equal weight when sample_weight is not provided. 		0.0
max_features max_features: {"sqrt", "log2", None}, int or float, default="sqrt"

The number of features to consider when looking for the best split:

- If int, then consider `max_features` features at each split.
- If float, then `max_features` is a fraction and
`max(1, int(max_features * n_features_in_))` features are considered at each
split.
- If "sqrt", then `max_features=sqrt(n_features)`.
- If "log2", then `max_features=log2(n_features)`.
- If None, then `max_features=n_features`.

.. versionchanged:: 1.1
The default of `max_features` changed from `"auto"` to `"sqrt"`.

Note: the search for a split does not stop until at least one
valid partition of the node samples is found, even if it requires to
effectively inspect more than ``max_features`` features. 		5
max_leaf_nodes max_leaf_nodes: int, default=None

Grow trees with ``max_leaf_nodes`` in best-first fashion.
Best nodes are defined as relative reduction in impurity.
If None then unlimited number of leaf nodes. 		None
min_impurity_decrease min_impurity_decrease: float, default=0.0

A node will be split if this split induces a decrease of the impurity
greater than or equal to this value.

The weighted impurity decrease equation is the following::

N_t / N * (impurity - N_t_R / N_t * right_impurity
- N_t_L / N_t * left_impurity)

where ``N`` is the total number of samples, ``N_t`` is the number of
samples at the current node, ``N_t_L`` is the number of samples in the
left child, and ``N_t_R`` is the number of samples in the right child.

``N``, ``N_t``, ``N_t_R`` and ``N_t_L`` all refer to the weighted sum,
if ``sample_weight`` is passed.

.. versionadded:: 0.19 		0.0
bootstrap bootstrap: bool, default=True

Whether bootstrap samples are used when building trees. If False, the
whole dataset is used to build each tree. 		True
oob_score oob_score: bool or callable, default=False

Whether to use out-of-bag samples to estimate the generalization score.
By default, :func:`~sklearn.metrics.accuracy_score` is used.
Provide a callable with signature `metric(y_true, y_pred)` to use a
custom metric. Only available if `bootstrap=True`.

For an illustration of out-of-bag (OOB) error estimation, see the example
:ref:`sphx_glr_auto_examples_ensemble_plot_ensemble_oob.py`. 		False
n_jobs n_jobs: int, default=None

The number of jobs to run in parallel. :meth:`fit`, :meth:`predict`,
:meth:`decision_path` and :meth:`apply` are all parallelized over the
trees. ``None`` means 1 unless in a :obj:`joblib.parallel_backend`
context. ``-1`` means using all processors. See :term:`Glossary
` for more details. 		None
random_state random_state: int, RandomState instance or None, default=None

Controls both the randomness of the bootstrapping of the samples used
when building trees (if ``bootstrap=True``) and the sampling of the
features to consider when looking for the best split at each node
(if ``max_features < n_features``).
See :term:`Glossary ` for details. 		42
verbose verbose: int, default=0

Controls the verbosity when fitting and predicting. 		0
warm_start warm_start: bool, default=False

When set to ``True``, reuse the solution of the previous call to fit
and add more estimators to the ensemble, otherwise, just fit a whole
new forest. See :term:`Glossary ` and
:ref:`tree_ensemble_warm_start` for details. 		False
class_weight class_weight: {"balanced", "balanced_subsample"}, dict or list of dicts, default=None

Weights associated with classes in the form ``{class_label: weight}``.
If not given, all classes are supposed to have weight one. For
multi-output problems, a list of dicts can be provided in the same
order as the columns of y.

Note that for multioutput (including multilabel) weights should be
defined for each class of every column in its own dict. For example,
for four-class multilabel classification weights should be
[{0: 1, 1: 1}, {0: 1, 1: 5}, {0: 1, 1: 1}, {0: 1, 1: 1}] instead of
[{1:1}, {2:5}, {3:1}, {4:1}].

The "balanced" mode uses the values of y to automatically adjust
weights inversely proportional to class frequencies in the input data
as ``n_samples / (n_classes * np.bincount(y))``

The "balanced_subsample" mode is the same as "balanced" except that
weights are computed based on the bootstrap sample for every tree
grown.

For multi-output, the weights of each column of y will be multiplied.

Note that these weights will be multiplied with sample_weight (passed
through the fit method) if sample_weight is specified. 		'balanced'
ccp_alpha ccp_alpha: non-negative float, default=0.0

Complexity parameter used for Minimal Cost-Complexity Pruning. The
subtree with the largest cost complexity that is smaller than
``ccp_alpha`` will be chosen. By default, no pruning is performed. See
:ref:`minimal_cost_complexity_pruning` for details. See
:ref:`sphx_glr_auto_examples_tree_plot_cost_complexity_pruning.py`
for an example of such pruning.

.. versionadded:: 0.22 		0.0
max_samples max_samples: int or float, default=None

If bootstrap is True, the number of samples to draw from X
to train each base estimator.

- If None (default), then draw `X.shape[0]` samples.
- If int, then draw `max_samples` samples.
- If float, then draw `max(round(n_samples * max_samples), 1)` samples. Thus,
`max_samples` should be in the interval `(0.0, 1.0]`.

.. versionadded:: 0.22 		None
monotonic_cst monotonic_cst: array-like of int of shape (n_features), default=None

Indicates the monotonicity constraint to enforce on each feature.
- 1: monotonic increase
- 0: no constraint
- -1: monotonic decrease

If monotonic_cst is None, no constraints are applied.

Monotonicity constraints are not supported for:
- multiclass classifications (i.e. when `n_classes > 2`),
- multioutput classifications (i.e. when `n_outputs_ > 1`),
- classifications trained on data with missing values.

The constraints hold over the probability of the positive class.

Read more in the :ref:`User Guide `.

.. versionadded:: 1.4 		None
In [51]:
# 7. Best parameters
print("Best Parameters:\n", random_search.best_params_)
print("Best average precision:", random_search.best_score_)

# 8. Final model
rf_model = random_search.best_estimator_

Best Parameters:
 {'model__max_depth': 14, 'model__max_features': 5, 'model__min_samples_leaf': 1, 'model__min_samples_split': 5, 'model__n_estimators': 157}
Best average precision: 0.9507978989011164
In [52]:
# Predict on validation set

# Predict class labels
y_pred_rf = rf_model.predict(X_val)

# Predict probabilities (for ROC-AUC)
y_val_proba_rf = rf_model.predict_proba(X_val)[:, 1]
y_train_proba_rf = rf_model.predict_proba(X_balanced)[:, 1]
In [53]:
# Accuracy
accuracy = accuracy_score(y_true = y_val, y_pred = y_pred_rf)
print("Accuracy:", accuracy)

# Classification report (precision, recall, F1-score)
print("Classification Report:\n", classification_report(y_true = y_val, y_pred = y_pred_rf))

# ROC-AUC
auc_score_train = roc_auc_score(y_balanced, y_train_proba_rf)
auc_score_val = roc_auc_score(y_val, y_val_proba_rf)
print(f"AUC training: {auc_score_train:.4f}")
print(f"AUC validation: {auc_score_val:.4f}")

# Confusion matrix
cm = confusion_matrix(y_val, y_pred_rf, labels=rf_model.classes_)
disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=rf_model.classes_)
disp.plot()
plt.show()

Accuracy: 0.8195278111244498
Classification Report:
               precision    recall  f1-score   support

           0       0.88      0.90      0.89      2023
           1       0.53      0.49      0.51       476

    accuracy                           0.82      2499
   macro avg       0.71      0.69      0.70      2499
weighted avg       0.81      0.82      0.82      2499

AUC training: 0.9876
AUC validation: 0.7872
In [54]:
# Feature importance
# Access the RandomForestClassifier model from the pipeline
rf_model_fitted = rf_model.named_steps['model']

importance_df = pd.DataFrame({
    "Feature": features,
    "Importance": rf_model_fitted.feature_importances_
}).sort_values(by="Importance", ascending=False)

importance_df

plt.figure(figsize=(20, 15))
sns.barplot(x='Importance', y='Feature', data=importance_df)
plt.title('Feature Importances from Random Forest')
plt.xlabel('Features', size = 20)
plt.ylabel('Importance', size = 20)
plt.xticks(rotation=45, ha='right', size = 20)
plt.tight_layout()
plt.show()

Probability calibration for the selected tree-based model¶

The neural-network section has been removed so that the notebook does not require TensorFlow. For credit scoring, the probability calibration is applied to the best-performing available tree-based model using a separate calibration set. This keeps calibration independent from model fitting and preserves the original validation and final test sets.

In [55]:
# Select the final tree-based model using validation ROC-AUC.
model_validation_auc = {
    "XGBoost": roc_auc_score(y_val, y_val_proba_xgb),
    "Random Forest": roc_auc_score(y_val, y_val_proba_rf)
}

print("Validation ROC-AUC:")
for name, score in model_validation_auc.items():
    print(f"{name}: {score:.4f}")

selected_model_name = max(model_validation_auc, key=model_validation_auc.get)
selected_model = xgb_model if selected_model_name == "XGBoost" else rf_model

print(f"\nSelected final model: {selected_model_name}")

Validation ROC-AUC:
XGBoost: 0.7658
Random Forest: 0.7872

Selected final model: Random Forest
In [56]:
# Obtain probabilities for the separate calibration and final test sets.
y_cal_proba_final = selected_model.predict_proba(X_cal)[:, 1]
y_val_proba_final = selected_model.predict_proba(X_val)[:, 1]
y_test_proba_final = selected_model.predict_proba(X_test)[:, 1]

Logistic probability calibration (Platt-style scaling)¶

In [57]:
# Fit the calibrator only on the separate calibration set.
# The calibration model learns how the model's raw probability scores
# correspond to the observed default frequency.
eps = 1e-6

def probability_to_logit(p):
    p = np.clip(p, eps, 1 - eps)
    return np.log(p / (1 - p))

calibration_model = LogisticRegression(solver="lbfgs")
calibration_model.fit(
    probability_to_logit(y_cal_proba_final).reshape(-1, 1),
    y_cal
)

# Apply calibration to validation and final test probabilities.
y_val_proba_calibrated = calibration_model.predict_proba(
    probability_to_logit(y_val_proba_final).reshape(-1, 1)
)[:, 1]

y_test_proba_calibrated = calibration_model.predict_proba(
    probability_to_logit(y_test_proba_final).reshape(-1, 1)
)[:, 1]

y_pred_final = (y_val_proba_calibrated >= 0.5).astype(int)
In [58]:
y_val_proba_calibrated
Out[58]:
array([0.16882286, 0.12793105, 0.05952449, ..., 0.13220038, 0.11071356,
       0.08740416], shape=(2499,))
In [59]:
# Evaluate the selected model using calibrated probabilities.
print("Selected model:", selected_model_name)
print("Accuracy:", accuracy_score(y_val, y_pred_final))
print("Classification Report:\n", classification_report(y_val, y_pred_final))
print(f"ROC-AUC validation: {roc_auc_score(y_val, y_val_proba_final):.4f}")
print(f"Brier score before calibration: {brier_score_loss(y_val, y_val_proba_final):.4f}")
print(f"Brier score after calibration:  {brier_score_loss(y_val, y_val_proba_calibrated):.4f}")

Selected model: Random Forest
Accuracy: 0.8407362945178071
Classification Report:
               precision    recall  f1-score   support

           0       0.85      0.97      0.91      2023
           1       0.69      0.29      0.41       476

    accuracy                           0.84      2499
   macro avg       0.77      0.63      0.66      2499
weighted avg       0.82      0.84      0.81      2499

ROC-AUC validation: 0.7872
Brier score before calibration: 0.1421
Brier score after calibration:  0.1225

Calibration curve¶

Calibration asks whether predicted probabilities correspond to observed default frequencies. For example, among customers assigned a PD of about 20%, a well-calibrated model should observe a default rate of approximately 20%.

In [63]:
# Compute calibration curves on the untouched final test set.
prob_true_raw, prob_pred_raw = calibration_curve(
    y_test,
    y_test_proba_final,
    n_bins=10,
    strategy="quantile"
)

prob_true_cal, prob_pred_cal = calibration_curve(
    y_test,
    y_test_proba_calibrated,
    n_bins=10,
    strategy="quantile"
)
In [64]:
# Visualize calibration before and after calibration.
plt.figure(figsize=(8, 6))
plt.plot([0, 1], [0, 1], "--", label="Perfect calibration")
plt.plot(
    prob_pred_raw,
    prob_true_raw,
    marker="o",
    label=f"Uncalibrated {selected_model_name}"
)
plt.plot(
    prob_pred_cal,
    prob_true_cal,
    marker="o",
    label=f"Calibrated {selected_model_name}"
)

plt.xlabel("Mean predicted probability")
plt.ylabel("Observed default rate")
plt.title("Calibration Curve")
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

print(f"Test Brier score before calibration: {brier_score_loss(y_test, y_test_proba_final):.4f}")
print(f"Test Brier score after calibration:  {brier_score_loss(y_test, y_test_proba_calibrated):.4f}")

Test Brier score before calibration: 0.1426
Test Brier score after calibration:  0.1208
In [65]:
# Compare average predicted PD with the observed default rate on the test set.
print(f"Mean predicted PD before calibration: {y_test_proba_final.mean():.4f}")
print(f"Mean predicted PD after calibration:  {y_test_proba_calibrated.mean():.4f}")
print(f"Observed default rate:               {y_test.mean():.4f}")

Mean predicted PD before calibration: 0.3273
Mean predicted PD after calibration:  0.1904
Observed default rate:               0.1903

The selected tree-based model is used as the final model. Its probability outputs are calibrated on a separate calibration set before being interpreted as Probability of Default (PD). The preprocessing/model pipeline and calibration model are saved for later integration.

In [66]:
# Save the selected model and calibration model.
joblib.dump(selected_model, "selected_credit_model.pkl")
Out[66]:
['selected_credit_model.pkl']
In [ ]:
# Save the probability calibrator.
joblib.dump(calibration_model, "pd_calibrator.pkl")
Out[ ]:
['pd_threshold.pkl']
In [96]:
# Save thresholdes
joblib.dump(pd_threshold, "pd_threshold.pkl")
joblib.dump(ead_threshold, "ead_threshold.pkl")
Out[96]:
['ead_threshold.pkl']

6. Expected loss and Risk management¶

6.1. Expected loss function¶

From the selected model, the distribution of expected loss will be estimated. We assume that in the event of default, customers are able to repay only 20% of the outstanding amount, implying a Loss Given Default (LGD) of 0.8.

The expected next-month bill is calculated as the product of the bill trend and the previous month’s bill. Exposure at Default (EAD) is defined as the minimum of the predicted next-month bill and the available credit limit. The Probability of Default (PD) used in expected loss is the calibrated output of the selected tree-based model.

In [68]:
# Calculate exposure at default (EAD)
df_train_clean["next_bill"] = (1+df_train_clean["average_bill_trend"]) * df_train_clean["Bill_amt6"]
df_train_clean["EAD"] = np.where(df_train_clean["next_bill"] > df_train_clean["LIMIT_BAL"], df_train_clean["LIMIT_BAL"], df_train_clean["next_bill"])
In [69]:
# Calculate calibrated probability of default (PD) for the portfolio.
# The selected model can directly transform the original feature matrix.
raw_pd_all = selected_model.predict_proba(X)[:, 1]

df_train_clean["PD"] = calibration_model.predict_proba(
    probability_to_logit(raw_pd_all).reshape(-1, 1)
)[:, 1]
In [70]:
df_train_clean["PD"]
Out[70]:
0        0.411320
1        0.096725
2        0.069068
3        0.060673
4        0.071265
           ...   
25242    0.201017
25243    0.171637
25244    0.218708
25245    0.053410
25246    0.158892
Name: PD, Length: 25247, dtype: float64
In [71]:
# Calculate expected loss (EL)
LGD = 0.8
df_train_clean["EL"] = df_train_clean["EAD"] * df_train_clean["PD"] * LGD
df_train_clean["EL"]
Out[71]:
0         5672.858535
1        22440.172502
2         2987.399169
3         3611.927799
4            0.267667
             ...     
25242    14473.240208
25243     2351.856315
25244     4307.178822
25245     8532.659665
25246      183.437534
Name: EL, Length: 25247, dtype: float64
In [72]:
p90_EL = df_train_clean["EL"].quantile(0.90)
p95_EL = df_train_clean["EL"].quantile(0.95)
p99_EL = df_train_clean["EL"].quantile(0.99)

plt.figure(figsize=(10,6))

plt.hist(df_train_clean["EL"], bins=50)

plt.axvline(
    p90_EL,
    linestyle='--',
    linewidth=2,
    label=f'90th percentile = {p90_EL:,.0f}'
)

plt.axvline(
    p95_EL,
    linestyle='-',
    linewidth=2,
    label=f'95th percentile = {p95_EL:,.0f}'
)

plt.axvline(
    p99_EL,
    linestyle=':',
    linewidth=2,
    label=f'99th percentile = {p99_EL:,.0f}'
)

plt.xlabel("Expected Loss")
plt.ylabel("Frequency")
plt.title("Expected Loss Distribution")
plt.legend()

plt.show()

The expected loss distribution is typically right-skewed because a small number of customers can combine high exposure with a high probability of default. The percentile values below are used to quantify tail risk and to define exposure-management scenarios.

6.2. Risk management strategy¶

Threshold and Probability of Default¶
In [73]:
# Check distribution of PD
p80_PD = df_train_clean["PD"].quantile(0.80)

plt.figure(figsize=(10,6))

plt.hist(df_train_clean["PD"], bins=50)

plt.axvline(
    p80_PD,
    linestyle='--',
    linewidth=2,
    label=f'80th percentile = {p80_PD:,.2f}'
)

plt.xlabel("Probability of Default (PD)")
plt.ylabel("Frequency")
plt.title("Probability of Default Distribution")
plt.legend()

plt.show()
In [74]:
# Evaluate the 80th-percentile PD threshold on the validation set.
p80_PD_val = np.quantile(y_val_proba_calibrated, 0.80)
y_val_p80 = (y_val_proba_calibrated >= p80_PD_val).astype(int)

detected = np.sum((y_val == 1) & (y_val_p80 == 1))
total = np.sum(y_val == 1)

print(f"80th-percentile PD threshold: {p80_PD_val:.3f}")
print(f"Defaults detected on validation set: {detected}/{total} ({detected/total:.2%})")

80th-percentile PD threshold: 0.251
Defaults detected on validation set: 249/476 (52.31%)

Thresholds should be selected using the validation set rather than the full dataset. This avoids using the same observations for threshold selection and performance reporting. The final test set is reserved for an unbiased assessment of the chosen threshold.

In [97]:
# Select the PD cutoff on the validation set.
threshold_grid = np.arange(0.05, 0.96, 0.01)
threshold_results = []

for threshold in threshold_grid:
    y_pred_threshold = (y_val_proba_calibrated >= threshold).astype(int)
    threshold_results.append({
        "Threshold": threshold,
        "F1": f1_score(y_val, y_pred_threshold),
        "Recall": recall_score(y_val, y_pred_threshold)
    })

threshold_results = pd.DataFrame(threshold_results)
print(threshold_results)

# F1 is the primary criterion; among thresholds close to the best F1,
# prefer the one with higher recall for the Default class.
best_f1 = threshold_results["F1"].max()
candidates = threshold_results[threshold_results["F1"] >= best_f1 - 0.01]
pd_threshold = candidates.loc[candidates["Recall"].idxmax(), "Threshold"]

print(f"Selected PD threshold: {pd_threshold:.2f}")
print(f"Validation F1: {threshold_results.loc[threshold_results['Threshold'] == pd_threshold, 'F1'].iloc[0]:.3f}")
print(f"Validation recall: {threshold_results.loc[threshold_results['Threshold'] == pd_threshold, 'Recall'].iloc[0]:.3f}")

    Threshold        F1    Recall
0        0.05  0.340762  0.995798
1        0.06  0.352367  0.985294
2        0.07  0.360552  0.960084
3        0.08  0.377611  0.949580
4        0.09  0.396413  0.928571
..        ...       ...       ...
86       0.91  0.008368  0.004202
87       0.92  0.004193  0.002101
88       0.93  0.000000  0.000000
89       0.94  0.000000  0.000000
90       0.95  0.000000  0.000000

[91 rows x 3 columns]
Selected PD threshold: 0.23
Validation F1: 0.507
Validation recall: 0.559

The PD cutoff is selected from the validation set using F1-score as the primary criterion and recall as a secondary criterion. This is preferable to choosing a fixed cutoff such as 0.45 without checking whether the predicted probabilities are calibrated. The selected cutoff is then evaluated once on the untouched test set.

In [76]:
# Final evaluation of the selected PD threshold on the untouched test set.
y_test_pred_threshold = (y_test_proba_calibrated >= pd_threshold).astype(int)
print("Test F1:", f1_score(y_test, y_test_pred_threshold))
print("Test recall:", recall_score(y_test, y_test_pred_threshold))
print("Test precision:", precision_score(y_test, y_test_pred_threshold))

Test F1: 0.5156468939747781
Test recall: 0.5714285714285714
Test precision: 0.4697872340425532
Loss Given Default¶
In [77]:
df_train_clean[df_train_clean["next_bill"] >= df_train_clean["LIMIT_BAL"]]["next_month_default"].sum()
Out[77]:
np.int64(496)

The number of customers reaching their credit limit is reported by the calculation above and helps assess how frequently the EAD constraint is binding.

In [84]:
# Check distribution of EAD
p90_EAD = df_train_clean["EAD"].quantile(0.90)
p95_EAD = df_train_clean["EAD"].quantile(0.95)
p99_EAD = df_train_clean["EAD"].quantile(0.99)
p50_EAD = df_train_clean["EAD"].quantile(0.50)

plt.figure(figsize=(10,6))

plt.hist(df_train_clean["EAD"], bins=50)

plt.axvline(
    p90_EAD,
    linestyle='--',
    linewidth=2,
    label=f'90th percentile = {p90_EAD:,.0f}'
)

plt.axvline(
    p95_EAD,
    linestyle='-',
    linewidth=2,
    label=f'95th percentile = {p95_EAD:,.0f}'
)

plt.axvline(
    p99_EAD,
    linestyle=':',
    linewidth=2,
    label=f'99th percentile = {p99_EAD:,.0f}'
)

plt.axvline(
    p50_EAD,
    linestyle='-.',
    linewidth=2,
    color='green',
    label=f'50th percentile = {p50_EAD:,.0f}'
)

plt.xlabel("Exposure at Default (EAD)")
plt.ylabel("Frequency")
plt.title("Exposure at Default (EAD) Distribution")
plt.legend()

plt.show()

Exposure at Default is typically right-skewed. The 95th percentile is used as an exposure threshold so that customers in the upper tail can be explicitly identified for risk-management purposes.

Strategy¶

In practice, since other components of expected loss cannot be directly controlled, we focus on managing credit limit as a lever to influence Exposure at Default (EAD).

The high-risk rule combines the data-driven calibrated PD cutoff selected on the validation set (pd_threshold) with the 95th-percentile EAD threshold. A customer is classified as high risk if the calibrated PD exceeds the selected cutoff or if EAD exceeds the exposure threshold.

In [85]:
# Classify the full portfolio for the risk-management simulation.
ead_threshold = df_train_clean["EAD"].quantile(0.95)

df_train_clean["class"] = np.where(
    (df_train_clean["PD"] >= pd_threshold) |
    (df_train_clean["EAD"] >= ead_threshold),
    "High Risk",
    "Low Risk"
)
In [80]:
# Share of customers classified as high risk
high_risk_share = (df_train_clean["class"] == "High Risk").mean()
print(f"High-risk share: {high_risk_share:.2%}")

High-risk share: 27.67%
In [81]:
# Number of customers in each risk group
df_train_clean["class"].value_counts()
Out[81]:
class
Low Risk     18261
High Risk     6986
Name: count, dtype: int64

The high-risk share is determined by the calibrated PD cutoff and the EAD threshold and should be reported from the calculation above rather than fixed in advance.

In [82]:
# Evaluate the classification rule on the untouched test set.
test_portfolio = df_train_clean.loc[X_test.index].copy()
test_portfolio["PD"] = y_test_proba_calibrated
test_portfolio["EAD"] = np.minimum(
    test_portfolio["next_bill"],
    test_portfolio["LIMIT_BAL"]
)

test_portfolio["class"] = np.where(
    (test_portfolio["PD"] >= pd_threshold) |
    (test_portfolio["EAD"] >= ead_threshold),
    "High Risk",
    "Low Risk"
)

default_capture_rate = (
    test_portfolio.loc[
        (test_portfolio["class"] == "High Risk") &
        (test_portfolio["next_month_default"] == 1)
    ].shape[0]
    / test_portfolio["next_month_default"].sum()
)

print(f"Default capture rate on test set: {default_capture_rate:.2%}")

Default capture rate on test set: 60.35%

The test-set default capture rate provides a more defensible estimate of how well the high-risk rule identifies customers who subsequently default. Credit-limit adjustments should still be treated as a risk-management intervention rather than an automatic consequence of model classification.

We then evaluate the impact of alternative EAD thresholds and exposure reductions on the simulated expected-loss distribution.

In [86]:
for i in np.arange(0.5, 1, 0.05):
    df_train_clean["new_LMIT"] = np.where(
        df_train_clean["class"] == "High Risk",
        df_train_clean["LIMIT_BAL"] * i,
        df_train_clean["LIMIT_BAL"]
    )
    df_train_clean["new_EAD"] = np.where(
        df_train_clean["next_bill"] > df_train_clean["new_LMIT"],
        df_train_clean["new_LMIT"],
        df_train_clean["next_bill"],
    )
    df_train_clean["new_EL"] = df_train_clean["new_EAD"] * df_train_clean["PD"] * LGD
    new_p95 = df_train_clean["new_EL"].quantile(0.95)
    new_p99 = df_train_clean["new_EL"].quantile(0.99)
    reduce_95 = (( new_p95 - p95_EL ) / p95_EL) * 100
    reduce_99 = (( new_p99 - p99_EL ) / p99_EL) * 100
    print(f"i = {i:.2f}, 95th percentile of Expected Loss changes: {reduce_95:.2f}%, 99th percentile of Expected Loss changes: {reduce_99:.2f}%")

i = 0.50, 95th percentile of Expected Loss changes: -28.71%, 99th percentile of Expected Loss changes: -34.08%
i = 0.55, 95th percentile of Expected Loss changes: -24.54%, 99th percentile of Expected Loss changes: -29.03%
i = 0.60, 95th percentile of Expected Loss changes: -20.40%, 99th percentile of Expected Loss changes: -23.96%
i = 0.65, 95th percentile of Expected Loss changes: -17.29%, 99th percentile of Expected Loss changes: -19.67%
i = 0.70, 95th percentile of Expected Loss changes: -14.45%, 99th percentile of Expected Loss changes: -17.22%
i = 0.75, 95th percentile of Expected Loss changes: -11.95%, 99th percentile of Expected Loss changes: -12.87%
i = 0.80, 95th percentile of Expected Loss changes: -8.88%, 99th percentile of Expected Loss changes: -9.04%
i = 0.85, 95th percentile of Expected Loss changes: -6.07%, 99th percentile of Expected Loss changes: -6.82%
i = 0.90, 95th percentile of Expected Loss changes: -3.06%, 99th percentile of Expected Loss changes: -3.95%
i = 0.95, 95th percentile of Expected Loss changes: -1.07%, 99th percentile of Expected Loss changes: -0.89%

The credit-limit reduction results should be interpreted as a scenario analysis. The recalculated 95th and 99th percentiles show how the assumed exposure reduction changes model-based expected loss, rather than proving that realized losses will fall by the same percentage.

In [89]:
i = 0.75
df_train_clean["new_LIMIT"] = np.where(
        df_train_clean["class"] == "High Risk",
        df_train_clean["LIMIT_BAL"] * i,
        df_train_clean["LIMIT_BAL"]
    )

df_train_clean["new_EAD"] = np.where(
    df_train_clean["next_bill"] > df_train_clean["new_LIMIT"],
    df_train_clean["new_LIMIT"],
    df_train_clean["next_bill"],
)
df_train_clean["new_EL"] = df_train_clean["new_EAD"] * df_train_clean["PD"] * LGD
new_p95 = df_train_clean["new_EL"].quantile(0.95)
new_p99 = df_train_clean["new_EL"].quantile(0.99)
print(f"i = {i:.2f}, 95th percentile of Expected Loss: {new_p95:.2f}, 99th percentile of Expected Loss: {new_p99:.2f}")
print(f"95th percentile of Expected Loss changes: {(( new_p95 - p95_EL ) / p95_EL) * 100:.2f}%, 99th percentile of Expected Loss changes: {(( new_p99 - p99_EL ) / p99_EL) * 100:.2f}%")

i = 0.75, 95th percentile of Expected Loss: 25664.26, 99th percentile of Expected Loss: 73928.49
95th percentile of Expected Loss changes: -11.95%, 99th percentile of Expected Loss changes: -12.87%

The 25% credit-limit reduction scenario should be interpreted by comparing the recalculated 95th and 99th percentiles of expected loss with the baseline. The percentage changes are scenario-based estimates and depend on the assumed LGD and calibrated PDs.

Another approach is to penalize customers with late payments by reducing their credit limits, as late payment behavior is one of the strongest indicators of future default. Under this strategy, the credit limit is reduced by 10% for each month in which a customer makes a late payment.

In [90]:
df_train_clean["late_payment_months"] = (
    df_train_clean[[f"pay_{i}" for i in [0, 2, 3, 4, 5, 6]]] >= 1
).sum(axis=1)
In [91]:
df_train_clean["new_limit_late"] = df_train_clean["LIMIT_BAL"] - df_train_clean["late_payment_months"] * 0.1 * df_train_clean["LIMIT_BAL"]
df_train_clean["new_limit_late"] = df_train_clean["new_limit_late"].clip(lower=0)  # Ensure limit doesn't go negative
In [92]:
df_train_clean["new_EAD_late"] = np.where(
    df_train_clean["next_bill"] > df_train_clean["new_limit_late"],
    df_train_clean["new_limit_late"],
    df_train_clean["next_bill"],
)
df_train_clean["new_EL2"] = df_train_clean["new_EAD_late"] * df_train_clean["PD"] * LGD
new_p95_late = df_train_clean["new_EL2"].quantile(0.95)
new_p99_late = df_train_clean["new_EL2"].quantile(0.99)
reduce_95_late = ((new_p95_late - p95_EL) / p95_EL) * 100
reduce_99_late = ((new_p99_late - p99_EL) / p99_EL) * 100
print(f"95th percentile of Expected Loss: {new_p95_late:.2f}, 99th percentile of Expected Loss: {new_p99_late:.2f}")
print(f"95th percentile of Expected Loss changes: {reduce_95_late:.2f}%, 99th percentile of Expected Loss changes: {reduce_99_late:.2f}%")

95th percentile of Expected Loss: 24851.63, 99th percentile of Expected Loss: 65591.43
95th percentile of Expected Loss changes: -14.74%, 99th percentile of Expected Loss changes: -22.70%

The loss-reduction results should be interpreted as a scenario analysis, not as evidence that reducing credit limits will necessarily reduce realized losses. The estimated impact depends on the PD calibration, LGD assumption, EAD definition, and behavioral response of customers.

Considering both approaches, gradual adjustments based on repeated late-payment behavior may be preferable because they are more directly linked to observed repayment behavior and can avoid relying solely on a model classification.

Finally, classification can be used to identify high-risk customers and generate warnings within the system for further monitoring, such as: “This customer is classified as High Risk. Consider reducing credit limit and applying stricter credit terms." However, actual credit limit reductions should primarily be based on repeated late payment behavior rather than solely on the model classification result.