信用卡评分模型开发

一.背景：

信用风险指的是交易对手未能履行约定合同中的义务造成经济损失的风险，即受信人不能履行还本付息的责任而使授信人的预期收益与实际收益发生偏离可能性，为金融风险的主要类型
借贷的评分卡，是一种以分数的形式来衡量风险几率的一种手段，一般来说分数越高，风险越小
信用风险计量体系包括主体评级模型和债项评级两个部分。主体评级包含以下四个内容：
1.申请者评级 2.行为评级模型3.催收评级模型4.欺诈评级模型

二. 数据来源：

本项目来源kaggle竞赛 Give Me Some Credit
2.1获取数据
观察现有的数据和指标

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
data = pd.read_csv(r'E:/平时自学相关资料/CreditScoreModel-master/cs-training.csv')
data.describe()

2.2数据预处理
2.2.1缺失值

null=data.isnull().sum().sort_values(ascending=False)
null_rate =data.isnull().sum().sort_values(ascending=False)/len(data)
pd.concat([null,null_rate],axis=1)

缺失值的处理：
①缺失值极多：缺失样本占比极高，直接舍弃，作为特征加入反而会引入噪声
②非连续特征值缺失适中：将NaN作为一个新类别，加入到类别特征中
③连续值缺失适中：如果缺失样本适中，考虑给定一个step，然后离散化，将NaN作为一个type加入属性类目中
④缺失值较少：利用填充方法进行处理，均值，众数填充。RandomForest模型拟合数据，进行填充。要求有一定相关性。否则会引入噪声
MonthIncome缺失19.8%和NumberofDependents 2.6%，对MonthIncome采用随机森林对数据插补，由于NumberofDependents缺失数据较少，可以直接删掉

from sklearn.ensemble  import  RandomForestRegressor
def set_missing(df):
    process_df = df.iloc[:,[6,1,2,3,4,5,7,8,9,10]]
    known = process_df[process_df.MonthlyIncome.notnull()].values
    unknown = process_df[process_df.MonthlyIncome.isnull()].values
    X = known[:, 1:]
    y = known[:, 0]
    rfr = RandomForestRegressor(random_state=0, n_estimators=200,max_depth=3,n_jobs=-1)
    rfr.fit(X,y)
    predicted = rfr.predict(unknown[:, 1:]).round(0)
    print(predicted)
    df.loc[(df.MonthlyIncome.isnull()), 'MonthlyIncome'] = predicted
    return df

data=set_missing(data)
data=data.dropna()
data = data.drop_duplicates()

2.2.2 异常值
age

plt.boxplot(x=data.age)
plt.xlabel('age')
plt.title('Age')
plt.show()

由箱线图可知，年龄小于等于0为异常值，需要删除

RevolvingUtilizationOfUnsecuredLines（可用额度比值）及DebtRatio（负债率）

fig =plt.figure(figsize=(6,6))
ax =fig.add_subplot()
ax.boxplot([data.RevolvingUtilizationOfUnsecuredLines,data.DebtRatio])
ax.set_xticklabels(['RevolvingUtilizationOfUnsecuredLines','DebtRatio'])
plt.title('UnsecuredLines and DebtRatio')
plt.show()

