{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":50160,"databundleVersionId":7921029,"sourceType":"competition"}],"dockerImageVersionId":30684,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<img src=\"https://i.ibb.co/BP7MsRh/header.png\" width=\"800\"/> \n\n## <center> Home Credit - Credit Risk Model Stability </center>\nWe aim to develop **a probability of default (PD) model** that upholds feature stability over time, consistent with the competition's standards. Our model is trained on \"depth of 0\" data, a snapshot reflecting borrower information at the loan application stage. We use a single-module model design and can conceptually divide the training data into three sub-categories: user profile, days past due (DPD), and credit bureau (CB) features.\n\nFor our predictive model, we have chosen **eXtreme Gradient Boosting (XGBoost)**, utilizing its logistic regression framework, as it closely resembles the original GBM (Gradient Boosting Machine) proposed by J. Friedman and has proven its utility in credit risk modeling. To address the issue of class imbalance, we use **a weighted log-likelihood** as our custom loss function.\n\nThe **modeling workflow** presented in this notebook is as follows:\n\n1. Data import and aggregation\n1. Model training via hyperparameter search\n1. Refitting the model with best hyperparameters\n1. Gathering model insights (evolution of PD over time, SHAP)\n1. Inference on the test data and submission file generation\n\n<div class=\"alert alert-block alert-warning\">\n<p class=\"first admonition-title\" style=\"font-weight: bold;\">Note</p>\n<p>Given the considerable size of the dataset and the notable class imbalance, exploring the undersampling of the non-default class represents a promising area for future research. This approach could reduce computational costs and broaden the range of viable training strategies for our models.</p>\n</div>\n\nAuthored by: www.github.com/deburky","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\ndataPath = \"/kaggle/input/home-credit-credit-risk-model-stability/\"\n\nlabels = pd.read_parquet(\n    dataPath + \"parquet_files/train/train_base.parquet\"\n)\n\n# Days Past Due and User Profile\nfeatures_dpd = [\n    \"maxdpdfrom6mto36m_3546853P\",\n    \"maxdpdlast12m_727P\",\n    \"maxdpdlast24m_143P\",\n    \"maxdpdlast3m_392P\",\n    \"maxdpdlast6m_474P\",\n]\n\nfeatures_user_profile = [\n    # debt characteristics\n    \"numinstls_657L\",\n    \"numinstlsallpaid_934L\",\n    \"pctinstlsallpaidlat10d_839L\",\n    \"totalsettled_863A\",\n    \"totaldebt_9A\",\n    \"currdebt_22A\",\n    \"credamount_770A\",\n    \"mobilephncnt_593L\",\n    \"homephncnt_628L\",\n    \"numactivecreds_622L\",\n    # application history\n    \"applicationcnt_361L\",\n    \"applications30d_658L\",\n    \"applicationscnt_1086L\",\n    \"applicationscnt_464L\",\n    \"applicationscnt_629L\",\n    # Additional features\n    \"avgdbddpdlast24m_3658932P\",\n    \"amtinstpaidbefduel24m_4187115A\",\n    \"maxdbddpdlast1m_3658939P\",\n    \"maxdbddpdtollast12m_3658940P\",\n    \"maxdbddpdtollast6m_4187119P\",\n    \"numinstpaidlate1d_3546852L\",\n]\n\ndpd_data = pd.read_parquet(\n    dataPath + \"parquet_files/train/train_static_0_1.parquet\",\n    columns=features_dpd + features_user_profile + [\"case_id\"],\n)\n\ntrain_data = labels.merge(\n    dpd_data,\n    on=\"case_id\",\n)\n\n# Credit Bureau\nfeatures_cb = [\n    \"numberofqueries_373L\",\n    \"days120_123L\",\n    \"days180_256L\",\n    \"days30_165L\",\n    \"days90_310L\",\n    \"days360_512L\",\n]\n\ncb_data = pd.read_parquet(\n    dataPath + \"parquet_files/train/train_static_cb_0.parquet\", \n    columns=features_cb + [\"case_id\"]\n)\n\ntrain_data = train_data.merge(cb_data, on=\"case_id\")\nprint(train_data.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-12T17:50:29.275746Z","iopub.execute_input":"2024-04-12T17:50:29.276628Z","iopub.status.idle":"2024-04-12T17:50:30.143864Z","shell.execute_reply.started":"2024-04-12T17:50:29.276587Z","shell.execute_reply":"2024-04-12T17:50:30.141687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Tuning hyperparameters\nIn this step we search for the best XGBoost parameters that optimize our metric of interest, Gini score, which is weighted by the difference between in-fold and out-of-fold variation using Optuna framework.\n\nTo further improve the performance of the model and address class imbalance, we use a custom loss function, the weighted log-likelihood. This approach allows us to assign a specific parameter, alpha ($\\alpha$), value to the negative (non-defaulted) class, effectively giving us a lever to adjust the model's sensitivity to the class of interest.\n\nBelow are some papers that discuss the weighted log-likelihood approach:\n\n- [Fithian W. and T. Hastie. 2014. Finite-Sample Equivalence in Statistical Model for Presence-Only Data.](https://arxiv.org/abs/1207.6950)\n\n- [Wang et al. 2021. Imbalance-XGBoost: Leveraging Weighted and Focal Losses for Binary Label-Imbalanced Classification with XGBoost.](https://arxiv.org/abs/1908.01672)","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX, y = (\n    train_data[features_dpd + features_cb + features_user_profile],\n    train_data[\"target\"],\n)\n\nix_train, ix_test = train_test_split(\n    X.index, stratify=y, test_size=0.3, random_state=42\n)","metadata":{"execution":{"iopub.status.busy":"2024-04-12T17:48:41.009366Z","iopub.execute_input":"2024-04-12T17:48:41.009895Z","iopub.status.idle":"2024-04-12T17:48:42.576900Z","shell.execute_reply.started":"2024-04-12T17:48:41.009857Z","shell.execute_reply":"2024-04-12T17:48:42.575418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import optuna\nimport numpy as np\nimport xgboost as xgb\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.model_selection import StratifiedKFold\nimport time\n\n\n# Weighted log-likelihood\ndef get_weighted_log_loss(alpha=0.25):\n    def weighted_loss_grad_hess(y_true, y_pred):\n        weights = np.where(y_true == 0.0, alpha, 1)\n        preds = 1.0 / (1.0 + np.exp(-y_pred))\n        grad = preds - y_true\n        hess = preds * (1.0 - preds)\n        return grad * weights, hess * weights\n\n    return weighted_loss_grad_hess\n\n\nstart_time = time.time()\n\n\n# Objective function with k-fold cross-validation\ndef objective(trial):\n    kfold = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n\n    param = {\n        \"verbosity\": 0,\n        \"objective\": get_weighted_log_loss(\n            alpha=trial.suggest_float(\"alpha\", 0.01, 0.99)\n        ),\n        \"n_estimators\": trial.suggest_int(\"n_estimators\", 50, 300),\n        \"learning_rate\": trial.suggest_float(\"learning_rate\", 0.01, 0.3),\n        \"max_depth\": trial.suggest_int(\"max_depth\", 2, 10),\n        \"base_score\": trial.suggest_float(\"base_score\", 0.01, 0.5),\n        \"missing\": np.nan,\n        \"random_state\": 42,\n    }\n\n    gini_scores = []\n\n    for train_index, test_index in kfold.split(X, y):\n        X_train, X_test = X.iloc[train_index], X.iloc[test_index]\n        y_train, y_test = y.iloc[train_index], y.iloc[test_index]\n\n        clf = xgb.XGBClassifier(**param)\n        clf.fit(X_train, y_train)\n\n        # Train scores\n        predictions_train = clf.predict_proba(X_train)[:, 1]\n        auc_score_train = roc_auc_score(y_train, predictions_train)\n        gini_score_train = 2 * auc_score_train - 1\n        train_error = 1 - gini_score_train\n\n        # Test scores\n        predictions_test = clf.predict_proba(X_test)[:, 1]\n        auc_score_test = roc_auc_score(y_test, predictions_test)\n        gini_score_test = 2 * auc_score_test - 1\n        test_error = 1 - gini_score_test\n\n        weight = abs(train_error / test_error)\n        gini_scores.append(gini_score_test * weight)\n\n    return np.mean(gini_scores)\n\n\n# Create the study object and optimize it\ntpe_sampler = optuna.samplers.TPESampler(seed=42)\nstudy = optuna.create_study(direction=\"maximize\", sampler=tpe_sampler)\nstudy.optimize(objective, n_trials=5)\n\nprint(f\"Best trial: {study.best_trial.params}\")\nprint(f\"Best score: {study.best_trial.value}\")\nprint(f\"Time taken: {time.time() - start_time} seconds\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Refit model with the best hyperparameters\n\nIn this step, we refit our model based on the best hyperparameters identified with Optuna and calculate train and test Gini scores.","metadata":{}},{"cell_type":"code","source":"import xgboost as xgb\nfrom sklearn.metrics import roc_auc_score\n\n# Weighted log-likelihood\ndef get_weighted_log_loss(alpha=0.25):\n    def weighted_loss_grad_hess(y_true, y_pred):\n        weights = np.where(y_true == 0.0, alpha, 1)\n        preds = 1.0 / (1.0 + np.exp(-y_pred))\n        grad = preds - y_true\n        hess = preds * (1.0 - preds)\n        return grad * weights, hess * weights\n\n    return weighted_loss_grad_hess\n\nbest_params = study.best_trial.params\n\nxgb_model = xgb.XGBClassifier(**best_params)\n\nxgb_model.fit(X.loc[ix_train], y.loc[ix_train])\n\npredictions_train = xgb_model.predict_proba(X.loc[ix_train])[:, 1]\ngini_train = 2 * roc_auc_score(y[ix_train], predictions_train) - 1\n\npredictions_test = xgb_model.predict_proba(X.loc[ix_test])[:, 1]\ngini_test = 2 * roc_auc_score(y[ix_test], predictions_test) - 1\n\nprint(f\"Train Gini: {gini_train:.2%}\")\nprint(f\"Test Gini: {gini_test:.2%}\")","metadata":{"execution":{"iopub.status.busy":"2024-04-12T17:49:44.841393Z","iopub.execute_input":"2024-04-12T17:49:44.842285Z","iopub.status.idle":"2024-04-12T17:50:07.088517Z","shell.execute_reply.started":"2024-04-12T17:49:44.842244Z","shell.execute_reply":"2024-04-12T17:50:07.087052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model insights\n\nHere we take a look at model's predictions over time and visualize feature importances using SHAP's beeswarm plot. Additionally, we look at the weights of contributing feature groups through SHAP values.","metadata":{}},{"cell_type":"code","source":"from matplotlib import pyplot as plt\n%config InlineBackend.figure_format = 'retina'\n\ntrain_data[\"prediction\"] = xgb_model.predict_proba(\n    train_data[features_dpd + features_cb + features_user_profile]\n)[:, 1]\ntrain_data[\"target\"].groupby(train_data[\"date_decision\"]).mean().plot(\n    kind=\"line\",\n    title=\"Through-the-economic-cycle split by week of origination\",\n    color=\"#f187ff\",\n    linewidth=0.8,\n    label=\"DR\",\n    rot=45,\n)\ntrain_data[\"prediction\"].groupby(train_data[\"date_decision\"]).mean().plot(\n    kind=\"line\", label=\"PD\", linewidth=0.8, color=\"#53cef9\", rot=45\n)\n\nplt.axhline(\n    y=train_data[\"target\"].mean(),\n    color=\"r\",\n    linestyle=\"--\",\n    linewidth=0.5,\n    label=\"Average DR\",\n)\nplt.legend()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shap\nfrom matplotlib import pyplot as plt\n\nshap.initjs()\n\n# SHAP explainer\nexplainer = shap.TreeExplainer(xgb_model)\nshap_values = explainer.shap_values(X.loc[ix_test])\n\n# Beeswarm plot\n_ = plt.figure(dpi=100)\nshap.summary_plot(\n    shap_values,\n    X.loc[ix_test],\n    plot_type=\"dot\",\n    max_display=10,\n    show=False,\n    cmap=\"cool\"\n)\nplt.gcf().set_size_inches(10,6)\nplt.tight_layout() \nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SHAP weights\ndef calculate_shap_weight(shap_values, feature_names, features_group):\n    \"\"\"Calculate the SHAP values weight for a given features group.\"\"\"\n    feature_indices = [\n        feature_names.get_loc(feature) for feature in features_group\n    ]\n    shap_values_group = shap_values[:, feature_indices]\n    return np.sum(np.abs(shap_values_group))\n\n\nfeature_names = X.columns\n\nfeature_groups = {\n    \"user_profile\": features_user_profile,\n    \"cb\": features_cb,\n    \"dpd\": features_dpd,\n}\n\n# Calculate the total SHAP value sum for normalization\ntotal_shap_sum = sum(\n    calculate_shap_weight(shap_values, feature_names, group)\n    for group in feature_groups.values()\n)\n\n# Calculate weights for each feature block\nweights = {}\nfor group_name, features_group in feature_groups.items():\n    weight = (\n        calculate_shap_weight(shap_values, feature_names, features_group)\n        / total_shap_sum\n    )\n    weights[group_name] = weight\n    print(f\"{group_name} weight: {weight:.2%}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model stability\n\nHere we follow the logic described in the guidelines by calculating the Gini scores by week of origination, falling rate, and the stability metric.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.linear_model import LinearRegression\n\nweek_nums = np.array(\n    list(train_data['WEEK_NUM'].unique())\n).reshape(-1, 1)\n\nginis = []\n\nfor value in week_nums:\n    subset = train_data.query(f'WEEK_NUM == {value}')\n    prediction, target = subset['prediction'], subset['target']\n    gini_score = 2 * roc_auc_score(target, prediction) - 1\n    ginis.append(gini_score)\n    \nginis = np.array(ginis).reshape(-1, 1)\n\n# Step 1: Fit a linear regression\nlin_reg = LinearRegression(fit_intercept=True)\nlin_reg.fit(week_nums, ginis)\na = lin_reg.coef_[0][0]\n\n# Step 2: Calculate a falling rate\nfalling_rate = min(0, a)\n\n# Step 3: Calculate the standard deviation of residuals\npredictions = lin_reg.predict(week_nums)\nresiduals = ginis - predictions\nstd_residuals = np.std(residuals)\n\n# Step 4: Compute the stability metric\nmean_ginis = np.mean(ginis)\nstability_metric = mean_ginis + 88.0 * falling_rate - 0.5 * std_residuals\n\n# Output the final stability metric\nprint(f\"Stability metric: {stability_metric:.2f}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test data\n\nHere we make predictions on the test data not used in model training and prepare a submission file.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nlabels_test = pd.read_csv(dataPath + \"csv_files/test/test_base.csv\")\n\ndpd_data_test = pd.read_csv(\n    dataPath + \"csv_files/test/test_static_0_0.csv\",\n    usecols=features_dpd + features_user_profile + [\"case_id\"],\n)\n\ntest_data = labels_test.merge(\n    dpd_data_test,\n    on=\"case_id\",\n    how='left'\n)\n\ncb_data_test = pd.read_csv(\n    dataPath + \"csv_files/test/test_static_cb_0.csv\", \n    usecols=features_cb + [\"case_id\"]\n)\n\ntest_data = test_data.merge(cb_data_test, on=\"case_id\", how='left')\nprint(test_data.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Make inference on the test data\ncase_ids = test_data[\"case_id\"].to_numpy()\nscore_prediction = xgb_model.predict_proba(\n    test_data[features_dpd + features_cb + features_user_profile]\n)[:, 1].astype(float)\n\n# Combine submission file and score prediction\nsubmission = pd.DataFrame(\n    {\n        \"case_id\": case_ids,\n        \"score\": score_prediction\n    }\n).set_index('case_id')\n\ndisplay(submission.head(5))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Save csv to the current working directory\nsubmission.to_csv(\"./submission.csv\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}