README.md
In andrija-djurovic/PDtoolkit: Collection of Tools for PD Rating Model Development and Validation
PDtoolkit 1.2.1
PDtoolkit provides collection of tools for probability of default (PD)
rating model development and validation. Keeping in mind the fact
that model development is highly iterative and repetitive process,
having standardized and automated tools for this purpose is of the
utmost importance for analysts. The main goal of this package is to
cover the most common steps of PD model development. As additional
contribution we attempted to add some functionalities which at the
moment of the package development were not presented in R package
ecosystem in area of credit risk modeling. The package available
functionalities are those that refer to univariate, bivariate,
multivariate analysis, calibration and validation. Along with
accompanied monobin and monobinShiny packages, PDtoolkit provides
functions which are suitable for different data transformation and
modeling tasks such as: imputations, monotonic binning of numeric risk
factors, binning of categorical risk factors, weights of evidence (WoE)
and information value (IV) calculations, WoE coding (replacement of risk
factors’ modalities with WoE values), risk factor clustering, area under
curve (AUC) calculation and others. Beside mentioned features,
package provides set of validation functions for testing homogeneity,
heterogeneity, discriminatory and predictive power of the model.
The following case study shows the usage of PDtoolkit package. The
study is based on publicly available German credit data set which is
available and downloaded from this
link, but also
distributed along with PDtoolkit package under the data frame loans.
Presented examples are simplified, but yet realistic, version of the PD
model development and validation. In the first row, simplification is
assumed in final model selection, master scale creation, calibration and
some validation steps, whereas in practice we usually have highly
iterative process with couple of available options. Anyhow, this does
not diminish the practical value of the presented examples. So, let’s
start with a concrete case study.
First, we will import the libraries needed for the following examples.
library(PDtoolkit)
library(monobin)
library(rpart)
If the packages are not already installed, run the following code before
the library import
install.packages(c("PDtoolkit", "monobin", "rpart")) in order to
install it from the CRAN, while
devtools::install_github("andrija-djurovic/PDtoolkit") to install
development version.
Then, let’s import and inspect the structure of the modeling data set -
loans.
data(loans)
str(loans)
## 'data.frame': 1000 obs. of 21 variables:
## $ Creditability : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Account Balance : chr "1" "1" "2" "1" ...
## $ Duration of Credit (month) : num 18 9 12 12 12 10 8 6 18 24 ...
## $ Payment Status of Previous Credit: chr "4" "4" "2" "4" ...
## $ Purpose : chr "2" "0" "9" "0" ...
## $ Credit Amount : num 1049 2799 841 2122 2171 ...
## $ Value Savings/Stocks : chr "1" "1" "2" "1" ...
## $ Length of current employment : chr "2" "3" "4" "3" ...
## $ Instalment per cent : chr "4" "2" "2" "3" ...
## $ Sex & Marital Status : chr "2" "3" "2" "3" ...
## $ Guarantors : chr "1" "1" "1" "1" ...
## $ Duration in Current address : chr "4" "2" "4" "2" ...
## $ Most valuable available asset : chr "2" "1" "1" "1" ...
## $ Age (years) : num 21 36 23 39 38 48 39 40 65 23 ...
## $ Concurrent Credits : chr "3" "3" "3" "3" ...
## $ Type of apartment : chr "1" "1" "1" "1" ...
## $ No of Credits at this Bank : chr "1" "2" "1" "2" ...
## $ Occupation : chr "3" "3" "2" "2" ...
## $ No of dependents : chr "1" "2" "1" "2" ...
## $ Telephone : chr "1" "1" "1" "1" ...
## $ Foreign Worker : chr "1" "1" "1" "2" ...
Note that variable (column) names have some special characters which is
not usually considered as good practice in real studies, but for this
example we will use it as provided originally in data source.
Additionally, in order to present all functionalities of the PDtoolkit
package, we will slightly modify modeling data set inserting some
missing values for certain variables as given below:
loans$"Credit Amount"[1:10] <- NA
loans$"Age (years)"[30:50] <- NA
loans$"Purpose"[800:820] <- NA
Usually the first step in model development is the univariate analysis.
For this purpose we will use univariate function from the PDtoolkit
package:
univariate(db = loans)
## rf rf.type bin.type bin cnt pct
## 1 Creditability numeric complete cases complete cases 1000 1.000
## 2 Account Balance character complete cases complete cases 1000 1.000
## 3 Duration of Credit (month) numeric complete cases complete cases 1000 1.000
## 4 Payment Status of Previous Credit character complete cases complete cases 1000 1.000
## 5 Purpose character special cases special cases 21 0.021
## 6 Purpose character complete cases complete cases 979 0.979
## 7 Credit Amount numeric special cases special cases 10 0.010
## 8 Credit Amount numeric complete cases complete cases 990 0.990
## 9 Value Savings/Stocks character complete cases complete cases 1000 1.000
## 10 Length of current employment character complete cases complete cases 1000 1.000
## 11 Instalment per cent character complete cases complete cases 1000 1.000
## 12 Sex & Marital Status character complete cases complete cases 1000 1.000
## 13 Guarantors character complete cases complete cases 1000 1.000
## 14 Duration in Current address character complete cases complete cases 1000 1.000
## 15 Most valuable available asset character complete cases complete cases 1000 1.000
## 16 Age (years) numeric special cases special cases 21 0.021
## 17 Age (years) numeric complete cases complete cases 979 0.979
## 18 Concurrent Credits character complete cases complete cases 1000 1.000
## 19 Type of apartment character complete cases complete cases 1000 1.000
## 20 No of Credits at this Bank character complete cases complete cases 1000 1.000
## 21 Occupation character complete cases complete cases 1000 1.000
## 22 No of dependents character complete cases complete cases 1000 1.000
## 23 Telephone character complete cases complete cases 1000 1.000
## 24 Foreign Worker character complete cases complete cases 1000 1.000
## cnt.unique min p1 p5 p25 p50 avg avg.se p75 p95 p99
## 1 2 0 0.00 0.00 0.00 0 0.30000 0.01449863 1.00 1.0 1.00
## 2 4 NA NA NA NA NA NA NA NA NA NA
## 3 33 4 6.00 6.00 12.00 18 20.90300 0.38133320 24.00 48.0 60.00
## 4 5 NA NA NA NA NA NA NA NA NA NA
## 5 1 NA NA NA NA NA NA NA NA NA NA
## 6 10 NA NA NA NA NA NA NA NA NA NA
## 7 1 Inf NA NA NA NA NaN NA NA NA NA
## 8 915 250 424.13 708.45 1371.25 2324 3283.24242 90.03256824 3978.25 9219.7 14194.29
## 9 5 NA NA NA NA NA NA NA NA NA NA
## 10 5 NA NA NA NA NA NA NA NA NA NA
## 11 4 NA NA NA NA NA NA NA NA NA NA
## 12 4 NA NA NA NA NA NA NA NA NA NA
## 13 3 NA NA NA NA NA NA NA NA NA NA
## 14 4 NA NA NA NA NA NA NA NA NA NA
## 15 4 NA NA NA NA NA NA NA NA NA NA
## 16 1 Inf NA NA NA NA NaN NA NA NA NA
## 17 53 19 20.00 22.00 27.00 33 35.52605 0.36315861 42.00 60.0 67.22
## 18 3 NA NA NA NA NA NA NA NA NA NA
## 19 3 NA NA NA NA NA NA NA NA NA NA
## 20 4 NA NA NA NA NA NA NA NA NA NA
## 21 4 NA NA NA NA NA NA NA NA NA NA
## 22 2 NA NA NA NA NA NA NA NA NA NA
## 23 2 NA NA NA NA NA NA NA NA NA NA
## 24 2 NA NA NA NA NA NA NA NA NA NA
## max neg pos cnt.outliers sc.ind
## 1 1 0 300 0 0
## 2 NA NA NA NA 0
## 3 72 0 1000 70 0
## 4 NA NA NA NA 0
## 5 NA NA NA NA 0
## 6 NA NA NA NA 0
## 7 -Inf NA NA 0 0
## 8 18424 0 990 72 0
## 9 NA NA NA NA 0
## 10 NA NA NA NA 0
## 11 NA NA NA NA 0
## 12 NA NA NA NA 0
## 13 NA NA NA NA 0
## 14 NA NA NA NA 0
## 15 NA NA NA NA 0
## 16 -Inf NA NA 0 0
## 17 75 0 979 23 0
## 18 NA NA NA NA 0
## 19 NA NA NA NA 0
## 20 NA NA NA NA 0
## 21 NA NA NA NA 0
## 22 NA NA NA NA 0
## 23 NA NA NA NA 0
## 24 NA NA NA NA 0
Based on the results we can see that univariate treats differently
so-called special and complete cases. For result details and additional
arguments of univariate function check help page ?univariate.
From the structure and univariate results, we conclude that there are
four numeric variables, while the rest are categorical. One of the
numeric variables Creditability presents binary indicator (0/1) of
default status which will serve as our dependent variable for PD model
development. The rest of the variables present potential risk factors
for the PD model. Additionally, we see that some risk factors
(Purpose, Credit Amount, Age (years)) have certain share of
special cases and that for numeric risk factors potential outliers are
identified. Sometimes, analysts want to impute certain values instead of
special cases (usually missing values) as well as for potential
outliers. Without going too much into deeper analysis for this case
study we can perform imputations using two functions from
PDtoolkitpackage: imp.sc for special cases and imp.outliers for
outliers.
#imputation for special cases for risk factors "Credit Amount", "Purpose"
imp.sc.res <- imp.sc(db = loans[, c("Credit Amount", "Purpose")],
method.num = "automatic",
p.val = 0.05)
names(imp.sc.res)
## [1] "db" "report"
#new risk factors with imputed values
head(imp.sc.res[[1]])
## Credit Amount Purpose
## 1 2324 2
## 2 2324 0
## 3 2324 9
## 4 2324 0
## 5 2324 0
## 6 2324 0
#imputation report
imp.sc.res[[2]]
## rf info imputation.method imputed.value imputation.num
## 1 Credit Amount Imputation completed. automatic - median 2324 10
## 2 Purpose Imputation completed. mode NA 21
## imputed.mode
## 1 <NA>
## 2 3
#replace Credit Amount and Purpose with new values
loans[, c("Credit Amount", "Purpose")] <- imp.sc.res[[1]]
colSums(is.na(loans))
## Creditability Account Balance
## 0 0
## Duration of Credit (month) Payment Status of Previous Credit
## 0 0
## Purpose Credit Amount
## 0 0
## Value Savings/Stocks Length of current employment
## 0 0
## Instalment per cent Sex & Marital Status
## 0 0
## Guarantors Duration in Current address
## 0 0
## Most valuable available asset Age (years)
## 0 21
## Concurrent Credits Type of apartment
## 0 0
## No of Credits at this Bank Occupation
## 0 0
## No of dependents Telephone
## 0 0
## Foreign Worker
## 0
#imputation for outliers for risk factor Credit Amount
imp.out.res <- imp.outliers(db = loans[, "Credit Amount", drop = FALSE],
method = "iqr",
range = 1.5)
#new risk factors with imputed values
head(imp.out.res[[1]])
## Credit Amount
## 1 2324
## 2 2324
## 3 2324
## 4 2324
## 5 2324
## 6 2324
#imputation report
imp.out.res[[2]]
## rf info imputation.method imputation.val.upper
## 1 Credit Amount Imputation completed. iqr 7865
## imputation.val.lower imputation.num.upper imputation.num.lower
## 1 250 73 0
#replace Credit Amount with new values
loans[, "Credit Amount"] <- imp.out.res[[1]]
Usually, when building PD models numeric risk factors are discretized,
so we will proceed next with that step. For the purpose of binning the
numeric risk factors, we will use one of the functions (ndr.bin) from
the monobin package. Details about this package can be found
here.
#identify numeric risk factors
num.rf <- sapply(loans, is.numeric)
num.rf <- names(num.rf)[!names(num.rf)%in%"Creditability" & num.rf]
#discretized numeric risk factors using ndr.bin from monobin package
loans[, num.rf] <- sapply(num.rf, function(x)
monobin::ndr.bin(x = loans[, x], y = loans[, "Creditability"])[[2]])
#check loans structure again and confirm binning was sucessful
str(loans)
## 'data.frame': 1000 obs. of 21 variables:
## $ Creditability : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Account Balance : chr "1" "1" "2" "1" ...
## $ Duration of Credit (month) : chr "03 [16,45)" "02 [8,16)" "02 [8,16)" "02 [8,16)" ...
## $ Payment Status of Previous Credit: chr "4" "4" "2" "4" ...
## $ Purpose : chr "2" "0" "9" "0" ...
## $ Credit Amount : chr "01 (-Inf,3914)" "01 (-Inf,3914)" "01 (-Inf,3914)" "01 (-Inf,3914)" ...
## $ Value Savings/Stocks : chr "1" "1" "2" "1" ...
## $ Length of current employment : chr "2" "3" "4" "3" ...
## $ Instalment per cent : chr "4" "2" "2" "3" ...
## $ Sex & Marital Status : chr "2" "3" "2" "3" ...
## $ Guarantors : chr "1" "1" "1" "1" ...
## $ Duration in Current address : chr "4" "2" "4" "2" ...
## $ Most valuable available asset : chr "2" "1" "1" "1" ...
## $ Age (years) : chr "01 (-Inf,26)" "03 [35,Inf)" "01 (-Inf,26)" "03 [35,Inf)" ...
## $ Concurrent Credits : chr "3" "3" "3" "3" ...
## $ Type of apartment : chr "1" "1" "1" "1" ...
## $ No of Credits at this Bank : chr "1" "2" "1" "2" ...
## $ Occupation : chr "3" "3" "2" "2" ...
## $ No of dependents : chr "1" "2" "1" "2" ...
## $ Telephone : chr "1" "1" "1" "1" ...
## $ Foreign Worker : chr "1" "1" "1" "2" ...
As we can see after the binning, all risk factors are categorical. Now
the modeling data set is ready for bivariate analysis. For that purpose
we will use the function bivariate from the PDtoolkitpackage.
bivariate(db = loans, target = "Creditability")
## $results
## rf bin no ng nb pct.o pct.g
## 1 Account Balance 1 274 139 135 0.274 0.198571429
## 2 Account Balance 2 269 164 105 0.269 0.234285714
## 3 Account Balance 3 63 49 14 0.063 0.070000000
## 4 Account Balance 4 394 348 46 0.394 0.497142857
## 5 Duration of Credit (month) 01 (-Inf,8) 87 78 9 0.087 0.111428571
## 6 Duration of Credit (month) 02 [8,16) 344 264 80 0.344 0.377142857
## 7 Duration of Credit (month) 03 [16,45) 499 328 171 0.499 0.468571429
## 8 Duration of Credit (month) 04 [45,Inf) 70 30 40 0.070 0.042857143
## 9 Payment Status of Previous Credit 0 40 15 25 0.040 0.021428571
## 10 Payment Status of Previous Credit 1 49 21 28 0.049 0.030000000
## 11 Payment Status of Previous Credit 2 530 361 169 0.530 0.515714286
## 12 Payment Status of Previous Credit 3 88 60 28 0.088 0.085714286
## 13 Payment Status of Previous Credit 4 293 243 50 0.293 0.347142857
## 14 Purpose 0 224 145 79 0.224 0.207142857
## 15 Purpose 1 102 86 16 0.102 0.122857143
## 16 Purpose 10 12 7 5 0.012 0.010000000
## 17 Purpose 2 177 123 54 0.177 0.175714286
## 18 Purpose 3 300 218 82 0.300 0.311428571
## 19 Purpose 4 12 8 4 0.012 0.011428571
## 20 Purpose 5 21 14 7 0.021 0.020000000
## 21 Purpose 6 49 28 21 0.049 0.040000000
## 22 Purpose 8 9 8 1 0.009 0.011428571
## 23 Purpose 9 94 63 31 0.094 0.090000000
## 24 Credit Amount 01 (-Inf,3914) 739 550 189 0.739 0.785714286
## 25 Credit Amount 02 [3914,Inf) 261 150 111 0.261 0.214285714
## 26 Value Savings/Stocks 1 603 386 217 0.603 0.551428571
## 27 Value Savings/Stocks 2 103 69 34 0.103 0.098571429
## 28 Value Savings/Stocks 3 63 52 11 0.063 0.074285714
## 29 Value Savings/Stocks 4 48 42 6 0.048 0.060000000
## 30 Value Savings/Stocks 5 183 151 32 0.183 0.215714286
## 31 Length of current employment 1 62 39 23 0.062 0.055714286
## 32 Length of current employment 2 172 102 70 0.172 0.145714286
## 33 Length of current employment 3 339 235 104 0.339 0.335714286
## 34 Length of current employment 4 174 135 39 0.174 0.192857143
## 35 Length of current employment 5 253 189 64 0.253 0.270000000
## 36 Instalment per cent 1 136 102 34 0.136 0.145714286
## 37 Instalment per cent 2 231 169 62 0.231 0.241428571
## 38 Instalment per cent 3 157 112 45 0.157 0.160000000
## 39 Instalment per cent 4 476 317 159 0.476 0.452857143
## 40 Sex & Marital Status 1 50 30 20 0.050 0.042857143
## 41 Sex & Marital Status 2 310 201 109 0.310 0.287142857
## 42 Sex & Marital Status 3 548 402 146 0.548 0.574285714
## 43 Sex & Marital Status 4 92 67 25 0.092 0.095714286
## 44 Guarantors 1 907 635 272 0.907 0.907142857
## 45 Guarantors 2 41 23 18 0.041 0.032857143
## 46 Guarantors 3 52 42 10 0.052 0.060000000
## 47 Duration in Current address 1 130 94 36 0.130 0.134285714
## 48 Duration in Current address 2 308 211 97 0.308 0.301428571
## 49 Duration in Current address 3 149 106 43 0.149 0.151428571
## 50 Duration in Current address 4 413 289 124 0.413 0.412857143
## 51 Most valuable available asset 1 282 222 60 0.282 0.317142857
## 52 Most valuable available asset 2 232 161 71 0.232 0.230000000
## 53 Most valuable available asset 3 332 230 102 0.332 0.328571429
## 54 Most valuable available asset 4 154 87 67 0.154 0.124285714
## 55 Age (years) 01 (-Inf,26) 187 108 79 0.187 0.154285714
## 56 Age (years) 02 [26,35) 349 238 111 0.349 0.340000000
## 57 Age (years) 03 [35,Inf) 443 335 108 0.443 0.478571429
## 58 Age (years) SC 21 19 2 0.021 0.027142857
## 59 Concurrent Credits 1 139 82 57 0.139 0.117142857
## 60 Concurrent Credits 2 47 28 19 0.047 0.040000000
## 61 Concurrent Credits 3 814 590 224 0.814 0.842857143
## 62 Type of apartment 1 179 109 70 0.179 0.155714286
## 63 Type of apartment 2 714 528 186 0.714 0.754285714
## 64 Type of apartment 3 107 63 44 0.107 0.090000000
## 65 No of Credits at this Bank 1 633 433 200 0.633 0.618571429
## 66 No of Credits at this Bank 2 333 241 92 0.333 0.344285714
## 67 No of Credits at this Bank 3 28 22 6 0.028 0.031428571
## 68 No of Credits at this Bank 4 6 4 2 0.006 0.005714286
## 69 Occupation 1 22 15 7 0.022 0.021428571
## 70 Occupation 2 200 144 56 0.200 0.205714286
## 71 Occupation 3 630 444 186 0.630 0.634285714
## 72 Occupation 4 148 97 51 0.148 0.138571429
## 73 No of dependents 1 845 591 254 0.845 0.844285714
## 74 No of dependents 2 155 109 46 0.155 0.155714286
## 75 Telephone 1 596 409 187 0.596 0.584285714
## 76 Telephone 2 404 291 113 0.404 0.415714286
## 77 Foreign Worker 1 963 667 296 0.963 0.952857143
## 78 Foreign Worker 2 37 33 4 0.037 0.047142857
## pct.b dr so sg sb dist.g dist.b woe iv.b
## 1 0.450000000 0.4927007 1000 700 300 0.198571429 0.450000000 -0.8180987057 0.2056933888604
## 2 0.350000000 0.3903346 1000 700 300 0.234285714 0.350000000 -0.4013917827 0.0464467634291
## 3 0.046666667 0.2222222 1000 700 300 0.070000000 0.046666667 0.4054651081 0.0094608525225
## 4 0.153333333 0.1167513 1000 700 300 0.497142857 0.153333333 1.1762632229 0.4044104985393
## 5 0.030000000 0.1034483 1000 700 300 0.111428571 0.030000000 1.3121863890 0.1068494631015
## 6 0.266666667 0.2325581 1000 700 300 0.377142857 0.266666667 0.3466246081 0.0382937662266
## 7 0.570000000 0.3426854 1000 700 300 0.468571429 0.570000000 -0.1959478085 0.0198747062913
## 8 0.133333333 0.5714286 1000 700 300 0.042857143 0.133333333 -1.1349799328 0.1026886605902
## 9 0.083333333 0.6250000 1000 700 300 0.021428571 0.083333333 -1.3581234842 0.0840743109238
## 10 0.093333333 0.5714286 1000 700 300 0.030000000 0.093333333 -1.1349799328 0.0718820624131
## 11 0.563333333 0.3188679 1000 700 300 0.515714286 0.563333333 -0.0883186170 0.0042056484275
## 12 0.093333333 0.3181818 1000 700 300 0.085714286 0.093333333 -0.0851578083 0.0006488213969
## 13 0.166666667 0.1706485 1000 700 300 0.347142857 0.166666667 0.7337405775 0.1324227042295
## 14 0.263333333 0.3526786 1000 700 300 0.207142857 0.263333333 -0.2400119704 0.0134863869101
## 15 0.053333333 0.1568627 1000 700 300 0.122857143 0.053333333 0.8344607136 0.0580148877093
## 16 0.016666667 0.4166667 1000 700 300 0.010000000 0.016666667 -0.5108256238 0.0034055041584
## 17 0.180000000 0.3050847 1000 700 300 0.175714286 0.180000000 -0.0240975516 0.0001032752211
## 18 0.273333333 0.2733333 1000 700 300 0.311428571 0.273333333 0.1304779551 0.0049705887671
## 19 0.013333333 0.3333333 1000 700 300 0.011428571 0.013333333 -0.1541506798 0.0002936203425
## 20 0.023333333 0.3333333 1000 700 300 0.020000000 0.023333333 -0.1541506798 0.0005138355994
## 21 0.070000000 0.4285714 1000 700 300 0.040000000 0.070000000 -0.5596157879 0.0167884736381
## 22 0.003333333 0.1111111 1000 700 300 0.011428571 0.003333333 1.2321436813 0.0099744964676
## 23 0.103333333 0.3297872 1000 700 300 0.090000000 0.103333333 -0.1381503385 0.0018420045131
## 24 0.630000000 0.2557510 1000 700 300 0.785714286 0.630000000 0.2208734028 0.0343931441471
## 25 0.370000000 0.4252874 1000 700 300 0.214285714 0.370000000 -0.5461927676 0.0850500166697
## 26 0.723333333 0.3598673 1000 700 300 0.551428571 0.723333333 -0.2713578445 0.0466477056434
## 27 0.113333333 0.3300971 1000 700 300 0.098571429 0.113333333 -0.1395518804 0.0020600515679
## 28 0.036666667 0.1746032 1000 700 300 0.074285714 0.036666667 0.7060505854 0.0265609505935
## 29 0.020000000 0.1250000 1000 700 300 0.060000000 0.020000000 1.0986122887 0.0439444915467
## 30 0.106666667 0.1748634 1000 700 300 0.215714286 0.106666667 0.7042460736 0.0767963575528
## 31 0.076666667 0.3709677 1000 700 300 0.055714286 0.076666667 -0.3192304302 0.0066886375849
## 32 0.233333333 0.4069767 1000 700 300 0.145714286 0.233333333 -0.4708202892 0.0412528253352
## 33 0.346666667 0.3067847 1000 700 300 0.335714286 0.346666667 -0.0321032454 0.0003516069733
## 34 0.130000000 0.2241379 1000 700 300 0.192857143 0.130000000 0.3944152719 0.0247918170922
## 35 0.213333333 0.2529644 1000 700 300 0.270000000 0.213333333 0.2355660713 0.0133487440411
## 36 0.113333333 0.2500000 1000 700 300 0.145714286 0.113333333 0.2513144283 0.0081378005348
## 37 0.206666667 0.2683983 1000 700 300 0.241428571 0.206666667 0.1554664695 0.0054043106061
## 38 0.150000000 0.2866242 1000 700 300 0.160000000 0.150000000 0.0645385211 0.0006453852114
## 39 0.530000000 0.3340336 1000 700 300 0.452857143 0.530000000 -0.1573002887 0.0121345937020
## 40 0.066666667 0.4000000 1000 700 300 0.042857143 0.066666667 -0.4418327523 0.0105198274352
## 41 0.363333333 0.3516129 1000 700 300 0.287142857 0.363333333 -0.2353408346 0.0179307302520
## 42 0.486666667 0.2664234 1000 700 300 0.574285714 0.486666667 0.1655476065 0.0145051236192
## 43 0.083333333 0.2717391 1000 700 300 0.095714286 0.083333333 0.1385189341 0.0017149963274
## 44 0.906666667 0.2998897 1000 700 300 0.907142857 0.906666667 0.0005250722 0.0000002500344
## 45 0.060000000 0.4390244 1000 700 300 0.032857143 0.060000000 -0.6021754024 0.0163447609210
## 46 0.033333333 0.1923077 1000 700 300 0.060000000 0.033333333 0.5877866649 0.0156743110641
## 47 0.120000000 0.2769231 1000 700 300 0.134285714 0.120000000 0.1124779834 0.0016068283347
## 48 0.323333333 0.3149351 1000 700 300 0.301428571 0.323333333 -0.0701507054 0.0015366344996
## 49 0.143333333 0.2885906 1000 700 300 0.151428571 0.143333333 0.0549411180 0.0004447614317
## 50 0.413333333 0.3002421 1000 700 300 0.412857143 0.413333333 -0.0011527379 0.0000005489228
## 51 0.200000000 0.2127660 1000 700 300 0.317142857 0.200000000 0.4610349593 0.0540069523708
## 52 0.236666667 0.3060345 1000 700 300 0.230000000 0.236666667 -0.0285733724 0.0001904891496
## 53 0.340000000 0.3072289 1000 700 300 0.328571429 0.340000000 -0.0341913647 0.0003907584543
## 54 0.223333333 0.4350649 1000 700 300 0.124285714 0.223333333 -0.5860823611 0.0580500624351
## 55 0.263333333 0.4224599 1000 700 300 0.154285714 0.263333333 -0.5346144857 0.0582984367772
## 56 0.370000000 0.3180516 1000 700 300 0.340000000 0.370000000 -0.0845573880 0.0025367216408
## 57 0.360000000 0.2437923 1000 700 300 0.478571429 0.360000000 0.2847014443 0.0337574569686
## 58 0.006666667 0.0952381 1000 700 300 0.027142857 0.006666667 1.4039939382 0.0287484473064
## 59 0.190000000 0.4100719 1000 700 300 0.117142857 0.190000000 -0.4836298810 0.0352358913269
## 60 0.063333333 0.4042553 1000 700 300 0.040000000 0.063333333 -0.4595323294 0.0107224210188
## 61 0.746666667 0.2751843 1000 700 300 0.842857143 0.746666667 0.1211786247 0.0116562296099
## 62 0.233333333 0.3910615 1000 700 300 0.155714286 0.233333333 -0.4044452202 0.0313926528066
## 63 0.620000000 0.2605042 1000 700 300 0.754285714 0.620000000 0.1960517496 0.0263269492328
## 64 0.146666667 0.4112150 1000 700 300 0.090000000 0.146666667 -0.4883527679 0.0276733235151
## 65 0.666666667 0.3159558 1000 700 300 0.618571429 0.666666667 -0.0748774989 0.0036012511391
## 66 0.306666667 0.2762763 1000 700 300 0.344285714 0.306666667 0.1157104961 0.0043529186611
## 67 0.020000000 0.2142857 1000 700 300 0.031428571 0.020000000 0.4519851237 0.0051655442713
## 68 0.006666667 0.3333333 1000 700 300 0.005714286 0.006666667 -0.1541506798 0.0001468101713
## 69 0.023333333 0.3181818 1000 700 300 0.021428571 0.023333333 -0.0851578083 0.0001622053492
## 70 0.186666667 0.2800000 1000 700 300 0.205714286 0.186666667 0.0971637485 0.0018507380658
## 71 0.620000000 0.2952381 1000 700 300 0.634285714 0.620000000 0.0227800283 0.0003254289762
## 72 0.170000000 0.3445946 1000 700 300 0.138571429 0.170000000 -0.2044125146 0.0064243933163
## 73 0.846666667 0.3005917 1000 700 300 0.844285714 0.846666667 -0.0028161100 0.0000067050238
## 74 0.153333333 0.2967742 1000 700 300 0.155714286 0.153333333 0.0154086254 0.0000366872032
## 75 0.623333333 0.3137584 1000 700 300 0.584285714 0.623333333 -0.0646913212 0.0025260420659
## 76 0.376666667 0.2797030 1000 700 300 0.415714286 0.376666667 0.0986375881 0.0038515629628
## 77 0.986666667 0.3073728 1000 700 300 0.952857143 0.986666667 -0.0348672688 0.0011788457545
## 78 0.013333333 0.1081081 1000 700 300 0.047142857 0.013333333 1.2629153400 0.0426985662558
## iv.s auc
## 1 0.66601150335 0.7077690
## 2 0.66601150335 0.7077690
## 3 0.66601150335 0.7077690
## 4 0.66601150335 0.7077690
## 5 0.26770659621 0.6241762
## 6 0.26770659621 0.6241762
## 7 0.26770659621 0.6241762
## 8 0.26770659621 0.6241762
## 9 0.29323354739 0.6268048
## 10 0.29323354739 0.6268048
## 11 0.29323354739 0.6268048
## 12 0.29323354739 0.6268048
## 13 0.29323354739 0.6268048
## 14 0.10939307333 0.5821500
## 15 0.10939307333 0.5821500
## 16 0.10939307333 0.5821500
## 17 0.10939307333 0.5821500
## 18 0.10939307333 0.5821500
## 19 0.10939307333 0.5821500
## 20 0.10939307333 0.5821500
## 21 0.10939307333 0.5821500
## 22 0.10939307333 0.5821500
## 23 0.10939307333 0.5821500
## 24 0.11944316082 0.5778571
## 25 0.11944316082 0.5778571
## 26 0.19600955690 0.5991429
## 27 0.19600955690 0.5991429
## 28 0.19600955690 0.5991429
## 29 0.19600955690 0.5991429
## 30 0.19600955690 0.5991429
## 31 0.08643363103 0.5808190
## 32 0.08643363103 0.5808190
## 33 0.08643363103 0.5808190
## 34 0.08643363103 0.5808190
## 35 0.08643363103 0.5808190
## 36 0.02632209005 0.5433833
## 37 0.02632209005 0.5433833
## 38 0.02632209005 0.5433833
## 39 0.02632209005 0.5433833
## 40 0.04467067763 0.5524238
## 41 0.04467067763 0.5524238
## 42 0.04467067763 0.5524238
## 43 0.04467067763 0.5524238
## 44 0.03201932202 0.5256524
## 45 0.03201932202 0.5256524
## 46 0.03201932202 0.5256524
## 47 0.00358877319 0.5161786
## 48 0.00358877319 0.5161786
## 49 0.00358877319 0.5161786
## 50 0.00358877319 0.5161786
## 51 0.11263826241 0.5853286
## 52 0.11263826241 0.5853286
## 53 0.11263826241 0.5853286
## 54 0.11263826241 0.5853286
## 55 0.12334106269 0.5890381
## 56 0.12334106269 0.5890381
## 57 0.12334106269 0.5890381
## 58 0.12334106269 0.5890381
## 59 0.05761454196 0.5481857
## 60 0.05761454196 0.5481857
## 61 0.05761454196 0.5481857
## 62 0.08539292555 0.5680619
## 63 0.08539292555 0.5680619
## 64 0.08539292555 0.5680619
## 65 0.01326652424 0.5260571
## 66 0.01326652424 0.5260571
## 67 0.01326652424 0.5260571
## 68 0.01326652424 0.5260571
## 69 0.00876276571 0.5214429
## 70 0.00876276571 0.5214429
## 71 0.00876276571 0.5214429
## 72 0.00876276571 0.5214429
## 73 0.00004339223 0.5011905
## 74 0.00004339223 0.5011905
## 75 0.00637760503 0.5195238
## 76 0.00637760503 0.5195238
## 77 0.04387741201 0.5169048
## 78 0.04387741201 0.5169048
##
## $info
## data frame with 0 columns and 0 rows
Based on results of the bivariate analysis we can see that risk factor
Purpose has 10 modalities, that Age (years) has share of special
cases 2.1% and that there are quite some other risk factors with
modality share less than 5%. In order to correct above potential issues,
we can further categorize risk factors Purpose, Foreign Worker
and Age (years) given the following criteria:
- each modality has at least 5% of the observations or minimum 30
observations,
- each modality has at least 1% of defaults or minimum 1 default,
- the maximum number of modalities cannot be greater than 5,
- if special cases exist to be merged with a closest modality (based
on default rate).
For this purpose we will use function cat.bin. Additional examples of
cat.bin function usage can be found in its help page ?cat.bin.
Note: Usually in the practice greater attention is paid on this step
and all risk factors are examined closer.
rf <- c("Purpose", "Age (years)", "Foreign Worker")
loans[, rf] <- sapply(rf, function(x)
cat.bin(x = loans[, x],
y = loans[, "Creditability"],
sc = c("SC", NA),
sc.merge = "closest",
min.pct.obs = 0.05,
min.avg.rate = 0.01,
max.groups = 5,
force.trend = "dr")[[2]])
#run bivariate analysis again
bivariate(db = loans[, c(rf, "Creditability")], target = "Creditability")
## $results
## rf bin no ng nb pct.o pct.g pct.b dr so
## 1 Purpose 1 [1,8] 111 94 17 0.111 0.1342857 0.05666667 0.1531532 1000
## 2 Purpose 2 [3] 300 218 82 0.300 0.3114286 0.27333333 0.2733333 1000
## 3 Purpose 3 [2] 177 123 54 0.177 0.1757143 0.18000000 0.3050847 1000
## 4 Purpose 4 [0,4,5,9] 351 230 121 0.351 0.3285714 0.40333333 0.3447293 1000
## 5 Purpose 5 [10,6] 61 35 26 0.061 0.0500000 0.08666667 0.4262295 1000
## 6 Age (years) 1 [03 [35,Inf)] 464 354 110 0.464 0.5057143 0.36666667 0.2370690 1000
## 7 Age (years) 2 [02 [26,35)] 349 238 111 0.349 0.3400000 0.37000000 0.3180516 1000
## 8 Age (years) 3 [01 (-Inf,26)] 187 108 79 0.187 0.1542857 0.26333333 0.4224599 1000
## 9 Foreign Worker 1 [2,1] 1000 700 300 1.000 1.0000000 1.00000000 0.3000000 1000
## sg sb dist.g dist.b woe iv.b iv.s auc
## 1 700 300 0.1342857 0.05666667 0.86278358 0.0669684396 0.1075375 0.5805190
## 2 700 300 0.3114286 0.27333333 0.13047796 0.0049705888 0.1075375 0.5805190
## 3 700 300 0.1757143 0.18000000 -0.02409755 0.0001032752 0.1075375 0.5805190
## 4 700 300 0.3285714 0.40333333 -0.20500910 0.0153268706 0.1075375 0.5805190
## 5 700 300 0.0500000 0.08666667 -0.55004634 0.0201683657 0.1075375 0.5805190
## 6 700 300 0.5057143 0.36666667 0.32151869 0.0447064079 0.1055416 0.5857476
## 7 700 300 0.3400000 0.37000000 -0.08455739 0.0025367216 0.1055416 0.5857476
## 8 700 300 0.1542857 0.26333333 -0.53461449 0.0582984368 0.1055416 0.5857476
## 9 700 300 1.0000000 1.00000000 0.00000000 0.0000000000 0.0000000 NA
##
## $info
## data frame with 0 columns and 0 rows
After additional binning based on above criteria and running again the
bivariate analysis, we see that risk factor Foreign Worker has only
one value and it can be discarded from the further analysis.
loans <- loans[, !names(loans)%in%"Foreign Worker"]
In case when too many risk factors are available, analysts have an
option to perform cluster analysis and shorten the list of potential
candidates for final PD model. For this purpose, PDtoolkit package
contains function rf.clustering. Since the most metrics needed for
distance matrix calculation require numeric values, we will first
replace risk factors’ modalities with Weights of Evidence (WoE) and then
perform clustering.
#replace modalities with WoE
woe.res <- replace.woe(db = loans, target = "Creditability")
head(woe.res[[1]])
## Creditability Account Balance Duration of Credit (month)
## 1 0 -0.8180987 -0.1959478
## 2 0 -0.8180987 0.3466246
## 3 0 -0.4013918 0.3466246
## 4 0 -0.8180987 0.3466246
## 5 0 -0.8180987 0.3466246
## 6 0 -0.8180987 0.3466246
## Payment Status of Previous Credit Purpose Credit Amount Value Savings/Stocks
## 1 0.73374058 -0.02409755 0.2208734 -0.2713578
## 2 0.73374058 -0.20500910 0.2208734 -0.2713578
## 3 -0.08831862 -0.20500910 0.2208734 -0.1395519
## 4 0.73374058 -0.20500910 0.2208734 -0.2713578
## 5 0.73374058 -0.20500910 0.2208734 -0.2713578
## 6 0.73374058 -0.20500910 0.2208734 -0.2713578
## Length of current employment Instalment per cent Sex & Marital Status Guarantors
## 1 -0.47082029 -0.15730029 -0.2353408 0.0005250722
## 2 -0.03210325 0.15546647 0.1655476 0.0005250722
## 3 0.39441527 0.15546647 -0.2353408 0.0005250722
## 4 -0.03210325 0.06453852 0.1655476 0.0005250722
## 5 -0.03210325 -0.15730029 0.1655476 0.0005250722
## 6 -0.47082029 0.25131443 0.1655476 0.0005250722
## Duration in Current address Most valuable available asset Age (years) Concurrent Credits
## 1 -0.001152738 -0.02857337 -0.5346145 0.1211786
## 2 -0.070150705 0.46103496 0.3215187 0.1211786
## 3 -0.001152738 0.46103496 -0.5346145 0.1211786
## 4 -0.070150705 0.46103496 0.3215187 0.1211786
## 5 -0.001152738 -0.02857337 0.3215187 -0.4836299
## 6 0.054941118 0.46103496 0.3215187 0.1211786
## Type of apartment No of Credits at this Bank Occupation No of dependents Telephone
## 1 -0.4044452 -0.0748775 0.02278003 -0.00281611 -0.06469132
## 2 -0.4044452 0.1157105 0.02278003 0.01540863 -0.06469132
## 3 -0.4044452 -0.0748775 0.09716375 -0.00281611 -0.06469132
## 4 -0.4044452 0.1157105 0.09716375 0.01540863 -0.06469132
## 5 0.1960517 0.1157105 0.09716375 -0.00281611 -0.06469132
## 6 -0.4044452 0.1157105 0.09716375 0.01540863 -0.06469132
loans.woe <- woe.res[[1]]
#perform cluster analysis using "common pearson" metric
#number of clusters is set to NA in order to enable automatic determination of optimal number of clusters using elbow method
rf.clustering(db = loans.woe, metric = "common pearson", k = NA)
## rf clusters dist.to.centroid
## 1 Creditability 1 0.0000000
## 2 Account Balance 2 0.0000000
## 4 Duration of Credit (month) 3 0.5788368
## 6 Occupation 3 0.6001961
## 3 Credit Amount 3 0.6892248
## 7 Type of apartment 3 0.7067914
## 5 Most valuable available asset 3 0.8086682
## 8 Payment Status of Previous Credit 4 0.0000000
## 9 Purpose 5 0.0000000
## 10 Value Savings/Stocks 6 0.0000000
## 11 Length of current employment 7 0.0000000
## 12 Instalment per cent 8 0.0000000
## 13 Sex & Marital Status 9 0.0000000
## 14 Guarantors 10 0.0000000
## 15 Duration in Current address 11 0.0000000
## 16 Age (years) 12 0.0000000
## 17 Concurrent Credits 13 0.0000000
## 18 No of Credits at this Bank 14 0.0000000
## 19 No of dependents 15 0.0000000
## 20 Telephone 16 0.0000000
Beside metrics that require numeric inputs, there is x2y metric that
works with both numeric and categorical data. This metric is especially
handy if analyst wants to perform clustering before any binning
procedures and to decrease the number of risk factors. More examples of
the clustering can be found in help page ?rf.clustering, while the
details about x2y metric are presented in this
link.
Due to the fact that for this example we don’t have too many risk
factors, we will not use the results of cluster analysis to shorten the
list of potential risk factors for the final PD model.
The next step is then to perform multivariate analysis and to select the
final PD model. In order to complete this task we will use function
stepMIV. This function uses the concept of Marginal Information Value
(MIV) in order to select risk factors in the model. Details can be found
on the help page of the function ?stepMIV.
multiv.res <- stepMIV(start.model = Creditability ~ 1,
miv.threshold = 0.02,
m.ch.p.val = 0.05,
coding = "WoE",
db = loans)
## Running iteration: 1
## Running iteration: 2
## Running iteration: 3
## Running iteration: 4
## Running iteration: 5
## Running iteration: 6
## Running iteration: 7
#print the output elements
names(multiv.res)
## [1] "model" "steps" "miv.iter" "warnings" "dev.db"
#print the final model and estimated coefficients
multiv.res$model
##
## Call: glm(formula = as.formula(mod.frm), family = "binomial", data = db)
##
## Coefficients:
## (Intercept) `Account Balance`
## -0.8386 -0.8516
## `Duration of Credit (month)` `Payment Status of Previous Credit`
## -0.9586 -0.7857
## `Value Savings/Stocks` Purpose
## -0.7754 -0.9566
## Guarantors
## -1.1917
##
## Degrees of Freedom: 999 Total (i.e. Null); 993 Residual
## Null Deviance: 1222
## Residual Deviance: 973.6 AIC: 987.6
summary(multiv.res$model)$coefficients
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8385682 0.08070531 -10.390496 2.739256e-25
## `Account Balance` -0.8515981 0.10201758 -8.347562 6.968627e-17
## `Duration of Credit (month)` -0.9586008 0.16294350 -5.883026 4.028332e-09
## `Payment Status of Previous Credit` -0.7856716 0.14969919 -5.248336 1.534793e-07
## `Value Savings/Stocks` -0.7754343 0.19240625 -4.030193 5.573104e-05
## Purpose -0.9565752 0.25273860 -3.784840 1.538075e-04
## Guarantors -1.1916767 0.42403712 -2.810312 4.949344e-03
#store the final model
pd.model <- multiv.res$model
#store the data set used for the model estimation
pd.db <- multiv.res$dev.db
#note that risk factors from the model are coded using WoE due to the fact that WoE coding method is used in stepMIV estimation
head(pd.db)
## Creditability Account Balance Duration of Credit (month)
## 1 0 -0.8180987 -0.1959478
## 2 0 -0.8180987 0.3466246
## 3 0 -0.4013918 0.3466246
## 4 0 -0.8180987 0.3466246
## 5 0 -0.8180987 0.3466246
## 6 0 -0.8180987 0.3466246
## Payment Status of Previous Credit Purpose Credit Amount Value Savings/Stocks
## 1 0.73374058 -0.02409755 01 (-Inf,3914) -0.2713578
## 2 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 3 -0.08831862 -0.20500910 01 (-Inf,3914) -0.1395519
## 4 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 5 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 6 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## Length of current employment Instalment per cent Sex & Marital Status Guarantors
## 1 2 4 2 0.0005250722
## 2 3 2 3 0.0005250722
## 3 4 2 2 0.0005250722
## 4 3 3 3 0.0005250722
## 5 3 4 3 0.0005250722
## 6 2 1 3 0.0005250722
## Duration in Current address Most valuable available asset Age (years)
## 1 4 2 3 [01 (-Inf,26)]
## 2 2 1 1 [03 [35,Inf)]
## 3 4 1 3 [01 (-Inf,26)]
## 4 2 1 1 [03 [35,Inf)]
## 5 4 2 1 [03 [35,Inf)]
## 6 3 1 1 [03 [35,Inf)]
## Concurrent Credits Type of apartment No of Credits at this Bank Occupation
## 1 3 1 1 3
## 2 3 1 2 3
## 3 3 1 1 2
## 4 3 1 2 2
## 5 1 2 2 2
## 6 3 1 2 2
## No of dependents Telephone
## 1 1 1
## 2 2 1
## 3 1 1
## 4 2 1
## 5 1 1
## 6 2 1
#calculate area under curve (AUC) for the final model
auc.model(predictions = predict(pd.model, type = "response", newdata = pd.db),
observed = pd.db$Creditability)
## [1] 0.7983738
In practice multivariate analysis is highly iterative process, so
analysts usually end up with a few candidate models. After that based on
different statistical and business validations of these models, one is
selected as the final. For the further process we will assume that the
model from the previous step is the final and proceed with its
validation. For that purpose, we will apply k-fold cross-validation and
validation based on bootstrapping. Both methods are implemented in
PDtoolkit package, concretely in functions kfold.vld and
boots.vld.
kfold.vld.res <- kfold.vld(model = pd.model, k = 10, seed = 1984)
kfold.vld.res
## $iter
## k no amse rmse auc
## 1 1 100 0.1630258 0.4037645 0.7712482
## 2 2 100 0.1374731 0.3707736 0.8635211
## 3 3 100 0.1908883 0.4369077 0.7076403
## 4 4 100 0.1179182 0.3433922 0.8726708
## 5 5 100 0.1880557 0.4336539 0.7256897
## 6 6 100 0.1637761 0.4046927 0.8037087
## 7 7 100 0.1400881 0.3742835 0.8508985
## 8 8 100 0.1827844 0.4275329 0.8596527
## 9 9 100 0.1910433 0.4370850 0.7354138
## 10 10 100 0.1617250 0.4021504 0.7536383
##
## $summary
## amse rmse auc
## 1 0.1636778 0.4034236 0.7944082
boots.vld.res <- boots.vld (model = pd.model, B = 10, seed = 1122)
head(boots.vld.res[[1]])
## B amse rmse auc
## 1 1 0.1555666 0.3944193 0.8057968
## 2 2 0.1620125 0.4025078 0.7944065
## 3 3 0.1555915 0.3944509 0.8093665
## 4 4 0.1705968 0.4130337 0.7788997
## 5 5 0.1491317 0.3861758 0.8268758
## 6 6 0.1559972 0.3949648 0.8033898
boots.vld.res[[2]]
## amse rmse auc
## 1 0.1603226 0.4002763 0.7975205
After validating performance metrics from the previous step, we would
like to test if there are some segments within modeling data set where
our final model significantly overestimates or underestimates observed
defaults. This can be achieved using function segment.vld, also
implemented in PDtoolkit package.
segment.vld(model = pd.model,
db = pd.db,
alpha = 0.05)
## $segment.model
## n= 1000
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 1000 160.705 -0.000000000000002700284 *
##
## $segment.testing
## info
## 1 No significant split of residuals.
##
## $segment.rules
## NULL
For more examples check help page ?segment.vld. After decision on the
final model, analysts usually transform model probabilities into scaled
scores and form the rating scale. Next example presents the way to
accomplish these two tasks using PDtoolkit and monobin package.
#extract model predictions
pd.db$probs <- unname(predict(pd.model, type = "response", newdata = pd.db))
#scale probabilities (note: higher score lower probability)
pd.db$score <- round(scaled.score(probs = pd.db$probs, score = 600, odd = 50/1, pdo = 20), 0)
hist(pd.db$score, col = "red", main = "Score distribution", xlab = "score")
#create ratings based on binning algorithm
pd.db$rating <- monobin::sts.bin(x = pd.db$score, y = pd.db$Creditability, y.type = "bina")[[2]]
#summarise rating scale
rs <- pd.db %>%
group_by(rating) %>%
summarise(no = n(),
nb = sum(Creditability),
ng = sum(1 - Creditability)) %>%
mutate(dr = nb / no)
rs
## # A tibble: 6 x 5
## rating no nb ng dr
## <chr> <int> <dbl> <dbl> <dbl>
## 1 01 (-Inf,473) 99 75 24 0.758
## 2 02 [473,482) 72 46 26 0.639
## 3 03 [482,506) 196 83 113 0.423
## 4 04 [506,515) 100 30 70 0.3
## 5 05 [515,547) 267 52 215 0.195
## 6 06 [547,Inf) 266 14 252 0.0526
Often default rate from the modeling data set differs from the so-called
portfolio central tendency (in practice known as long-run average).
Without going deeper into reasons and direction of this difference we
will assume that central tendency for this case study is 27% and
re-calibrate our rating scale. In order to do so, we will use
rs.calibration function from PDtoolkit package. For details on
calibration options check its help page ?rs.calibration.
rs$pd <- rs.calibration(rs = rs,
dr = "dr",
w = "no",
ct = 0.27,
min.pd = 0.05,
method = "log.odds.ab")[[1]]
rs
## # A tibble: 6 x 6
## rating no nb ng dr pd
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 01 (-Inf,473) 99 75 24 0.758 0.740
## 2 02 [473,482) 72 46 26 0.639 0.607
## 3 03 [482,506) 196 83 113 0.423 0.375
## 4 04 [506,515) 100 30 70 0.3 0.252
## 5 05 [515,547) 267 52 215 0.195 0.154
## 6 06 [547,Inf) 266 14 252 0.0526 0.05
#check if calibration is performed as expected
sum(rs$pd * rs$no) / sum(rs$no)
## [1] 0.27
#plot main rating scale characteristics
bp <- barplot(rs$no, names.arg = rs$rating, cex.names = 0.8, las = 2, ylim = c(0, max(rs$no) * 1.2),
main = "Number of observations and default rate per rating")
xlim <- c(floor(min(bp)), ceiling(max(bp)))
par(new = TRUE)
plot(bp, rs$pd, type = "l", col = "red", xaxt = "n", yaxt = "n",
xlab = "", ylab = "", xlim = xlim, ylim = c(0, 1))
text(x = bp, y = rs$pd, label = paste0(round(100 * rs$pd, 2), "%"), col = "red",
cex = 0.6, pos = 3)
#bring PDs to the modeling data set and calculate AUC
pd.db <- merge(pd.db, rs[, c("rating", "pd")], by = "rating", all.x = TRUE)
auc.model(predictions = pd.db$pd, observed = pd.db$Creditability)
## [1] 0.79465
After finalization of PD model development, analysts have to perform
regular validation of the model performance. Model validation is
performed at least on yearly basis on so-called application portfolio.
PDtoolkit has set of functions that are useful for model validation
and next steps present its usage. Before implementation of any of the
validation tests, let’s simulate application portfolio.
set.seed(112233)
app.port <- pd.db[sample(1:nrow(pd.db), 500, replace = FALSE), ]
Now, let’s perform heterogeneity test of the rating scale in order to
check if the model still discriminates well between rating grades. This
validation procedure applies test of two proportions on adjacent rating
grades and check if there is statistical significant difference between
observed default rates. It is always one-tail test, but direction
depends on direction of default rate per rating (increasing or
decreasing).
heterogeneity(app.port = app.port,
def.ind = "Creditability",
rating = "rating",
alpha = 0.05)
## rating no nb dr p.val alpha
## 1 01 (-Inf,473) 47 36 0.76595745 NA 0.05
## 2 02 [473,482) 35 21 0.60000000 0.05319294 0.05
## 3 03 [482,506) 93 37 0.39784946 0.02028868 0.05
## 4 04 [506,515) 53 15 0.28301887 0.08176299 0.05
## 5 05 [515,547) 131 22 0.16793893 0.03889198 0.05
## 6 06 [547,Inf) 141 11 0.07801418 0.01161535 0.05
## res
## 1 <NA>
## 2 H0: DR(02 [473,482)) >= DR(01 (-Inf,473))
## 3 H1: DR(03 [482,506)) < DR(02 [473,482))
## 4 H0: DR(04 [506,515)) >= DR(03 [482,506))
## 5 H1: DR(05 [515,547)) < DR(04 [506,515))
## 6 H1: DR(06 [547,Inf)) < DR(05 [515,547))
Based on the results we can see that there are two pairs of rating
grades where we cannot conclude that observed default rate in the second
rating grade is less than observed default rate in the previous rating
grade. These two pairs are:
- 01 (-Inf,473) - 02 [473,482) and
- 03 [482,506) - 04 [506,515).
Usually in practice the further investigation is performed in order to
better understand source of poor discrimination between these two rating
pairs. Due to the fact that application portfolio in this example
does not have enough observations to perform homogeneity test, the help
page of the ?homogeneity function contains examples, details and
design behind the implemented test. Closely related validation
procedure to heterogeneity and homogeneity is validation of
discriminatory power of the PD model. Namely, analysts usually have to
validate and make sure that discriminatory power of the PD model in use
did not deteriorate significantly from the moment of its development. In
order to test the discriminatory power of the PD model in use, we will
apply dp.testing function for the PDtoolkit package.
#calculate again the AUC for the development sample
auc.dev <- auc.model(predictions = pd.db$pd, observed = pd.db$Creditability)
auc.dev
## [1] 0.79465
#test possible AUC deterioration on the application portfolio
dp.testing(app.port = app.port,
def.ind = "Creditability",
pdc = "pd",
auc.test = auc.dev,
alternative = "less",
alpha = 0.05)
## auc auc.test estimate auc.se test.stat p.val alpha res
## 1 0.7821032 0.79465 -0.01254677 0.0247617 -0.5067005 0.3061825 0.05 H0: AUC >= AUC test
So, from the validation results, we can conclude that slight decrease of
AUC in application portfolio in comparison to development sample cannot
be considered as statistically significant for significance level of 5%.
Beside the discriminatory power, analysts are also interested in testing
the predictive ability of PD model in use. The PDtoolkit package has
pp.testing function that can be used for this purpose. Four tests are
supported in this validation function: three of them refer to testing on
the rating grade level (binomial, Jeffreys, z-score), while one refers
to the complete rating scale (Hosmer-Lemeshow test). The null hypothesis
for the binomial, Jeffreys, and z-score test is that the calibrated PDs are
not underestimated, i.e. they are not significantly lower than the observed
default rate, while for the Hosmer-Lemeshow test is that the calibrated
PD is the true one.
#summarise application portfolio to rating grade level
rs.ap <- app.port %>%
group_by(rating) %>%
summarise(no = n(),
nb = sum(Creditability),
ng = sum(1 - Creditability)) %>%
mutate(dr = nb / no)
#bring calibrated pd as a based for predictive power testing
rs.ap <- merge(rs[, c("rating", "pd")], rs.ap, by = "rating", all.x = TRUE)
rs.ap
## rating pd no nb ng dr
## 1 01 (-Inf,473) 0.7400826 47 36 11 0.76595745
## 2 02 [473,482) 0.6069830 35 21 14 0.60000000
## 3 03 [482,506) 0.3750350 93 37 56 0.39784946
## 4 04 [506,515) 0.2516151 53 15 38 0.28301887
## 5 05 [515,547) 0.1537853 131 22 109 0.16793893
## 6 06 [547,Inf) 0.0500000 141 11 130 0.07801418
#perform predictive power testing
pp.testing(rating.label = rs.ap$rating,
pdc = rs.ap$pd,
no = rs.ap$no,
nb = rs.ap$nb,
alpha = 0.05)
## rating no nb odr pdc alpha binomial binomial.res jeffreys
## 1 01 (-Inf,473) 47 36 0.76595745 0.7400826 0.05 0.4159984 H0: ODR <= PDC 0.35159815
## 2 02 [473,482) 35 21 0.60000000 0.6069830 0.05 0.6059111 H0: ODR <= PDC 0.53854844
## 3 03 [482,506) 93 37 0.39784946 0.3750350 0.05 0.3613128 H0: ODR <= PDC 0.32231637
## 4 04 [506,515) 53 15 0.28301887 0.2516151 0.05 0.3479606 H0: ODR <= PDC 0.29276184
## 5 05 [515,547) 131 22 0.16793893 0.1537853 0.05 0.3621511 H0: ODR <= PDC 0.31876653
## 6 06 [547,Inf) 141 11 0.07801418 0.0500000 0.05 0.0966087 H0: ODR <= PDC 0.07085215
## jeffreys.res zscore zscore.res hosmer.lemeshow hosmer.lemeshow.res
## 1 H0: ODR <= PDC 0.34293973 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 2 H0: ODR <= PDC 0.53370338 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 3 H0: ODR <= PDC 0.32475205 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 4 H0: ODR <= PDC 0.29914802 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 5 H0: ODR <= PDC 0.32669329 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 6 H0: ODR <= PDC 0.06346722 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
As we can see from the results, all three tests rejected the hypothesis
of underestimation of the calibrated PD in comparison to the observed
default rate and Hosmer-Lemeshow test that PD is not the true one.
Sometimes the tests that serve the same purpose can produce opposite
results (lead to opposite conclusions) leaving analysts in dilemma about
overall conclusion. For that purpose we developed the function power
which implements simple framework of power analysis of the statistical
tests. More concretely, this function calculates the power of the four
tests used for predictive ability testing. Assumption of the power
analysis is that observed default rate is the real default rate, so the
ultimate goal of this analysis is to quantify ability of each
statistical test to detect the difference between the calibrated and
real default rate.
power(rating.label = rs.ap$rating,
pdc = rs.ap$pd,
no = rs.ap$no,
nb = rs.ap$nb,
alpha = 0.05,
sim.num = 1000,
seed = 2211)
## $interval.estimator
## rating no nb odr pdc binomial jeffreys zscore
## 1 01 (-Inf,473) 47 36 0.76595745 0.7400826 0.051 0.120 0.120
## 2 02 [473,482) 35 21 0.60000000 0.6069830 0.016 0.053 0.053
## 3 03 [482,506) 93 37 0.39784946 0.3750350 0.087 0.125 0.125
## 4 04 [506,515) 53 15 0.28301887 0.2516151 0.088 0.142 0.142
## 5 05 [515,547) 131 22 0.16793893 0.1537853 0.095 0.095 0.141
## 6 06 [547,Inf) 141 11 0.07801418 0.0500000 0.287 0.391 0.391
##
## $hosmer.lemeshow
## rating
## 1 01 (-Inf,473) + 02 [473,482) + 03 [482,506) + 04 [506,515) + 05 [515,547) + 06 [547,Inf)
## hosmer.lemeshow
## 1 0.268
For this application portfolio characteristic we can conclude that in
average higher power of detecting the real underestimation of calibrated
PD have Jeffreys and z-score in comparison to the binomial test. Since
the Hosmer-Lemeshow test is the only test applied to the complete rating
scale, there is no other test to be compared with. For
additional package features that are implemented but not presented in
the above examples, please, check help(package = "PDtoolkit").
andrija-djurovic/PDtoolkit documentation built on Jan. 29, 2024, 9:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
PDtoolkit 1.2.1
PDtoolkit provides collection of tools for probability of default (PD)
rating model development and validation. Keeping in mind the fact
that model development is highly iterative and repetitive process,
having standardized and automated tools for this purpose is of the
utmost importance for analysts. The main goal of this package is to
cover the most common steps of PD model development. As additional
contribution we attempted to add some functionalities which at the
moment of the package development were not presented in R package
ecosystem in area of credit risk modeling. The package available
functionalities are those that refer to univariate, bivariate,
multivariate analysis, calibration and validation. Along with
accompanied monobin and monobinShiny packages, PDtoolkit provides
functions which are suitable for different data transformation and
modeling tasks such as: imputations, monotonic binning of numeric risk
factors, binning of categorical risk factors, weights of evidence (WoE)
and information value (IV) calculations, WoE coding (replacement of risk
factors’ modalities with WoE values), risk factor clustering, area under
curve (AUC) calculation and others. Beside mentioned features,
package provides set of validation functions for testing homogeneity,
heterogeneity, discriminatory and predictive power of the model.
The following case study shows the usage of PDtoolkit package. The
study is based on publicly available German credit data set which is
available and downloaded from this
link, but also
distributed along with PDtoolkit package under the data frame loans.
Presented examples are simplified, but yet realistic, version of the PD
model development and validation. In the first row, simplification is
assumed in final model selection, master scale creation, calibration and
some validation steps, whereas in practice we usually have highly
iterative process with couple of available options. Anyhow, this does
not diminish the practical value of the presented examples. So, let’s
start with a concrete case study.
First, we will import the libraries needed for the following examples.
library(PDtoolkit)
library(monobin)
library(rpart)
If the packages are not already installed, run the following code before
the library import
install.packages(c("PDtoolkit", "monobin", "rpart")) in order to
install it from the CRAN, while
devtools::install_github("andrija-djurovic/PDtoolkit") to install
development version.
Then, let’s import and inspect the structure of the modeling data set -
loans.
data(loans)
str(loans)
## 'data.frame': 1000 obs. of 21 variables:
## $ Creditability : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Account Balance : chr "1" "1" "2" "1" ...
## $ Duration of Credit (month) : num 18 9 12 12 12 10 8 6 18 24 ...
## $ Payment Status of Previous Credit: chr "4" "4" "2" "4" ...
## $ Purpose : chr "2" "0" "9" "0" ...
## $ Credit Amount : num 1049 2799 841 2122 2171 ...
## $ Value Savings/Stocks : chr "1" "1" "2" "1" ...
## $ Length of current employment : chr "2" "3" "4" "3" ...
## $ Instalment per cent : chr "4" "2" "2" "3" ...
## $ Sex & Marital Status : chr "2" "3" "2" "3" ...
## $ Guarantors : chr "1" "1" "1" "1" ...
## $ Duration in Current address : chr "4" "2" "4" "2" ...
## $ Most valuable available asset : chr "2" "1" "1" "1" ...
## $ Age (years) : num 21 36 23 39 38 48 39 40 65 23 ...
## $ Concurrent Credits : chr "3" "3" "3" "3" ...
## $ Type of apartment : chr "1" "1" "1" "1" ...
## $ No of Credits at this Bank : chr "1" "2" "1" "2" ...
## $ Occupation : chr "3" "3" "2" "2" ...
## $ No of dependents : chr "1" "2" "1" "2" ...
## $ Telephone : chr "1" "1" "1" "1" ...
## $ Foreign Worker : chr "1" "1" "1" "2" ...
Note that variable (column) names have some special characters which is
not usually considered as good practice in real studies, but for this
example we will use it as provided originally in data source.
Additionally, in order to present all functionalities of the PDtoolkit
package, we will slightly modify modeling data set inserting some
missing values for certain variables as given below:
loans$"Credit Amount"[1:10] <- NA
loans$"Age (years)"[30:50] <- NA
loans$"Purpose"[800:820] <- NA
Usually the first step in model development is the univariate analysis.
For this purpose we will use univariate function from the PDtoolkit
package:
univariate(db = loans)
## rf rf.type bin.type bin cnt pct
## 1 Creditability numeric complete cases complete cases 1000 1.000
## 2 Account Balance character complete cases complete cases 1000 1.000
## 3 Duration of Credit (month) numeric complete cases complete cases 1000 1.000
## 4 Payment Status of Previous Credit character complete cases complete cases 1000 1.000
## 5 Purpose character special cases special cases 21 0.021
## 6 Purpose character complete cases complete cases 979 0.979
## 7 Credit Amount numeric special cases special cases 10 0.010
## 8 Credit Amount numeric complete cases complete cases 990 0.990
## 9 Value Savings/Stocks character complete cases complete cases 1000 1.000
## 10 Length of current employment character complete cases complete cases 1000 1.000
## 11 Instalment per cent character complete cases complete cases 1000 1.000
## 12 Sex & Marital Status character complete cases complete cases 1000 1.000
## 13 Guarantors character complete cases complete cases 1000 1.000
## 14 Duration in Current address character complete cases complete cases 1000 1.000
## 15 Most valuable available asset character complete cases complete cases 1000 1.000
## 16 Age (years) numeric special cases special cases 21 0.021
## 17 Age (years) numeric complete cases complete cases 979 0.979
## 18 Concurrent Credits character complete cases complete cases 1000 1.000
## 19 Type of apartment character complete cases complete cases 1000 1.000
## 20 No of Credits at this Bank character complete cases complete cases 1000 1.000
## 21 Occupation character complete cases complete cases 1000 1.000
## 22 No of dependents character complete cases complete cases 1000 1.000
## 23 Telephone character complete cases complete cases 1000 1.000
## 24 Foreign Worker character complete cases complete cases 1000 1.000
## cnt.unique min p1 p5 p25 p50 avg avg.se p75 p95 p99
## 1 2 0 0.00 0.00 0.00 0 0.30000 0.01449863 1.00 1.0 1.00
## 2 4 NA NA NA NA NA NA NA NA NA NA
## 3 33 4 6.00 6.00 12.00 18 20.90300 0.38133320 24.00 48.0 60.00
## 4 5 NA NA NA NA NA NA NA NA NA NA
## 5 1 NA NA NA NA NA NA NA NA NA NA
## 6 10 NA NA NA NA NA NA NA NA NA NA
## 7 1 Inf NA NA NA NA NaN NA NA NA NA
## 8 915 250 424.13 708.45 1371.25 2324 3283.24242 90.03256824 3978.25 9219.7 14194.29
## 9 5 NA NA NA NA NA NA NA NA NA NA
## 10 5 NA NA NA NA NA NA NA NA NA NA
## 11 4 NA NA NA NA NA NA NA NA NA NA
## 12 4 NA NA NA NA NA NA NA NA NA NA
## 13 3 NA NA NA NA NA NA NA NA NA NA
## 14 4 NA NA NA NA NA NA NA NA NA NA
## 15 4 NA NA NA NA NA NA NA NA NA NA
## 16 1 Inf NA NA NA NA NaN NA NA NA NA
## 17 53 19 20.00 22.00 27.00 33 35.52605 0.36315861 42.00 60.0 67.22
## 18 3 NA NA NA NA NA NA NA NA NA NA
## 19 3 NA NA NA NA NA NA NA NA NA NA
## 20 4 NA NA NA NA NA NA NA NA NA NA
## 21 4 NA NA NA NA NA NA NA NA NA NA
## 22 2 NA NA NA NA NA NA NA NA NA NA
## 23 2 NA NA NA NA NA NA NA NA NA NA
## 24 2 NA NA NA NA NA NA NA NA NA NA
## max neg pos cnt.outliers sc.ind
## 1 1 0 300 0 0
## 2 NA NA NA NA 0
## 3 72 0 1000 70 0
## 4 NA NA NA NA 0
## 5 NA NA NA NA 0
## 6 NA NA NA NA 0
## 7 -Inf NA NA 0 0
## 8 18424 0 990 72 0
## 9 NA NA NA NA 0
## 10 NA NA NA NA 0
## 11 NA NA NA NA 0
## 12 NA NA NA NA 0
## 13 NA NA NA NA 0
## 14 NA NA NA NA 0
## 15 NA NA NA NA 0
## 16 -Inf NA NA 0 0
## 17 75 0 979 23 0
## 18 NA NA NA NA 0
## 19 NA NA NA NA 0
## 20 NA NA NA NA 0
## 21 NA NA NA NA 0
## 22 NA NA NA NA 0
## 23 NA NA NA NA 0
## 24 NA NA NA NA 0
Based on the results we can see that univariate treats differently
so-called special and complete cases. For result details and additional
arguments of univariate function check help page ?univariate.
From the structure and univariate results, we conclude that there are
four numeric variables, while the rest are categorical. One of the
numeric variables Creditability presents binary indicator (0/1) of
default status which will serve as our dependent variable for PD model
development. The rest of the variables present potential risk factors
for the PD model. Additionally, we see that some risk factors
(Purpose, Credit Amount, Age (years)) have certain share of
special cases and that for numeric risk factors potential outliers are
identified. Sometimes, analysts want to impute certain values instead of
special cases (usually missing values) as well as for potential
outliers. Without going too much into deeper analysis for this case
study we can perform imputations using two functions from
PDtoolkitpackage: imp.sc for special cases and imp.outliers for
outliers.
#imputation for special cases for risk factors "Credit Amount", "Purpose"
imp.sc.res <- imp.sc(db = loans[, c("Credit Amount", "Purpose")],
method.num = "automatic",
p.val = 0.05)
names(imp.sc.res)
## [1] "db" "report"
#new risk factors with imputed values
head(imp.sc.res[[1]])
## Credit Amount Purpose
## 1 2324 2
## 2 2324 0
## 3 2324 9
## 4 2324 0
## 5 2324 0
## 6 2324 0
#imputation report
imp.sc.res[[2]]
## rf info imputation.method imputed.value imputation.num
## 1 Credit Amount Imputation completed. automatic - median 2324 10
## 2 Purpose Imputation completed. mode NA 21
## imputed.mode
## 1 <NA>
## 2 3
#replace Credit Amount and Purpose with new values
loans[, c("Credit Amount", "Purpose")] <- imp.sc.res[[1]]
colSums(is.na(loans))
## Creditability Account Balance
## 0 0
## Duration of Credit (month) Payment Status of Previous Credit
## 0 0
## Purpose Credit Amount
## 0 0
## Value Savings/Stocks Length of current employment
## 0 0
## Instalment per cent Sex & Marital Status
## 0 0
## Guarantors Duration in Current address
## 0 0
## Most valuable available asset Age (years)
## 0 21
## Concurrent Credits Type of apartment
## 0 0
## No of Credits at this Bank Occupation
## 0 0
## No of dependents Telephone
## 0 0
## Foreign Worker
## 0
#imputation for outliers for risk factor Credit Amount
imp.out.res <- imp.outliers(db = loans[, "Credit Amount", drop = FALSE],
method = "iqr",
range = 1.5)
#new risk factors with imputed values
head(imp.out.res[[1]])
## Credit Amount
## 1 2324
## 2 2324
## 3 2324
## 4 2324
## 5 2324
## 6 2324
#imputation report
imp.out.res[[2]]
## rf info imputation.method imputation.val.upper
## 1 Credit Amount Imputation completed. iqr 7865
## imputation.val.lower imputation.num.upper imputation.num.lower
## 1 250 73 0
#replace Credit Amount with new values
loans[, "Credit Amount"] <- imp.out.res[[1]]
Usually, when building PD models numeric risk factors are discretized,
so we will proceed next with that step. For the purpose of binning the
numeric risk factors, we will use one of the functions (ndr.bin) from
the monobin package. Details about this package can be found
here.
#identify numeric risk factors
num.rf <- sapply(loans, is.numeric)
num.rf <- names(num.rf)[!names(num.rf)%in%"Creditability" & num.rf]
#discretized numeric risk factors using ndr.bin from monobin package
loans[, num.rf] <- sapply(num.rf, function(x)
monobin::ndr.bin(x = loans[, x], y = loans[, "Creditability"])[[2]])
#check loans structure again and confirm binning was sucessful
str(loans)
## 'data.frame': 1000 obs. of 21 variables:
## $ Creditability : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Account Balance : chr "1" "1" "2" "1" ...
## $ Duration of Credit (month) : chr "03 [16,45)" "02 [8,16)" "02 [8,16)" "02 [8,16)" ...
## $ Payment Status of Previous Credit: chr "4" "4" "2" "4" ...
## $ Purpose : chr "2" "0" "9" "0" ...
## $ Credit Amount : chr "01 (-Inf,3914)" "01 (-Inf,3914)" "01 (-Inf,3914)" "01 (-Inf,3914)" ...
## $ Value Savings/Stocks : chr "1" "1" "2" "1" ...
## $ Length of current employment : chr "2" "3" "4" "3" ...
## $ Instalment per cent : chr "4" "2" "2" "3" ...
## $ Sex & Marital Status : chr "2" "3" "2" "3" ...
## $ Guarantors : chr "1" "1" "1" "1" ...
## $ Duration in Current address : chr "4" "2" "4" "2" ...
## $ Most valuable available asset : chr "2" "1" "1" "1" ...
## $ Age (years) : chr "01 (-Inf,26)" "03 [35,Inf)" "01 (-Inf,26)" "03 [35,Inf)" ...
## $ Concurrent Credits : chr "3" "3" "3" "3" ...
## $ Type of apartment : chr "1" "1" "1" "1" ...
## $ No of Credits at this Bank : chr "1" "2" "1" "2" ...
## $ Occupation : chr "3" "3" "2" "2" ...
## $ No of dependents : chr "1" "2" "1" "2" ...
## $ Telephone : chr "1" "1" "1" "1" ...
## $ Foreign Worker : chr "1" "1" "1" "2" ...
As we can see after the binning, all risk factors are categorical. Now
the modeling data set is ready for bivariate analysis. For that purpose
we will use the function bivariate from the PDtoolkitpackage.
bivariate(db = loans, target = "Creditability")
## $results
## rf bin no ng nb pct.o pct.g
## 1 Account Balance 1 274 139 135 0.274 0.198571429
## 2 Account Balance 2 269 164 105 0.269 0.234285714
## 3 Account Balance 3 63 49 14 0.063 0.070000000
## 4 Account Balance 4 394 348 46 0.394 0.497142857
## 5 Duration of Credit (month) 01 (-Inf,8) 87 78 9 0.087 0.111428571
## 6 Duration of Credit (month) 02 [8,16) 344 264 80 0.344 0.377142857
## 7 Duration of Credit (month) 03 [16,45) 499 328 171 0.499 0.468571429
## 8 Duration of Credit (month) 04 [45,Inf) 70 30 40 0.070 0.042857143
## 9 Payment Status of Previous Credit 0 40 15 25 0.040 0.021428571
## 10 Payment Status of Previous Credit 1 49 21 28 0.049 0.030000000
## 11 Payment Status of Previous Credit 2 530 361 169 0.530 0.515714286
## 12 Payment Status of Previous Credit 3 88 60 28 0.088 0.085714286
## 13 Payment Status of Previous Credit 4 293 243 50 0.293 0.347142857
## 14 Purpose 0 224 145 79 0.224 0.207142857
## 15 Purpose 1 102 86 16 0.102 0.122857143
## 16 Purpose 10 12 7 5 0.012 0.010000000
## 17 Purpose 2 177 123 54 0.177 0.175714286
## 18 Purpose 3 300 218 82 0.300 0.311428571
## 19 Purpose 4 12 8 4 0.012 0.011428571
## 20 Purpose 5 21 14 7 0.021 0.020000000
## 21 Purpose 6 49 28 21 0.049 0.040000000
## 22 Purpose 8 9 8 1 0.009 0.011428571
## 23 Purpose 9 94 63 31 0.094 0.090000000
## 24 Credit Amount 01 (-Inf,3914) 739 550 189 0.739 0.785714286
## 25 Credit Amount 02 [3914,Inf) 261 150 111 0.261 0.214285714
## 26 Value Savings/Stocks 1 603 386 217 0.603 0.551428571
## 27 Value Savings/Stocks 2 103 69 34 0.103 0.098571429
## 28 Value Savings/Stocks 3 63 52 11 0.063 0.074285714
## 29 Value Savings/Stocks 4 48 42 6 0.048 0.060000000
## 30 Value Savings/Stocks 5 183 151 32 0.183 0.215714286
## 31 Length of current employment 1 62 39 23 0.062 0.055714286
## 32 Length of current employment 2 172 102 70 0.172 0.145714286
## 33 Length of current employment 3 339 235 104 0.339 0.335714286
## 34 Length of current employment 4 174 135 39 0.174 0.192857143
## 35 Length of current employment 5 253 189 64 0.253 0.270000000
## 36 Instalment per cent 1 136 102 34 0.136 0.145714286
## 37 Instalment per cent 2 231 169 62 0.231 0.241428571
## 38 Instalment per cent 3 157 112 45 0.157 0.160000000
## 39 Instalment per cent 4 476 317 159 0.476 0.452857143
## 40 Sex & Marital Status 1 50 30 20 0.050 0.042857143
## 41 Sex & Marital Status 2 310 201 109 0.310 0.287142857
## 42 Sex & Marital Status 3 548 402 146 0.548 0.574285714
## 43 Sex & Marital Status 4 92 67 25 0.092 0.095714286
## 44 Guarantors 1 907 635 272 0.907 0.907142857
## 45 Guarantors 2 41 23 18 0.041 0.032857143
## 46 Guarantors 3 52 42 10 0.052 0.060000000
## 47 Duration in Current address 1 130 94 36 0.130 0.134285714
## 48 Duration in Current address 2 308 211 97 0.308 0.301428571
## 49 Duration in Current address 3 149 106 43 0.149 0.151428571
## 50 Duration in Current address 4 413 289 124 0.413 0.412857143
## 51 Most valuable available asset 1 282 222 60 0.282 0.317142857
## 52 Most valuable available asset 2 232 161 71 0.232 0.230000000
## 53 Most valuable available asset 3 332 230 102 0.332 0.328571429
## 54 Most valuable available asset 4 154 87 67 0.154 0.124285714
## 55 Age (years) 01 (-Inf,26) 187 108 79 0.187 0.154285714
## 56 Age (years) 02 [26,35) 349 238 111 0.349 0.340000000
## 57 Age (years) 03 [35,Inf) 443 335 108 0.443 0.478571429
## 58 Age (years) SC 21 19 2 0.021 0.027142857
## 59 Concurrent Credits 1 139 82 57 0.139 0.117142857
## 60 Concurrent Credits 2 47 28 19 0.047 0.040000000
## 61 Concurrent Credits 3 814 590 224 0.814 0.842857143
## 62 Type of apartment 1 179 109 70 0.179 0.155714286
## 63 Type of apartment 2 714 528 186 0.714 0.754285714
## 64 Type of apartment 3 107 63 44 0.107 0.090000000
## 65 No of Credits at this Bank 1 633 433 200 0.633 0.618571429
## 66 No of Credits at this Bank 2 333 241 92 0.333 0.344285714
## 67 No of Credits at this Bank 3 28 22 6 0.028 0.031428571
## 68 No of Credits at this Bank 4 6 4 2 0.006 0.005714286
## 69 Occupation 1 22 15 7 0.022 0.021428571
## 70 Occupation 2 200 144 56 0.200 0.205714286
## 71 Occupation 3 630 444 186 0.630 0.634285714
## 72 Occupation 4 148 97 51 0.148 0.138571429
## 73 No of dependents 1 845 591 254 0.845 0.844285714
## 74 No of dependents 2 155 109 46 0.155 0.155714286
## 75 Telephone 1 596 409 187 0.596 0.584285714
## 76 Telephone 2 404 291 113 0.404 0.415714286
## 77 Foreign Worker 1 963 667 296 0.963 0.952857143
## 78 Foreign Worker 2 37 33 4 0.037 0.047142857
## pct.b dr so sg sb dist.g dist.b woe iv.b
## 1 0.450000000 0.4927007 1000 700 300 0.198571429 0.450000000 -0.8180987057 0.2056933888604
## 2 0.350000000 0.3903346 1000 700 300 0.234285714 0.350000000 -0.4013917827 0.0464467634291
## 3 0.046666667 0.2222222 1000 700 300 0.070000000 0.046666667 0.4054651081 0.0094608525225
## 4 0.153333333 0.1167513 1000 700 300 0.497142857 0.153333333 1.1762632229 0.4044104985393
## 5 0.030000000 0.1034483 1000 700 300 0.111428571 0.030000000 1.3121863890 0.1068494631015
## 6 0.266666667 0.2325581 1000 700 300 0.377142857 0.266666667 0.3466246081 0.0382937662266
## 7 0.570000000 0.3426854 1000 700 300 0.468571429 0.570000000 -0.1959478085 0.0198747062913
## 8 0.133333333 0.5714286 1000 700 300 0.042857143 0.133333333 -1.1349799328 0.1026886605902
## 9 0.083333333 0.6250000 1000 700 300 0.021428571 0.083333333 -1.3581234842 0.0840743109238
## 10 0.093333333 0.5714286 1000 700 300 0.030000000 0.093333333 -1.1349799328 0.0718820624131
## 11 0.563333333 0.3188679 1000 700 300 0.515714286 0.563333333 -0.0883186170 0.0042056484275
## 12 0.093333333 0.3181818 1000 700 300 0.085714286 0.093333333 -0.0851578083 0.0006488213969
## 13 0.166666667 0.1706485 1000 700 300 0.347142857 0.166666667 0.7337405775 0.1324227042295
## 14 0.263333333 0.3526786 1000 700 300 0.207142857 0.263333333 -0.2400119704 0.0134863869101
## 15 0.053333333 0.1568627 1000 700 300 0.122857143 0.053333333 0.8344607136 0.0580148877093
## 16 0.016666667 0.4166667 1000 700 300 0.010000000 0.016666667 -0.5108256238 0.0034055041584
## 17 0.180000000 0.3050847 1000 700 300 0.175714286 0.180000000 -0.0240975516 0.0001032752211
## 18 0.273333333 0.2733333 1000 700 300 0.311428571 0.273333333 0.1304779551 0.0049705887671
## 19 0.013333333 0.3333333 1000 700 300 0.011428571 0.013333333 -0.1541506798 0.0002936203425
## 20 0.023333333 0.3333333 1000 700 300 0.020000000 0.023333333 -0.1541506798 0.0005138355994
## 21 0.070000000 0.4285714 1000 700 300 0.040000000 0.070000000 -0.5596157879 0.0167884736381
## 22 0.003333333 0.1111111 1000 700 300 0.011428571 0.003333333 1.2321436813 0.0099744964676
## 23 0.103333333 0.3297872 1000 700 300 0.090000000 0.103333333 -0.1381503385 0.0018420045131
## 24 0.630000000 0.2557510 1000 700 300 0.785714286 0.630000000 0.2208734028 0.0343931441471
## 25 0.370000000 0.4252874 1000 700 300 0.214285714 0.370000000 -0.5461927676 0.0850500166697
## 26 0.723333333 0.3598673 1000 700 300 0.551428571 0.723333333 -0.2713578445 0.0466477056434
## 27 0.113333333 0.3300971 1000 700 300 0.098571429 0.113333333 -0.1395518804 0.0020600515679
## 28 0.036666667 0.1746032 1000 700 300 0.074285714 0.036666667 0.7060505854 0.0265609505935
## 29 0.020000000 0.1250000 1000 700 300 0.060000000 0.020000000 1.0986122887 0.0439444915467
## 30 0.106666667 0.1748634 1000 700 300 0.215714286 0.106666667 0.7042460736 0.0767963575528
## 31 0.076666667 0.3709677 1000 700 300 0.055714286 0.076666667 -0.3192304302 0.0066886375849
## 32 0.233333333 0.4069767 1000 700 300 0.145714286 0.233333333 -0.4708202892 0.0412528253352
## 33 0.346666667 0.3067847 1000 700 300 0.335714286 0.346666667 -0.0321032454 0.0003516069733
## 34 0.130000000 0.2241379 1000 700 300 0.192857143 0.130000000 0.3944152719 0.0247918170922
## 35 0.213333333 0.2529644 1000 700 300 0.270000000 0.213333333 0.2355660713 0.0133487440411
## 36 0.113333333 0.2500000 1000 700 300 0.145714286 0.113333333 0.2513144283 0.0081378005348
## 37 0.206666667 0.2683983 1000 700 300 0.241428571 0.206666667 0.1554664695 0.0054043106061
## 38 0.150000000 0.2866242 1000 700 300 0.160000000 0.150000000 0.0645385211 0.0006453852114
## 39 0.530000000 0.3340336 1000 700 300 0.452857143 0.530000000 -0.1573002887 0.0121345937020
## 40 0.066666667 0.4000000 1000 700 300 0.042857143 0.066666667 -0.4418327523 0.0105198274352
## 41 0.363333333 0.3516129 1000 700 300 0.287142857 0.363333333 -0.2353408346 0.0179307302520
## 42 0.486666667 0.2664234 1000 700 300 0.574285714 0.486666667 0.1655476065 0.0145051236192
## 43 0.083333333 0.2717391 1000 700 300 0.095714286 0.083333333 0.1385189341 0.0017149963274
## 44 0.906666667 0.2998897 1000 700 300 0.907142857 0.906666667 0.0005250722 0.0000002500344
## 45 0.060000000 0.4390244 1000 700 300 0.032857143 0.060000000 -0.6021754024 0.0163447609210
## 46 0.033333333 0.1923077 1000 700 300 0.060000000 0.033333333 0.5877866649 0.0156743110641
## 47 0.120000000 0.2769231 1000 700 300 0.134285714 0.120000000 0.1124779834 0.0016068283347
## 48 0.323333333 0.3149351 1000 700 300 0.301428571 0.323333333 -0.0701507054 0.0015366344996
## 49 0.143333333 0.2885906 1000 700 300 0.151428571 0.143333333 0.0549411180 0.0004447614317
## 50 0.413333333 0.3002421 1000 700 300 0.412857143 0.413333333 -0.0011527379 0.0000005489228
## 51 0.200000000 0.2127660 1000 700 300 0.317142857 0.200000000 0.4610349593 0.0540069523708
## 52 0.236666667 0.3060345 1000 700 300 0.230000000 0.236666667 -0.0285733724 0.0001904891496
## 53 0.340000000 0.3072289 1000 700 300 0.328571429 0.340000000 -0.0341913647 0.0003907584543
## 54 0.223333333 0.4350649 1000 700 300 0.124285714 0.223333333 -0.5860823611 0.0580500624351
## 55 0.263333333 0.4224599 1000 700 300 0.154285714 0.263333333 -0.5346144857 0.0582984367772
## 56 0.370000000 0.3180516 1000 700 300 0.340000000 0.370000000 -0.0845573880 0.0025367216408
## 57 0.360000000 0.2437923 1000 700 300 0.478571429 0.360000000 0.2847014443 0.0337574569686
## 58 0.006666667 0.0952381 1000 700 300 0.027142857 0.006666667 1.4039939382 0.0287484473064
## 59 0.190000000 0.4100719 1000 700 300 0.117142857 0.190000000 -0.4836298810 0.0352358913269
## 60 0.063333333 0.4042553 1000 700 300 0.040000000 0.063333333 -0.4595323294 0.0107224210188
## 61 0.746666667 0.2751843 1000 700 300 0.842857143 0.746666667 0.1211786247 0.0116562296099
## 62 0.233333333 0.3910615 1000 700 300 0.155714286 0.233333333 -0.4044452202 0.0313926528066
## 63 0.620000000 0.2605042 1000 700 300 0.754285714 0.620000000 0.1960517496 0.0263269492328
## 64 0.146666667 0.4112150 1000 700 300 0.090000000 0.146666667 -0.4883527679 0.0276733235151
## 65 0.666666667 0.3159558 1000 700 300 0.618571429 0.666666667 -0.0748774989 0.0036012511391
## 66 0.306666667 0.2762763 1000 700 300 0.344285714 0.306666667 0.1157104961 0.0043529186611
## 67 0.020000000 0.2142857 1000 700 300 0.031428571 0.020000000 0.4519851237 0.0051655442713
## 68 0.006666667 0.3333333 1000 700 300 0.005714286 0.006666667 -0.1541506798 0.0001468101713
## 69 0.023333333 0.3181818 1000 700 300 0.021428571 0.023333333 -0.0851578083 0.0001622053492
## 70 0.186666667 0.2800000 1000 700 300 0.205714286 0.186666667 0.0971637485 0.0018507380658
## 71 0.620000000 0.2952381 1000 700 300 0.634285714 0.620000000 0.0227800283 0.0003254289762
## 72 0.170000000 0.3445946 1000 700 300 0.138571429 0.170000000 -0.2044125146 0.0064243933163
## 73 0.846666667 0.3005917 1000 700 300 0.844285714 0.846666667 -0.0028161100 0.0000067050238
## 74 0.153333333 0.2967742 1000 700 300 0.155714286 0.153333333 0.0154086254 0.0000366872032
## 75 0.623333333 0.3137584 1000 700 300 0.584285714 0.623333333 -0.0646913212 0.0025260420659
## 76 0.376666667 0.2797030 1000 700 300 0.415714286 0.376666667 0.0986375881 0.0038515629628
## 77 0.986666667 0.3073728 1000 700 300 0.952857143 0.986666667 -0.0348672688 0.0011788457545
## 78 0.013333333 0.1081081 1000 700 300 0.047142857 0.013333333 1.2629153400 0.0426985662558
## iv.s auc
## 1 0.66601150335 0.7077690
## 2 0.66601150335 0.7077690
## 3 0.66601150335 0.7077690
## 4 0.66601150335 0.7077690
## 5 0.26770659621 0.6241762
## 6 0.26770659621 0.6241762
## 7 0.26770659621 0.6241762
## 8 0.26770659621 0.6241762
## 9 0.29323354739 0.6268048
## 10 0.29323354739 0.6268048
## 11 0.29323354739 0.6268048
## 12 0.29323354739 0.6268048
## 13 0.29323354739 0.6268048
## 14 0.10939307333 0.5821500
## 15 0.10939307333 0.5821500
## 16 0.10939307333 0.5821500
## 17 0.10939307333 0.5821500
## 18 0.10939307333 0.5821500
## 19 0.10939307333 0.5821500
## 20 0.10939307333 0.5821500
## 21 0.10939307333 0.5821500
## 22 0.10939307333 0.5821500
## 23 0.10939307333 0.5821500
## 24 0.11944316082 0.5778571
## 25 0.11944316082 0.5778571
## 26 0.19600955690 0.5991429
## 27 0.19600955690 0.5991429
## 28 0.19600955690 0.5991429
## 29 0.19600955690 0.5991429
## 30 0.19600955690 0.5991429
## 31 0.08643363103 0.5808190
## 32 0.08643363103 0.5808190
## 33 0.08643363103 0.5808190
## 34 0.08643363103 0.5808190
## 35 0.08643363103 0.5808190
## 36 0.02632209005 0.5433833
## 37 0.02632209005 0.5433833
## 38 0.02632209005 0.5433833
## 39 0.02632209005 0.5433833
## 40 0.04467067763 0.5524238
## 41 0.04467067763 0.5524238
## 42 0.04467067763 0.5524238
## 43 0.04467067763 0.5524238
## 44 0.03201932202 0.5256524
## 45 0.03201932202 0.5256524
## 46 0.03201932202 0.5256524
## 47 0.00358877319 0.5161786
## 48 0.00358877319 0.5161786
## 49 0.00358877319 0.5161786
## 50 0.00358877319 0.5161786
## 51 0.11263826241 0.5853286
## 52 0.11263826241 0.5853286
## 53 0.11263826241 0.5853286
## 54 0.11263826241 0.5853286
## 55 0.12334106269 0.5890381
## 56 0.12334106269 0.5890381
## 57 0.12334106269 0.5890381
## 58 0.12334106269 0.5890381
## 59 0.05761454196 0.5481857
## 60 0.05761454196 0.5481857
## 61 0.05761454196 0.5481857
## 62 0.08539292555 0.5680619
## 63 0.08539292555 0.5680619
## 64 0.08539292555 0.5680619
## 65 0.01326652424 0.5260571
## 66 0.01326652424 0.5260571
## 67 0.01326652424 0.5260571
## 68 0.01326652424 0.5260571
## 69 0.00876276571 0.5214429
## 70 0.00876276571 0.5214429
## 71 0.00876276571 0.5214429
## 72 0.00876276571 0.5214429
## 73 0.00004339223 0.5011905
## 74 0.00004339223 0.5011905
## 75 0.00637760503 0.5195238
## 76 0.00637760503 0.5195238
## 77 0.04387741201 0.5169048
## 78 0.04387741201 0.5169048
##
## $info
## data frame with 0 columns and 0 rows
Based on results of the bivariate analysis we can see that risk factor
Purpose has 10 modalities, that Age (years) has share of special
cases 2.1% and that there are quite some other risk factors with
modality share less than 5%. In order to correct above potential issues,
we can further categorize risk factors Purpose, Foreign Worker
and Age (years) given the following criteria:
- each modality has at least 5% of the observations or minimum 30 observations,
- each modality has at least 1% of defaults or minimum 1 default,
- the maximum number of modalities cannot be greater than 5,
- if special cases exist to be merged with a closest modality (based on default rate).
For this purpose we will use function cat.bin. Additional examples of
cat.bin function usage can be found in its help page ?cat.bin.
Note: Usually in the practice greater attention is paid on this step
and all risk factors are examined closer.
rf <- c("Purpose", "Age (years)", "Foreign Worker")
loans[, rf] <- sapply(rf, function(x)
cat.bin(x = loans[, x],
y = loans[, "Creditability"],
sc = c("SC", NA),
sc.merge = "closest",
min.pct.obs = 0.05,
min.avg.rate = 0.01,
max.groups = 5,
force.trend = "dr")[[2]])
#run bivariate analysis again
bivariate(db = loans[, c(rf, "Creditability")], target = "Creditability")
## $results
## rf bin no ng nb pct.o pct.g pct.b dr so
## 1 Purpose 1 [1,8] 111 94 17 0.111 0.1342857 0.05666667 0.1531532 1000
## 2 Purpose 2 [3] 300 218 82 0.300 0.3114286 0.27333333 0.2733333 1000
## 3 Purpose 3 [2] 177 123 54 0.177 0.1757143 0.18000000 0.3050847 1000
## 4 Purpose 4 [0,4,5,9] 351 230 121 0.351 0.3285714 0.40333333 0.3447293 1000
## 5 Purpose 5 [10,6] 61 35 26 0.061 0.0500000 0.08666667 0.4262295 1000
## 6 Age (years) 1 [03 [35,Inf)] 464 354 110 0.464 0.5057143 0.36666667 0.2370690 1000
## 7 Age (years) 2 [02 [26,35)] 349 238 111 0.349 0.3400000 0.37000000 0.3180516 1000
## 8 Age (years) 3 [01 (-Inf,26)] 187 108 79 0.187 0.1542857 0.26333333 0.4224599 1000
## 9 Foreign Worker 1 [2,1] 1000 700 300 1.000 1.0000000 1.00000000 0.3000000 1000
## sg sb dist.g dist.b woe iv.b iv.s auc
## 1 700 300 0.1342857 0.05666667 0.86278358 0.0669684396 0.1075375 0.5805190
## 2 700 300 0.3114286 0.27333333 0.13047796 0.0049705888 0.1075375 0.5805190
## 3 700 300 0.1757143 0.18000000 -0.02409755 0.0001032752 0.1075375 0.5805190
## 4 700 300 0.3285714 0.40333333 -0.20500910 0.0153268706 0.1075375 0.5805190
## 5 700 300 0.0500000 0.08666667 -0.55004634 0.0201683657 0.1075375 0.5805190
## 6 700 300 0.5057143 0.36666667 0.32151869 0.0447064079 0.1055416 0.5857476
## 7 700 300 0.3400000 0.37000000 -0.08455739 0.0025367216 0.1055416 0.5857476
## 8 700 300 0.1542857 0.26333333 -0.53461449 0.0582984368 0.1055416 0.5857476
## 9 700 300 1.0000000 1.00000000 0.00000000 0.0000000000 0.0000000 NA
##
## $info
## data frame with 0 columns and 0 rows
After additional binning based on above criteria and running again the
bivariate analysis, we see that risk factor Foreign Worker has only
one value and it can be discarded from the further analysis.
loans <- loans[, !names(loans)%in%"Foreign Worker"]
In case when too many risk factors are available, analysts have an
option to perform cluster analysis and shorten the list of potential
candidates for final PD model. For this purpose, PDtoolkit package
contains function rf.clustering. Since the most metrics needed for
distance matrix calculation require numeric values, we will first
replace risk factors’ modalities with Weights of Evidence (WoE) and then
perform clustering.
#replace modalities with WoE
woe.res <- replace.woe(db = loans, target = "Creditability")
head(woe.res[[1]])
## Creditability Account Balance Duration of Credit (month)
## 1 0 -0.8180987 -0.1959478
## 2 0 -0.8180987 0.3466246
## 3 0 -0.4013918 0.3466246
## 4 0 -0.8180987 0.3466246
## 5 0 -0.8180987 0.3466246
## 6 0 -0.8180987 0.3466246
## Payment Status of Previous Credit Purpose Credit Amount Value Savings/Stocks
## 1 0.73374058 -0.02409755 0.2208734 -0.2713578
## 2 0.73374058 -0.20500910 0.2208734 -0.2713578
## 3 -0.08831862 -0.20500910 0.2208734 -0.1395519
## 4 0.73374058 -0.20500910 0.2208734 -0.2713578
## 5 0.73374058 -0.20500910 0.2208734 -0.2713578
## 6 0.73374058 -0.20500910 0.2208734 -0.2713578
## Length of current employment Instalment per cent Sex & Marital Status Guarantors
## 1 -0.47082029 -0.15730029 -0.2353408 0.0005250722
## 2 -0.03210325 0.15546647 0.1655476 0.0005250722
## 3 0.39441527 0.15546647 -0.2353408 0.0005250722
## 4 -0.03210325 0.06453852 0.1655476 0.0005250722
## 5 -0.03210325 -0.15730029 0.1655476 0.0005250722
## 6 -0.47082029 0.25131443 0.1655476 0.0005250722
## Duration in Current address Most valuable available asset Age (years) Concurrent Credits
## 1 -0.001152738 -0.02857337 -0.5346145 0.1211786
## 2 -0.070150705 0.46103496 0.3215187 0.1211786
## 3 -0.001152738 0.46103496 -0.5346145 0.1211786
## 4 -0.070150705 0.46103496 0.3215187 0.1211786
## 5 -0.001152738 -0.02857337 0.3215187 -0.4836299
## 6 0.054941118 0.46103496 0.3215187 0.1211786
## Type of apartment No of Credits at this Bank Occupation No of dependents Telephone
## 1 -0.4044452 -0.0748775 0.02278003 -0.00281611 -0.06469132
## 2 -0.4044452 0.1157105 0.02278003 0.01540863 -0.06469132
## 3 -0.4044452 -0.0748775 0.09716375 -0.00281611 -0.06469132
## 4 -0.4044452 0.1157105 0.09716375 0.01540863 -0.06469132
## 5 0.1960517 0.1157105 0.09716375 -0.00281611 -0.06469132
## 6 -0.4044452 0.1157105 0.09716375 0.01540863 -0.06469132
loans.woe <- woe.res[[1]]
#perform cluster analysis using "common pearson" metric
#number of clusters is set to NA in order to enable automatic determination of optimal number of clusters using elbow method
rf.clustering(db = loans.woe, metric = "common pearson", k = NA)
## rf clusters dist.to.centroid
## 1 Creditability 1 0.0000000
## 2 Account Balance 2 0.0000000
## 4 Duration of Credit (month) 3 0.5788368
## 6 Occupation 3 0.6001961
## 3 Credit Amount 3 0.6892248
## 7 Type of apartment 3 0.7067914
## 5 Most valuable available asset 3 0.8086682
## 8 Payment Status of Previous Credit 4 0.0000000
## 9 Purpose 5 0.0000000
## 10 Value Savings/Stocks 6 0.0000000
## 11 Length of current employment 7 0.0000000
## 12 Instalment per cent 8 0.0000000
## 13 Sex & Marital Status 9 0.0000000
## 14 Guarantors 10 0.0000000
## 15 Duration in Current address 11 0.0000000
## 16 Age (years) 12 0.0000000
## 17 Concurrent Credits 13 0.0000000
## 18 No of Credits at this Bank 14 0.0000000
## 19 No of dependents 15 0.0000000
## 20 Telephone 16 0.0000000
Beside metrics that require numeric inputs, there is x2y metric that
works with both numeric and categorical data. This metric is especially
handy if analyst wants to perform clustering before any binning
procedures and to decrease the number of risk factors. More examples of
the clustering can be found in help page ?rf.clustering, while the
details about x2y metric are presented in this
link.
Due to the fact that for this example we don’t have too many risk factors, we will not use the results of cluster analysis to shorten the list of potential risk factors for the final PD model.
The next step is then to perform multivariate analysis and to select the
final PD model. In order to complete this task we will use function
stepMIV. This function uses the concept of Marginal Information Value
(MIV) in order to select risk factors in the model. Details can be found
on the help page of the function ?stepMIV.
multiv.res <- stepMIV(start.model = Creditability ~ 1,
miv.threshold = 0.02,
m.ch.p.val = 0.05,
coding = "WoE",
db = loans)
## Running iteration: 1
## Running iteration: 2
## Running iteration: 3
## Running iteration: 4
## Running iteration: 5
## Running iteration: 6
## Running iteration: 7
#print the output elements
names(multiv.res)
## [1] "model" "steps" "miv.iter" "warnings" "dev.db"
#print the final model and estimated coefficients
multiv.res$model
##
## Call: glm(formula = as.formula(mod.frm), family = "binomial", data = db)
##
## Coefficients:
## (Intercept) `Account Balance`
## -0.8386 -0.8516
## `Duration of Credit (month)` `Payment Status of Previous Credit`
## -0.9586 -0.7857
## `Value Savings/Stocks` Purpose
## -0.7754 -0.9566
## Guarantors
## -1.1917
##
## Degrees of Freedom: 999 Total (i.e. Null); 993 Residual
## Null Deviance: 1222
## Residual Deviance: 973.6 AIC: 987.6
summary(multiv.res$model)$coefficients
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8385682 0.08070531 -10.390496 2.739256e-25
## `Account Balance` -0.8515981 0.10201758 -8.347562 6.968627e-17
## `Duration of Credit (month)` -0.9586008 0.16294350 -5.883026 4.028332e-09
## `Payment Status of Previous Credit` -0.7856716 0.14969919 -5.248336 1.534793e-07
## `Value Savings/Stocks` -0.7754343 0.19240625 -4.030193 5.573104e-05
## Purpose -0.9565752 0.25273860 -3.784840 1.538075e-04
## Guarantors -1.1916767 0.42403712 -2.810312 4.949344e-03
#store the final model
pd.model <- multiv.res$model
#store the data set used for the model estimation
pd.db <- multiv.res$dev.db
#note that risk factors from the model are coded using WoE due to the fact that WoE coding method is used in stepMIV estimation
head(pd.db)
## Creditability Account Balance Duration of Credit (month)
## 1 0 -0.8180987 -0.1959478
## 2 0 -0.8180987 0.3466246
## 3 0 -0.4013918 0.3466246
## 4 0 -0.8180987 0.3466246
## 5 0 -0.8180987 0.3466246
## 6 0 -0.8180987 0.3466246
## Payment Status of Previous Credit Purpose Credit Amount Value Savings/Stocks
## 1 0.73374058 -0.02409755 01 (-Inf,3914) -0.2713578
## 2 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 3 -0.08831862 -0.20500910 01 (-Inf,3914) -0.1395519
## 4 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 5 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## 6 0.73374058 -0.20500910 01 (-Inf,3914) -0.2713578
## Length of current employment Instalment per cent Sex & Marital Status Guarantors
## 1 2 4 2 0.0005250722
## 2 3 2 3 0.0005250722
## 3 4 2 2 0.0005250722
## 4 3 3 3 0.0005250722
## 5 3 4 3 0.0005250722
## 6 2 1 3 0.0005250722
## Duration in Current address Most valuable available asset Age (years)
## 1 4 2 3 [01 (-Inf,26)]
## 2 2 1 1 [03 [35,Inf)]
## 3 4 1 3 [01 (-Inf,26)]
## 4 2 1 1 [03 [35,Inf)]
## 5 4 2 1 [03 [35,Inf)]
## 6 3 1 1 [03 [35,Inf)]
## Concurrent Credits Type of apartment No of Credits at this Bank Occupation
## 1 3 1 1 3
## 2 3 1 2 3
## 3 3 1 1 2
## 4 3 1 2 2
## 5 1 2 2 2
## 6 3 1 2 2
## No of dependents Telephone
## 1 1 1
## 2 2 1
## 3 1 1
## 4 2 1
## 5 1 1
## 6 2 1
#calculate area under curve (AUC) for the final model
auc.model(predictions = predict(pd.model, type = "response", newdata = pd.db),
observed = pd.db$Creditability)
## [1] 0.7983738
In practice multivariate analysis is highly iterative process, so
analysts usually end up with a few candidate models. After that based on
different statistical and business validations of these models, one is
selected as the final. For the further process we will assume that the
model from the previous step is the final and proceed with its
validation. For that purpose, we will apply k-fold cross-validation and
validation based on bootstrapping. Both methods are implemented in
PDtoolkit package, concretely in functions kfold.vld and
boots.vld.
kfold.vld.res <- kfold.vld(model = pd.model, k = 10, seed = 1984)
kfold.vld.res
## $iter
## k no amse rmse auc
## 1 1 100 0.1630258 0.4037645 0.7712482
## 2 2 100 0.1374731 0.3707736 0.8635211
## 3 3 100 0.1908883 0.4369077 0.7076403
## 4 4 100 0.1179182 0.3433922 0.8726708
## 5 5 100 0.1880557 0.4336539 0.7256897
## 6 6 100 0.1637761 0.4046927 0.8037087
## 7 7 100 0.1400881 0.3742835 0.8508985
## 8 8 100 0.1827844 0.4275329 0.8596527
## 9 9 100 0.1910433 0.4370850 0.7354138
## 10 10 100 0.1617250 0.4021504 0.7536383
##
## $summary
## amse rmse auc
## 1 0.1636778 0.4034236 0.7944082
boots.vld.res <- boots.vld (model = pd.model, B = 10, seed = 1122)
head(boots.vld.res[[1]])
## B amse rmse auc
## 1 1 0.1555666 0.3944193 0.8057968
## 2 2 0.1620125 0.4025078 0.7944065
## 3 3 0.1555915 0.3944509 0.8093665
## 4 4 0.1705968 0.4130337 0.7788997
## 5 5 0.1491317 0.3861758 0.8268758
## 6 6 0.1559972 0.3949648 0.8033898
boots.vld.res[[2]]
## amse rmse auc
## 1 0.1603226 0.4002763 0.7975205
After validating performance metrics from the previous step, we would
like to test if there are some segments within modeling data set where
our final model significantly overestimates or underestimates observed
defaults. This can be achieved using function segment.vld, also
implemented in PDtoolkit package.
segment.vld(model = pd.model,
db = pd.db,
alpha = 0.05)
## $segment.model
## n= 1000
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 1000 160.705 -0.000000000000002700284 *
##
## $segment.testing
## info
## 1 No significant split of residuals.
##
## $segment.rules
## NULL
For more examples check help page ?segment.vld. After decision on the
final model, analysts usually transform model probabilities into scaled
scores and form the rating scale. Next example presents the way to
accomplish these two tasks using PDtoolkit and monobin package.
#extract model predictions
pd.db$probs <- unname(predict(pd.model, type = "response", newdata = pd.db))
#scale probabilities (note: higher score lower probability)
pd.db$score <- round(scaled.score(probs = pd.db$probs, score = 600, odd = 50/1, pdo = 20), 0)
hist(pd.db$score, col = "red", main = "Score distribution", xlab = "score")
#create ratings based on binning algorithm
pd.db$rating <- monobin::sts.bin(x = pd.db$score, y = pd.db$Creditability, y.type = "bina")[[2]]
#summarise rating scale
rs <- pd.db %>%
group_by(rating) %>%
summarise(no = n(),
nb = sum(Creditability),
ng = sum(1 - Creditability)) %>%
mutate(dr = nb / no)
rs
## # A tibble: 6 x 5
## rating no nb ng dr
## <chr> <int> <dbl> <dbl> <dbl>
## 1 01 (-Inf,473) 99 75 24 0.758
## 2 02 [473,482) 72 46 26 0.639
## 3 03 [482,506) 196 83 113 0.423
## 4 04 [506,515) 100 30 70 0.3
## 5 05 [515,547) 267 52 215 0.195
## 6 06 [547,Inf) 266 14 252 0.0526
Often default rate from the modeling data set differs from the so-called
portfolio central tendency (in practice known as long-run average).
Without going deeper into reasons and direction of this difference we
will assume that central tendency for this case study is 27% and
re-calibrate our rating scale. In order to do so, we will use
rs.calibration function from PDtoolkit package. For details on
calibration options check its help page ?rs.calibration.
rs$pd <- rs.calibration(rs = rs,
dr = "dr",
w = "no",
ct = 0.27,
min.pd = 0.05,
method = "log.odds.ab")[[1]]
rs
## # A tibble: 6 x 6
## rating no nb ng dr pd
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 01 (-Inf,473) 99 75 24 0.758 0.740
## 2 02 [473,482) 72 46 26 0.639 0.607
## 3 03 [482,506) 196 83 113 0.423 0.375
## 4 04 [506,515) 100 30 70 0.3 0.252
## 5 05 [515,547) 267 52 215 0.195 0.154
## 6 06 [547,Inf) 266 14 252 0.0526 0.05
#check if calibration is performed as expected
sum(rs$pd * rs$no) / sum(rs$no)
## [1] 0.27
#plot main rating scale characteristics
bp <- barplot(rs$no, names.arg = rs$rating, cex.names = 0.8, las = 2, ylim = c(0, max(rs$no) * 1.2),
main = "Number of observations and default rate per rating")
xlim <- c(floor(min(bp)), ceiling(max(bp)))
par(new = TRUE)
plot(bp, rs$pd, type = "l", col = "red", xaxt = "n", yaxt = "n",
xlab = "", ylab = "", xlim = xlim, ylim = c(0, 1))
text(x = bp, y = rs$pd, label = paste0(round(100 * rs$pd, 2), "%"), col = "red",
cex = 0.6, pos = 3)
#bring PDs to the modeling data set and calculate AUC
pd.db <- merge(pd.db, rs[, c("rating", "pd")], by = "rating", all.x = TRUE)
auc.model(predictions = pd.db$pd, observed = pd.db$Creditability)
## [1] 0.79465
After finalization of PD model development, analysts have to perform
regular validation of the model performance. Model validation is
performed at least on yearly basis on so-called application portfolio.
PDtoolkit has set of functions that are useful for model validation
and next steps present its usage. Before implementation of any of the
validation tests, let’s simulate application portfolio.
set.seed(112233)
app.port <- pd.db[sample(1:nrow(pd.db), 500, replace = FALSE), ]
Now, let’s perform heterogeneity test of the rating scale in order to check if the model still discriminates well between rating grades. This validation procedure applies test of two proportions on adjacent rating grades and check if there is statistical significant difference between observed default rates. It is always one-tail test, but direction depends on direction of default rate per rating (increasing or decreasing).
heterogeneity(app.port = app.port,
def.ind = "Creditability",
rating = "rating",
alpha = 0.05)
## rating no nb dr p.val alpha
## 1 01 (-Inf,473) 47 36 0.76595745 NA 0.05
## 2 02 [473,482) 35 21 0.60000000 0.05319294 0.05
## 3 03 [482,506) 93 37 0.39784946 0.02028868 0.05
## 4 04 [506,515) 53 15 0.28301887 0.08176299 0.05
## 5 05 [515,547) 131 22 0.16793893 0.03889198 0.05
## 6 06 [547,Inf) 141 11 0.07801418 0.01161535 0.05
## res
## 1 <NA>
## 2 H0: DR(02 [473,482)) >= DR(01 (-Inf,473))
## 3 H1: DR(03 [482,506)) < DR(02 [473,482))
## 4 H0: DR(04 [506,515)) >= DR(03 [482,506))
## 5 H1: DR(05 [515,547)) < DR(04 [506,515))
## 6 H1: DR(06 [547,Inf)) < DR(05 [515,547))
Based on the results we can see that there are two pairs of rating grades where we cannot conclude that observed default rate in the second rating grade is less than observed default rate in the previous rating grade. These two pairs are:
- 01 (-Inf,473) - 02 [473,482) and
- 03 [482,506) - 04 [506,515).
Usually in practice the further investigation is performed in order to
better understand source of poor discrimination between these two rating
pairs. Due to the fact that application portfolio in this example
does not have enough observations to perform homogeneity test, the help
page of the ?homogeneity function contains examples, details and
design behind the implemented test. Closely related validation
procedure to heterogeneity and homogeneity is validation of
discriminatory power of the PD model. Namely, analysts usually have to
validate and make sure that discriminatory power of the PD model in use
did not deteriorate significantly from the moment of its development. In
order to test the discriminatory power of the PD model in use, we will
apply dp.testing function for the PDtoolkit package.
#calculate again the AUC for the development sample
auc.dev <- auc.model(predictions = pd.db$pd, observed = pd.db$Creditability)
auc.dev
## [1] 0.79465
#test possible AUC deterioration on the application portfolio
dp.testing(app.port = app.port,
def.ind = "Creditability",
pdc = "pd",
auc.test = auc.dev,
alternative = "less",
alpha = 0.05)
## auc auc.test estimate auc.se test.stat p.val alpha res
## 1 0.7821032 0.79465 -0.01254677 0.0247617 -0.5067005 0.3061825 0.05 H0: AUC >= AUC test
So, from the validation results, we can conclude that slight decrease of
AUC in application portfolio in comparison to development sample cannot
be considered as statistically significant for significance level of 5%.
Beside the discriminatory power, analysts are also interested in testing
the predictive ability of PD model in use. The PDtoolkit package has
pp.testing function that can be used for this purpose. Four tests are
supported in this validation function: three of them refer to testing on
the rating grade level (binomial, Jeffreys, z-score), while one refers
to the complete rating scale (Hosmer-Lemeshow test). The null hypothesis
for the binomial, Jeffreys, and z-score test is that the calibrated PDs are
not underestimated, i.e. they are not significantly lower than the observed
default rate, while for the Hosmer-Lemeshow test is that the calibrated
PD is the true one.
#summarise application portfolio to rating grade level
rs.ap <- app.port %>%
group_by(rating) %>%
summarise(no = n(),
nb = sum(Creditability),
ng = sum(1 - Creditability)) %>%
mutate(dr = nb / no)
#bring calibrated pd as a based for predictive power testing
rs.ap <- merge(rs[, c("rating", "pd")], rs.ap, by = "rating", all.x = TRUE)
rs.ap
## rating pd no nb ng dr
## 1 01 (-Inf,473) 0.7400826 47 36 11 0.76595745
## 2 02 [473,482) 0.6069830 35 21 14 0.60000000
## 3 03 [482,506) 0.3750350 93 37 56 0.39784946
## 4 04 [506,515) 0.2516151 53 15 38 0.28301887
## 5 05 [515,547) 0.1537853 131 22 109 0.16793893
## 6 06 [547,Inf) 0.0500000 141 11 130 0.07801418
#perform predictive power testing
pp.testing(rating.label = rs.ap$rating,
pdc = rs.ap$pd,
no = rs.ap$no,
nb = rs.ap$nb,
alpha = 0.05)
## rating no nb odr pdc alpha binomial binomial.res jeffreys
## 1 01 (-Inf,473) 47 36 0.76595745 0.7400826 0.05 0.4159984 H0: ODR <= PDC 0.35159815
## 2 02 [473,482) 35 21 0.60000000 0.6069830 0.05 0.6059111 H0: ODR <= PDC 0.53854844
## 3 03 [482,506) 93 37 0.39784946 0.3750350 0.05 0.3613128 H0: ODR <= PDC 0.32231637
## 4 04 [506,515) 53 15 0.28301887 0.2516151 0.05 0.3479606 H0: ODR <= PDC 0.29276184
## 5 05 [515,547) 131 22 0.16793893 0.1537853 0.05 0.3621511 H0: ODR <= PDC 0.31876653
## 6 06 [547,Inf) 141 11 0.07801418 0.0500000 0.05 0.0966087 H0: ODR <= PDC 0.07085215
## jeffreys.res zscore zscore.res hosmer.lemeshow hosmer.lemeshow.res
## 1 H0: ODR <= PDC 0.34293973 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 2 H0: ODR <= PDC 0.53370338 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 3 H0: ODR <= PDC 0.32475205 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 4 H0: ODR <= PDC 0.29914802 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 5 H0: ODR <= PDC 0.32669329 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
## 6 H0: ODR <= PDC 0.06346722 H0: ODR <= PDC 0.7851538 H0: PDC is TRUE
As we can see from the results, all three tests rejected the hypothesis
of underestimation of the calibrated PD in comparison to the observed
default rate and Hosmer-Lemeshow test that PD is not the true one.
Sometimes the tests that serve the same purpose can produce opposite
results (lead to opposite conclusions) leaving analysts in dilemma about
overall conclusion. For that purpose we developed the function power
which implements simple framework of power analysis of the statistical
tests. More concretely, this function calculates the power of the four
tests used for predictive ability testing. Assumption of the power
analysis is that observed default rate is the real default rate, so the
ultimate goal of this analysis is to quantify ability of each
statistical test to detect the difference between the calibrated and
real default rate.
power(rating.label = rs.ap$rating,
pdc = rs.ap$pd,
no = rs.ap$no,
nb = rs.ap$nb,
alpha = 0.05,
sim.num = 1000,
seed = 2211)
## $interval.estimator
## rating no nb odr pdc binomial jeffreys zscore
## 1 01 (-Inf,473) 47 36 0.76595745 0.7400826 0.051 0.120 0.120
## 2 02 [473,482) 35 21 0.60000000 0.6069830 0.016 0.053 0.053
## 3 03 [482,506) 93 37 0.39784946 0.3750350 0.087 0.125 0.125
## 4 04 [506,515) 53 15 0.28301887 0.2516151 0.088 0.142 0.142
## 5 05 [515,547) 131 22 0.16793893 0.1537853 0.095 0.095 0.141
## 6 06 [547,Inf) 141 11 0.07801418 0.0500000 0.287 0.391 0.391
##
## $hosmer.lemeshow
## rating
## 1 01 (-Inf,473) + 02 [473,482) + 03 [482,506) + 04 [506,515) + 05 [515,547) + 06 [547,Inf)
## hosmer.lemeshow
## 1 0.268
For this application portfolio characteristic we can conclude that in
average higher power of detecting the real underestimation of calibrated
PD have Jeffreys and z-score in comparison to the binomial test. Since
the Hosmer-Lemeshow test is the only test applied to the complete rating
scale, there is no other test to be compared with. For
additional package features that are implemented but not presented in
the above examples, please, check help(package = "PDtoolkit").
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.