(NumberOfTime30-59DaysPastDueNotWorse）逾期30天数；（NumberOfTimes90DaysLate）逾期90天数；（NumberOfTime60-89DaysPastDueNotWorse）逾期60天数

fig =plt.figure(figsize=(6,6))
ax =fig.add_subplot()
ax.boxplot([data['NumberOfTime30-59DaysPastDueNotWorse'],data['NumberOfTimes90DaysLate'],data['NumberOfTime60-89DaysPastDueNotWorse']])
ax.set_xticklabels(['30-90','90','60-89'])
plt.title('PastDueNotWorse')
plt.grid()
plt.show()

由图可知，均存在两个异常点，需要删除

data=data[data['age']>0]
data = data[data['NumberOfTime30-59DaysPastDueNotWorse'] < 90]

2.2.3 单变量分析
好坏客户的比值

data.SeriousDlqin2yrs.value_counts(normalize=True).plot.bar()
print("坏客户的比例%s"%[data.SeriousDlqin2yrs.value_counts(normalize=True)[1]/data.SeriousDlqin2yrs.value_counts(normalize=True)[0]])
plt.show()

2.2.4多变量的分析

fig = plt.figure(figsize=(7,6))
ax=sns.heatmap(data.corr(),annot=True,fmt='.2f')
bottom,top =ax.get_ylim()
ax.set_ylim()
ax.set_ylim(bottom+0.5,top-0.5)
plt.show()

由上图看出，各变量之间的相关性很小，不存在共线性可能

三.特征工程

特征选择非常重要，好的特征能够构造出好的模型，通过Python来实现，这里采用信用卡评分模型常用的IV值筛选。具体的IV值及WOE计算方法后面解释。
3.1 特征分箱
特征分箱指的是将连续变量离散化或将多状态的离散变量合并成少状态。离散特征的增加和减少都容易，易于模型的快速迭代，离散后的特征对异常数据有很强的鲁棒性，能够减少未离散化之前异常值对模型的干扰，离散化后可以进行特征交叉。本文的模型为逻辑回归，逻辑回归为广义线性模型，表达能力受限；单变量离散化为N个后，每个变量有单独权重，相当于为模型引入非线性，提升模型表达能力，加大拟合，同时降低了过拟合风险
特征分箱常用的有以下几种方法：有监督的Best-KS，ChiMerge，无监督的包括等频，等距，聚类。根据数据特征采用不同的分箱方式。这里采用聚类，等距，代码如下

from sklearn.model_selection import train_test_split
Y=data['SeriousDlqin2yrs']
X=data.iloc[:,1:]
X_train,X_test,Y_train,Y_test = train_test_split(X,Y,test_size=0.3,random_state=1)
train = pd.concat([Y_train,X_train],axis=1)
test = pd.concat([Y_test,X_test],axis=1)
train.to_csv('TrainData.csv',index=False)
test.to_csv('TestData.csv',index=False)

将数据分为训练集和测试集，便于后续验证，数据分箱代码如下：

#聚类分箱
from scipy.stats import spearmanr
def best_cut(X,Y,n):
    r= 0
    cutx =[]
    while abs(r)<1:
        d1= pd.DataFrame({'X':X,'Y':Y})
        cut= pd.qcut(X,n)
        d2= d1.groupby(cut)
        r,q = spearmanr(d2.mean().X,d2.mean().Y)
        n = n-1
    cutx.append(float('-inf'))
    for i in range(1,n+1):
        qua= X.quantile(i/(n+1))
        cutx.append(qua)
    cutx.append(float('inf'))    
    return cut,cutx

train_data =pd.read_csv('TrainData.csv')
cut1,cutx1=best_cut(train_data.RevolvingUtilizationOfUnsecuredLines,train_data.SeriousDlqin2yrs,10)
cut2,cutx2= best_cut(train_data.age,train_data.SeriousDlqin2yrs,10)
cut4,cutx4=best_cut(train_data.DebtRatio,train_data.SeriousDlqin2yrs,10)
cut5,cutx5=best_cut(train_data.MonthlyIncome,train_data.SeriousDlqin2yrs,10)

将可用额度比值，年龄，负债率，月收入进行聚类分箱，其余进行连续变量的离散化

pinf =float('inf')
ninf= float('-inf')
cutx3 = [ninf, 0, 1, 3, 5, pinf]
cutx6 = [ninf, 1, 2, 3, 5, pinf]
cutx7 = [ninf, 0, 1, 3, 5, pinf]
cutx8 = [ninf, 0,1,2, 3, pinf]
cutx9 = [ninf, 0, 1, 3, pinf]
cutx10 = [ninf, 0, 1, 2, 3, 5, pinf]
cut3=pd.cut(train_data["NumberOfTime30-59DaysPastDueNotWorse"],cutx3,labels=False)
cut6=pd.cut(train_data["NumberOfOpenCreditLinesAndLoans"],cutx6,labels=False)
cut7=pd.cut(train_data["NumberOfTimes90DaysLate"],cutx7,labels=False)
cut8=pd.cut(train_data["NumberRealEstateLoansOrLines"],cutx8,labels=False)
cut9=pd.cut(train_data["NumberOfTime60-89DaysPastDueNotWorse"],cutx9,labels=False)
cut10=pd.cut(train_data["NumberOfDependents"],cutx10,labels=False)

3.2 计算WOE（weight of evidence）
WOE的全称是“Weight of Evidence”，及证据权重，WOE是对原始自变量的一种编码形式 ，要对一个变量进行WOE编码，需要首先对这个变量进行 分组处理（也叫离散化、分箱等）分组后，对于第i组，WOE的计算公式如下：

pyi：这组中响应客户（风险模型中，对应的是违约客户，总之，指的是模型中预测变量取值为“是”或者1的个体）占所有样本中所有响应客户的比例

pni： 这组中未响应客户占样本中所有未响应客户的比例
yi：这个组中响应客户的数量
ni：这个组中未响应客户的数量
yT：是样本中所有响应的数量
nT：样本中所有未响应的数量

从上述公式中看出，WOE表示的实际上是“ 当前分组中响应客户占所有响应客户的比例”和“ 当前分组中没有响应的客户占所有没有响应的客户的比例”的差异
WOE越大，这种差异越大，这个分组里样本响应可能性就越大，WOE越小，差异越小。计算WOE代码

def cal_woe(df,cut):
    bad_sum = df.SeriousDlqin2yrs.sum()
    good_sum = df.SeriousDlqin2yrs.count()-bad_sum
    group = df.groupby(cut)['SeriousDlqin2yrs']
    cut_bad =group.sum()
    cut_good =group.count()-cut_bad
    woe= np.log((cut_bad/bad_sum)/(cut_good/good_sum))
    return woe

3.3 IV的概念
IV全称是information Value,中文是信息价值，或者信息量，IV这一指标就是用来衡量自变量预测能力（衡量自变量对目标变量影响程度的指标），类似的指标还有信息增益，基尼系数等。
假设在一个分类问题中，目标变量的类别有两类：Y1，Y2。对于一个待预测的个体A，要判断A属于Y1还是Y2，需要一定的信息，假定信息总量是I,而这些信息蕴含在自变量C1,C2,C3…,Cn中，那么对于其中的一个变量Ci来说，其蕴含的信息越多，那么它对于判断A属于Y1还是Y2的贡献就越大，Ci的信息价值就越大，Ci的IV就越大，它就越应该进入到入模变量

3.4 IV的计算

有了一个变量各分组的iv值，我们就可以计算整个变量的iv值，方法很简单，就是把各分组的IV相加：

其中，n为变量分组个数（即分箱后有多少组)

