<jupyter_start><jupyter_text>Home Credit Default Risk Prediction**<NAME>**This is a case study on Home Credit Default Risk Prediction using data from Home Credit Group and various machine learning algorithms. Techniques on feature selection and principal component analysis are also considered. In the end, XGBoost algorithm is selected. The prediction AUC score is 0.75.   Загрузка данных<jupyter_code>!pip install kaggle
from google.colab import files
files.upload()
! mkdir ~/.kaggle
! cp kaggle.json ~/.kaggle/
! chmod 600 ~/.kaggle/kaggle.json
! kaggle competitions download -c 'home-credit-default-risk'
! unzip application_train.csv.zip application_train.csv
! unzip application_test.csv.zip application_test.csv
! unzip bureau.csv.zip bureau.csv<jupyter_output>Archive:  application_train.csv.zip
  inflating: application_train.csv   
Archive:  application_test.csv.zip
  inflating: application_test.csv    
Archive:  bureau.csv.zip
  inflating: bureau.csv<jupyter_text>Импортируем библиотеки<jupyter_code>import numpy as np
import pandas as pd
#import matplotlib.pyplot as plt
import matplotlib, matplotlib.pyplot as plt

from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.feature_selection import RFECV

from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier

from sklearn.model_selection import cross_val_score, GridSearchCV

from sklearn.metrics import confusion_matrix, classification_report, roc_auc_score
import xgboost as xgb
from tqdm import tqdm
from collections import defaultdict<jupyter_output><empty_output><jupyter_text>Explore and process bureau data<jupyter_code>bureau = pd.read_csv('bureau.csv')<jupyter_output><empty_output><jupyter_text>Shape<jupyter_code>bureau.shape<jupyter_output><empty_output><jupyter_text>Categorical features in bureau<jupyter_code>cat_features = bureau.select_dtypes(include='object').columns.tolist()
cat_features<jupyter_output><empty_output><jupyter_text>Label encode the categorical features<jupyter_code>bureau_transform = bureau.copy()
encoded_features = defaultdict(dict)
for col in cat_features:
    label_encoder = LabelEncoder()
    uniq_labels = bureau_transform[col].unique()
    bureau_transform[col] = label_encoder.fit_transform(bureau_transform[col])
    encoded_features[col] = dict(zip(uniq_labels,label_encoder.transform(uniq_labels)))
encoded_features
#feature_encode = onehot_encoder.get_feature_names(cat_features)
#bureau_cat_encode.columns = feature_encode
#bureau_cat_encode<jupyter_output><empty_output><jupyter_text>Extract numerical features and combine into the encoded categorical data<jupyter_code>#num_features = ['SK_ID_CURR', 'DAYS_CREDIT', 'CREDIT_DAY_OVERDUE', 'DAYS_CREDIT_ENDDATE', 'CNT_CREDIT_PROLONG', 'AMT_CREDIT_SUM', 'AMT_CREDIT_SUM_OVERDUE', 'DAYS_CREDIT_UPDATE', 'DAYS_CREDIT']
#bureau_num = bureau[num_features] 
#bureau_proc = pd.concat([bureau_num, bureau_cat_encode], axis=1)
#print(bureau_proc.shape)
#print(bureau_num.shape)
#print(bureau_cat_encode.shape)<jupyter_output><empty_output><jupyter_text>Group sum and mean bureau data by loan ID<jupyter_code>bureau_proc_sum = bureau_transform.groupby(['SK_ID_CURR']).sum()
bureau_proc_mean = bureau_transform.groupby(['SK_ID_CURR']).mean()<jupyter_output><empty_output><jupyter_text>Rename columns<jupyter_code>new_collist_sum = []
for col in bureau_proc_sum.columns:
    if col == 'SK_ID_CURR':
        continue
    new_collist_sum.append(col + ' SUMBYID')
colname_dict = dict(zip(bureau_proc_sum.columns, new_collist_sum))
bureau_proc_sum = bureau_proc_sum.rename(columns=colname_dict)
new_collist_mean = []
for col in bureau_proc_mean.columns:
    if col == 'SK_ID_CURR':
        continue
    new_collist_mean.append(col + ' MEANBYID')