def get_iv(df,cut,cutwoe):
    bad_attribute = df.groupby(cut).sum().SeriousDlqin2yrs/df.SeriousDlqin2yrs.sum() 
    good_attribute =(df.groupby(cut).count().SeriousDlqin2yrs-df.groupby(cut).sum().SeriousDlqin2yrs)/(df.SeriousDlqin2yrs.count()-df.SeriousDlqin2yrs.sum())
    iv= (bad_attribute-good_attribute)*cutwoe
    return np.sum(iv)

cut1_iv=get_iv(train_data,cut1,cut1_woe)
cut2_iv=get_iv(train_data,cut2,cut2_woe)
cut3_iv=get_iv(train_data,cut3,cut3_woe)
cut4_iv=get_iv(train_data,cut4,cut4_woe)
cut5_iv=get_iv(train_data,cut5,cut5_woe)
cut6_iv=get_iv(train_data,cut6,cut6_woe)
cut7_iv=get_iv(train_data,cut7,cut7_woe)
cut8_iv=get_iv(train_data,cut8,cut8_woe)
cut9_iv=get_iv(train_data,cut9,cut9_woe)
cut10_iv=get_iv(train_data,cut10,cut10_woe)
cut_iv=pd.DataFrame([cut1_iv,cut2_iv,cut3_iv,cut4_iv,cut5_iv,cut6_iv,cut7_iv,cut8_iv,cut9_iv,cut10_iv],columns=['IV'])

cut_iv.plot.bar(figsize=(8,6))
iv_list=[cut1_iv,cut2_iv,cut3_iv,cut4_iv,cut5_iv,cut6_iv,cut7_iv,cut8_iv,cut9_iv,cut10_iv]
for a,b in zip(np.arange(10),iv_list):
    plt.text(a,b+0.01,'%.3f'%b,ha='center',va='bottom',fontsize=10)
plt.show()

IV值大小的判断，过高的IV可能会有潜在的风险，只能用于二分类的目标变量，IV小于0.02 没有预测性 不用；IV 0.02 to 0.1 弱预测；IV 0.1 to 0.2 一定预测；IV 0.2+ 高预测；将IV值小于0.2的变量删除。

train_data.drop(['SeriousDlqin2yrs','DebtRatio','MonthlyIncome', 'NumberOfOpenCreditLinesAndLoans','NumberRealEstateLoansOrLines','NumberOfDependents'],axis=1,inplace=True)
test= pd.read_csv('TestData.csv')
test.drop(['SeriousDlqin2yrs','DebtRatio','MonthlyIncome', 'NumberOfOpenCreditLinesAndLoans','NumberRealEstateLoansOrLines','NumberOfDependents'],axis=1,inplace=True)

3.5 WOE值替换
目的是减少逻辑回归的自变量处理量

def repalce_woe(series,cutx,woe):
    i=0
    list1=[]
    while i<len(series):
        value = series[i]
        m = len(cutx)-2
        n = len(cutx)-2
        while m>=0:
            if cutx[m]<=value:
                m=m-1
            else:
                m=m-1
                n=n-1
        list1.append(woe[n])
        i=i+1
    return list1   

train_data['RevolvingUtilizationOfUnsecuredLines'] = replace_woe(train_data['RevolvingUtilizationOfUnsecuredLines'].values, cutx1, woex1)
train_data['age'] = Series(replace_woe(train_data['age'].values, cutx2, woex2))
train_data['NumberOfTime30-59DaysPastDueNotWorse'] = replace_woe(train_data['NumberOfTime30-59DaysPastDueNotWorse'].values, cutx3, woex3)
train_data['NumberOfTimes90DaysLate'] = replace_woe(train_data['NumberOfTimes90DaysLate'].values, cutx7, woex7)
train_data['NumberOfTime60-89DaysPastDueNotWorse'] = replace_woe(train_data['NumberOfTime60-89DaysPastDueNotWorse'].values, cutx9, woex9)
train_data.to_csv('WoeData.csv', index=False)

test_data= pd.read_csv('TestData.csv')
test_data['RevolvingUtilizationOfUnsecuredLines'] = replace_woe(test_data['RevolvingUtilizationOfUnsecuredLines'].values, cutx1, woex1)
test_data['age'] = replace_woe(test_data['age'].values, cutx2, woex2)
test_data['NumberOfTime30-59DaysPastDueNotWorse'] = replace_woe(test_data['NumberOfTime30-59DaysPastDueNotWorse'].values, cutx3, woex3)
test_data['NumberOfTimes90DaysLate'] = replace_woe(test_data['NumberOfTimes90DaysLate'].values, cutx7, woex7)
test_data['NumberOfTime60-89DaysPastDueNotWorse'] =replace_woe(test_data['NumberOfTime60-89DaysPastDueNotWorse'].values, cutx9, woex9)
test_data.to_csv('TestWoeData.csv',index=False)

四.逻辑回归模型

from sklearn.linear_model import LogisticRegression
data = pd.read_csv('WoeData.csv')
Y=data['SeriousDlqin2yrs']
X=data.iloc[:,1:]
log= LogisticRegression(solver='lbfgs')
log.fit(X,Y)
test = pd.read_csv('TestWoeData.csv')
Y_test = test['SeriousDlqin2yrs']
X_test = test.iloc[:,1:]
Y_pre1=log.decision_function(X_test)

通过auc和roc_curve 来判断模型好坏