colname_dict_mean = dict(zip(bureau_proc_mean.columns, new_collist_mean))
bureau_proc_mean = bureau_proc_mean.rename(columns=colname_dict_mean)<jupyter_output><empty_output><jupyter_text>MERGE mean and sum bureau data<jupyter_code>bureau_proc_final = pd.merge(bureau_proc_sum, bureau_proc_mean, how='left', on='SK_ID_CURR')
bureau_proc_final.to_csv('bureau_proc_final.csv', index=True)
bureau_proc_final = pd.read_csv('bureau_proc_final.csv')<jupyter_output><empty_output><jupyter_text>Process application_train/test<jupyter_code>application_train = pd.read_csv('application_train.csv')
application_test = pd.read_csv('application_test.csv')<jupyter_output><empty_output><jupyter_text>Check if 'SK_ID_CURR' in the train data<jupyter_code>'SK_ID_CURR' in application_train.columns<jupyter_output><empty_output><jupyter_text>Merge bureau_proc_final to application data<jupyter_code>train_merged = pd.merge(application_train, bureau_proc_final, how='left', on='SK_ID_CURR')
test_merged = pd.merge(application_test, bureau_proc_final, how='left', on='SK_ID_CURR')
print(application_train.shape)
print(train_merged.shape)
print(application_test.shape)
print(test_merged.shape)<jupyter_output>(307511, 122)
(307511, 154)
(48744, 121)
(48744, 153)<jupyter_text>Drop ID columns for training<jupyter_code>train_merged.drop('SK_ID_CURR', axis=1, inplace=True)
test_merged.drop('SK_ID_CURR', axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>DROP SOME FEATURES<jupyter_code>features_to_drop = ['EXT_SOURCE_1', 'EXT_SOURCE_2', 'EXT_SOURCE_3', 'FLAG_MOBIL', 'FLAG_CONT_MOBILE']
for feature in train_merged.columns:
    if 'FLAG_DOCUMENT' in feature:
        features_to_drop.append(feature)
train_merged.drop(features_to_drop, axis = 1, inplace = True)
test_merged.drop(features_to_drop, axis = 1, inplace = True)<jupyter_output><empty_output><jupyter_text>ADD NEW FEATURES<jupyter_code>train_merged['PRICE/CREDIT'] = train_merged['AMT_GOODS_PRICE']/train_merged['AMT_CREDIT']
test_merged['PRICE/CREDIT'] = test_merged['AMT_GOODS_PRICE']/test_merged['AMT_CREDIT']<jupyter_output><empty_output><jupyter_text>Extract numerical and categorical features for further processing: impute, scaling<jupyter_code>num_features_1 = train_merged.select_dtypes(include='int64').columns.tolist()
num_features_1.remove('TARGET')
num_features_2 = train_merged.select_dtypes(include='float64').columns.tolist()
cat_features = train_merged.select_dtypes(include='object').columns.tolist()
num_features = num_features_1 + num_features_2
features = num_features + cat_features<jupyter_output><empty_output><jupyter_text>EDA<jupyter_code>import seaborn as sns
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Обработка данных Заполняем пропущенные значеня в вещественных фичах<jupyter_code>for feature in num_features:
    fill_value = train_merged[feature].mean()
    train_merged[feature].fillna(fill_value, inplace = True)
    test_merged[feature].fillna(fill_value, inplace = True)
np.sum(train_merged[num_features].isna().sum() != 0)<jupyter_output><empty_output><jupyter_text>Заполняем пропущенные значения в категориальных фичах и делаем лейбл энкодинг<jupyter_code>for feature in cat_features:
    train_merged[feature].fillna('missing', inplace = True)
    test_merged[feature].fillna('missing', inplace = True)
    label_encoder = LabelEncoder()
    uniq_labels = train_merged[feature].unique()
    train_merged[feature] = label_encoder.fit_transform(train_merged[feature])

    encoded_features[feature] = dict(zip(uniq_labels,label_encoder.transform(uniq_labels)))
    test_merged[feature] = label_encoder.transform(test_merged[feature])
train_merged[cat_features].isna().sum()<jupyter_output><empty_output><jupyter_text>Сохраняем данные и делим выделяем таргер<jupyter_code>#train_merged['AGE/EXPERIENCE'] = train_merged['DAYS_BIRTH']/train_merged['DAYS_EMPLOYED']
#test_merged['AGE/EXPERIENCE'] = test_merged['DAYS_BIRTH']/test_merged['DAYS_EMPLOYED']
train_merged.to_csv('train_merged.csv', index= False)
test_merged.to_csv('test_merged.csv', index = False)

X_train = train_merged.drop(['TARGET'], axis = 1)
X_test = test_merged.copy()
y_train = train_merged.TARGET.values
print('X_train shape:', X_train.shape)<jupyter_output>X_train shape: (307511, 128)<jupyter_text>TRAINING difficult Model<jupyter_code>from lightgbm import LGBMClassifier
from xgboost import XGBClassifier
difficult_model_2 = LGBMClassifier(random_state=42, eval_metric="auc", max_depth=4,learning_rate=0.277, gamma=0.382)

difficult_model_2.fit(X_train, y_train)
#difficult_model = LGBMClassifier(random_state=42, eval_metric="auc", max_depth=4,learning_rate=0.277, gamma=0.382)

difficult_model = XGBClassifier(objective="binary:logistic", random_state=42, eval_metric="auc", 
                              max_depth=4,learning_rate=0.277, gamma=0.382)

difficult_model.fit(X_train, y_train)
import sklearn.metrics as metrics
def plot_AUC(model,X_test,y_test):
    probs = model.predict_proba(X_test)
    preds = probs[:,1]
    fpr, tpr, threshold = metrics.roc_curve(y_test, preds)
    roc_auc = metrics.auc(fpr, tpr)

    plt.title('Receiver Operating Characteristic')
    plt.plot(fpr, tpr, 'b', label = 'AUC = %0.2f' % roc_auc)
    plt.legend(loc = 'lower right')
    plt.plot([0, 1], [0, 1],'r--')
    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.ylabel('True Positive Rate')
    plt.xlabel('False Positive Rate')
    plt.show()
    return roc_auc*2-1
plot_AUC(difficult_model, X_train, y_train)
plot_AUC(difficult_model_2, X_train, y_train)
dif_model_predicted = [p[1] for p in difficult_model.predict_proba(X_train)]
dif_model_2_predicted = [p[1] for p in difficult_model_2.predict_proba(X_train)]<jupyter_output><empty_output><jupyter_text>Train logreg<jupyter_code>from sklearn.linear_model import LogisticRegression
feature_imp = sorted([(f,i) for f,i in zip(X_train.columns, difficult_model.feature_importances_)], key = lambda x: x[1], reverse=True)
selected_features = [f[0] for f in feature_imp[:20]]
logreg_model = LogisticRegression(random_state=42, max_iter=10000)
logreg_model.fit(X_train[selected_features], y_train)
plot_AUC(logreg_model, X_train[selected_features], y_train)<jupyter_output><empty_output><jupyter_text>Train riskpy<jupyter_code>!pip install git+https://github.com/Falconwatch/riskpy.git@master
from riskpy.modeling.binning import Binner
from riskpy.graphs.graphs import roc

import warnings
warnings.filterwarnings('ignore')
data_for_riskpy = X_train.copy()
data_for_riskpy['target'] = y_train
binner  = Binner()
fitted_bins = binner.fit(data_for_riskpy, target = 'target')
data_for_riskpy_trainsformed = binner.transform(data_for_riskpy)
fitted_columns = [fb._name+'_woe' for fb in fitted_bins if len(fb._gaps)>1]
logreg_riskpy = LogisticRegression( max_iter =1000)
logreg_riskpy.fit(data_for_riskpy_trainsformed[fitted_columns], y_train)
riskpy_predictions = [p[1] for p in logreg_riskpy.predict_proba(data_for_riskpy_trainsformed[fitted_columns])]
roc(y_train, riskpy_predictions)<jupyter_output>Gini: 0.41082720436687925<jupyter_text>Оценка экономического эффекта<jupyter_code>data_analysis = X_train.copy()

logreg_predictions = [p[1] for p in logreg_model.predict_proba(X_train[selected_features])]

#data_analysis['logreg_pd'] = logreg_predictions
data_analysis['logreg_pd'] = riskpy_predictions
data_analysis['diff_pd'] = dif_model_predicted
data_analysis['diff_2_pd'] = dif_model_2_predicted
data_analysis['target'] = y_train<jupyter_output><empty_output><jupyter_text>Оценка по RWA<jupyter_code>from scipy.stats import norm
import math

def get_rwa(my_pd, my_lgd, my_ead, my_r=0.03):
    k = (norm.cdf((norm.ppf(my_pd) + np.sqrt(my_r) * norm.ppf(0.999))/(np.sqrt(1-my_r))) - my_pd) * my_lgd
    rwa = k * 12.5 * my_ead
    return rwa

def eval_model(dataset, pd_column = 'pd_ttc', default_column = 'default_12m', ead_column = 'balance_at_report_dt'):
    g = roc_auc_score(dataset[default_column], dataset[pd_column])*2-1
    rwa_arr = get_rwa(dataset[pd_column], 0.72, dataset[ead_column])
    #rwa_ead_arr = rwa_arr / dataset[ead_column]
    rwa = rwa_arr.sum()/dataset[ead_column].sum()
    return g, rwa

def model_enchancment_profit(eads, rwa_base, rwa_adv):
  return np.sum(eads) * (rwa_base - rwa_adv)
def rwa_two_models(base_pd, adv_pd):
  print("Базовая модель")
  base_model = eval_model(data_analysis, pd_column=base_pd, default_column='target', ead_column='AMT_CREDIT')
  base_gini = base_model[0]
  base_rwa_rub = base_model[1]
  base_rwa = data_analysis['AMT_CREDIT'].sum() * base_rwa_rub
  print('Gini: ', base_gini)
  print('RWA на рубль: ', base_rwa_rub)
  print('RWA всего: ',base_rwa)
  print('-------------------------------')

  print("Продвинутая модель")
  advanced_model = eval_model(data_analysis, pd_column=adv_pd, default_column='target', ead_column='AMT_CREDIT')
  adv_gini = advanced_model[0]
  adv_rwa_rub = advanced_model[1]
  adv_rwa = data_analysis['AMT_CREDIT'].sum() * adv_rwa_rub
  print('Gini: ', adv_gini)
  print('RWA на рубль: ', adv_rwa_rub)
  print('RWA всего: ', adv_rwa)
  print('-------------------------------')

  print("Экономический эффект от улучшения модели")
  print(model_enchancment_profit(data_analysis['AMT_CREDIT'], base_rwa_rub, adv_rwa_rub))  
  return {'base_gini':base_gini, 'base_rwa_rub':base_rwa_rub, 'base_rwa':base_rwa, 
          'adv_gini':adv_gini, 'adv_rwa_rub':adv_rwa_rub, 'adv_rwa':adv_rwa}
print("База - xgb, продвинутый  - lgbm")
dif_models_comp = rwa_two_models('diff_pd', 'diff_2_pd')
dif_models_comp
print("База - логрег, продвинутый  - xgb")
logreg_xgb_comp = rwa_two_models('logreg_pd', 'diff_pd')
logreg_xgb_comp<jupyter_output>База - логрег, продвинутый  - xgb
Базовая модель
Gini:  0.41082720436687925
RWA на рубль:  0.8424178044155498
RWA всего:  155179327425.76343
-------------------------------
Продвинутая модель
Gini:  0.5324742200590096
RWA на рубль:  0.8373997355922147
RWA всего:  154254963599.52454
-------------------------------
Экономический эффект от улучшения модели
924363826.238904<jupyter_text>Динамика<jupyter_code>k=0.3

def improve_logreg(k):
  selected_index = np.random.choice(data_analysis.index.values, int(data_analysis.shape[0]*k), replace=False)
  old_pds = data_analysis['logreg_pd'].copy()
  old_pds.loc[selected_index] = data_analysis.loc[selected_index, 'diff_pd']
  data_analysis['improved_pd'] = old_pds
  new_model = eval_model(data_analysis, pd_column='improved_pd', default_column='target', ead_column='AMT_CREDIT')
  return new_model
improvement_log = list()
for k in tqdm(range(101)):
  improvement_log.append(improve_logreg(k/100))
improvement_ginis = [l[0] for l in improvement_log]
improvement_rwa_rub = [l[1] for l in improvement_log]
sns.scatterplot(x = improvement_ginis, y = improvement_rwa_rub)
plt.xlabel('Gini')
plt.ylabel('RWA на рубль Exposure')
plt.plot(improvement_ginis)
plt.plot(improvement_rwa_rub)<jupyter_output><empty_output><jupyter_text>Оценка по доходности<jupyter_code>def get_profit(pds, defaults, margins, eads, lgds = 0.72, th=0.04):
    predicted_class = (pds>th).astype(int)
    tn, fp, fn, tp = confusion_matrix(defaults, predicted_class).ravel()
    profit = np.sum(tn*margins - fp*margins - fn*eads*lgds)/np.sum(eads)
    return profit
def analize_profits(model_pd):
  profits = list()
  for t in range(50):
    th = t/100
    profit = get_profit(pds = data_analysis[model_pd], 
                defaults = data_analysis['target'], 
                margins = data_analysis['AMT_CREDIT']*0.1, 
                eads = data_analysis['AMT_CREDIT'],
                th=th)
    profits.append(profit)
  return profits
profits_logreg = analize_profits('logreg_pd')
profits_xgb = analize_profits('diff_pd')

plt.figure(figsize=[8,6])
plt.plot(profits_xgb, label="Доходность продвинутой модели")
plt.plot(profits_logreg, label="Доходность базовой модели")
plt.xlabel('Порог отсечения (%)')
plt.ylabel('Доходность (рублей/кредит)')
plt.title("Сравнение доходностей по двум моделям")
plt.legend()
profit_difference = np.array(profits_xgb) - np.array(profits_logreg)

plt.figure(figsize=[8,6])
plt.plot(profit_difference, label="Доходность продвинутой модели")
plt.xlabel('Порог отсечения (%)')
plt.ylabel('Доходность (рублей)')
plt.title("Доходность от продвинутой модели")
plt.legend()
def profit_comparisson(base_profits, adv_profits, total_count):
  base_best_profit = np.max(base_profits)
  adv_best_profit = np.max(adv_profits)
  base_total_profit = base_best_profit * total_count
  adv_total_profit = adv_best_profit * total_count

  return {'base_best_profit':base_best_profit, 
          'base_total_profit':base_total_profit, 
          'adv_best_profit':adv_best_profit, 
          'adv_total_profit':adv_total_profit}

logreg_xgb_profits_comp = profit_comparisson(profits_logreg, profits_xgb, X_train.shape[0])
logreg_xgb_profits_comp<jupyter_output><empty_output><jupyter_text>Таблица сравнения логрега и бустинга<jupyter_code>comparisson = pd.DataFrame()
comparisson.loc['gini', 'Логрег']=logreg_xgb_comp['base_gini']
comparisson.loc['gini', 'Бустинг']=logreg_xgb_comp['adv_gini']

comparisson.loc['Максимальная средняя доходность с кредита', 'Логрег']=logreg_xgb_profits_comp['base_best_profit']
comparisson.loc['Максимальная средняя доходность с кредита', 'Бустинг']=logreg_xgb_profits_comp['adv_best_profit']

comparisson.loc['Максимальный доход с потока (млн. руб)', 'Логрег']=round(logreg_xgb_profits_comp['base_total_profit']/1000000, 0)
comparisson.loc['Максимальный доход с потока (млн. руб)', 'Бустинг']=round(logreg_xgb_profits_comp['adv_total_profit']/1000000, 0)

comparisson.loc['RWA на рубль EAD', 'Логрег']=logreg_xgb_comp['base_rwa_rub']
comparisson.loc['RWA на рубль EAD', 'Бустинг']=logreg_xgb_comp['adv_rwa_rub']

comparisson.loc['RWA (млн. руб)', 'Логрег']=round(logreg_xgb_comp['base_rwa']/1000000)
comparisson.loc['RWA (млн. руб)', 'Бустинг']=round(logreg_xgb_comp['adv_rwa']/1000000)
comparisson['Изменение'] = comparisson['Бустинг'] - comparisson['Логрег']
comparisson['Относительное изменение (%)'] = comparisson['Изменение'] / comparisson['Логрег'] *100
pd.set_option('display.float_format', lambda x: '%.3f' % x)
comparisson<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12190173283614827, 'fraction_numerical': 0.01990568107041018, 'mean_word_length': 7.373278236914601, 'pattern_counts': {'":': 0, '<': 113, '<?xml version=': 0, '>': 115, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '238758', 'license': 'MIT', 'n_tokens_mistral': 7192, 'n_tokens_neox': 6937, 'n_words': 1065, 'path': 'Examples/Economic_effect_estimation_vs_riskpy.ipynb', 'repo_name': 'Sterals/INTERPRETATOR', 'size': '18238'}