from sklearn.metrics import roc_curve, auc
plt.rcParams['font.sans-serif']=['SimHei']
fpr,tpr,threshold = roc_curve(Y_test,Y_pre1)
rocauc = auc(fpr, tpr)
plt.plot(fpr, tpr, label='AUC = %0.2f' % rocauc)
plt.legend(loc='best')
plt.plot([0, 1], [0, 1], 'r--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel(r'真正率')
plt.xlabel(r'假正率')
plt.show()

AUC值为0.85，说明该模型的预测效果还是不错的，正确率较高。

五. 信用评分模型

说明
由逻辑回归的基本原理，我们将客户违约的概率表示为p，则正常的概率为1-p。因此，可以得到：

此时，客户违约的概率p可表示为：

评分卡设定的分值刻度可以通过将分值表示为比率对数的线性表达式来定义，即可表示为下式：

其中，A和B是常数，式中的负号可以使得违约概率越低，得分越高。通常情况下，这是分值的理想变动方向，即高分值代表低风险，低分值代表高风险。逻辑回归模型计算比率如下所示：

其中，用建模参数拟合模型可以得到模型参数β0 ,β1…βn.
常数A,B值可以通过将两个已知或假设的分值带入到计算得到，通常情况下，需要设定两个假设：
1.给某个特定的比率设定特定的预期分值
2.确定比率翻番的分数(PDO）
假设设定评分卡刻度使得比率为{1:20}（违约正常比）时的分值为600分，PDO为20分，代入式中求得：

评分卡刻度参数A和B确定以后，就可以计算比率和违约概率，以及对应的分值，通常将A称为补偿，常数B称为刻度。
则评分卡的分值可表达为：

式中：变量x 1 …x n是出现在最终模型中的自变量，即为入模指标。由于此时所有变量都用WOE转换进行了转换，可以将这些自变量中的每一个都写(β i ω ij )δ ij 的形式

式中ω ij 为第i行第j个变量的WOE，为已知变量；β i 为逻辑回归方程中的系数，为已知变量；δ ij 为二元变量，表示变量i是否取第j个值。上式可重新表示为：

此式即为最终评分卡公式。如果x 1 …x n变量取不同行并计算其WOE值，式中表示的标准评分卡格式，如表3.20所示

5.1 建立信用评分
这里初始评分设定为600，违约率初始为1:20，当违约上升1倍时，PDO变化为20，则可以进行计算

def get_score(woe,coe,factor):
    scores=[]
    for w in woe:
        score=-round(w*coe*factor,0)
        scores.append(score)
    return scores

def compute_score(series,cut,score):
    list1 = []
    i = 0
    while i < len(series):
        value = series[i]
        j = len(cut) - 2
        m = len(cut) - 2
        while j >= 0:
            if value >= cut[j]:
                j = -1
            else:
                j -= 1
                m -= 1
        list1.append(score[m])
        i += 1
    return list1

log.coef_
log.intercept_

coe=[0.65211317,0.49102952,1.0209049,1.7698133,1.04021588]
B= 20/np.log(2)
A= 600 + B*np.log(120)
base_point = round(A- B*(log.intercept_),0)
# 各项部分分数
x1 = get_score(cut1_woe,coe[0],B) #'RevolvingUtilizationOfUnsecuredLines'
x2 = get_score(cut2_woe,coe[1],B)
#'age'
x3 = get_score(cut3_woe,coe[2],B)
#'NumberOfTime30-59DaysPastDueNotWorse'
x7 = get_score(cut7_woe,coe[3],B)
#'NumberOfTimes90DaysLate'
x9 = get_score(cut9_woe,coe[4],B)
#'NumberOfTime60-89DaysPastDueNotWorse'
print(x1,x2, x3, x7, x9)

不同变量，不同区间的得分如上图所示，在测试集进行总得分的计算，如下

test1 = pd.read_csv('TestData.csv')
test1['BaseScore']=np.zeros(len(test1))+baseScore
test1['x1'] = compute_score(test1['RevolvingUtilizationOfUnsecuredLines'], cutx1, x1)
test1['x2'] = compute_score(test1['age'], cutx2, x2)
test1['x3'] = compute_score(test1['NumberOfTime30-59DaysPastDueNotWorse'], cutx3, x3)
test1['x7'] = compute_score(test1['NumberOfTimes90DaysLate'], cutx7, x7)
test1['x9'] = compute_score(test1['NumberOfTime60-89DaysPastDueNotWorse'], cutx9, x9)
test1['Score'] = test1['x1'] + test1['x2'] + test1['x3'] + test1['x7'] +test1['x9']  + baseScore

这样完成了整个评分卡的建立，分数越高，违约的风险越低。反之，风险越高。