Escuela Politécnica Nacional


Facultad de Ciencias

Modelos de Riesgo

Integrantes:

Fecha: 08 de enero de 2021

Profesora: Dra. Adriana Uquillas

Trabajo N° 1: Modelo de Credit Scoring y Grupos Homogéneos

Introducción

Las entidades bancarias ofrecen préstamos a sus clientes, tomando en cuenta varia información sobre los mismos; así, lo deseable para el banco, es que sus clientes sean puntuales con sus pagos. Es por ello que, los bancos analizan la información para obtener una mejor visión sobre qué personas son puntuales y cuáles no, para así desarrollar metodologías estratégicas y, que el banco no pierda su dinero.

Entonces, para obtener una visión sobre los clientes, se puede realiza lo que se conoce como Modelo de Calificación Crediticia (Credit Scoring), siendo esta una herramienta de gestión de riesgos que evalúa la solvencia crediticia de un solicitante de préstamo, mediante la estimación de su probabilidad de incumplimiento según datos históricos. Utiliza herramientas numéricas para clasificar los casos de pedidos utilizando datos integrados en un valor único que intenta medir el riesgo o la solvencia crediticia.

En el presente trabajo, desarrollaremos el modelo de credit scoring por medio de los datos de la base bank-additional-full, están relacionados con campañas de marketing directo de una institución bancaria portuguesa. Estas campañas de marketing se basaron en llamadas telefónicas. A menudo, se requería más de un contacto con el mismo cliente, para poder acceder si el producto (depósito bancario a plazo) estaría (‘yes’) o no (‘no’) suscrito.

Así, una vez obtenidas las probabilidades de mora de los individuos, por medio del modelo de regresión logística donde se consideraron las variables idiosincráticas, sistémicas y políticas; ahora, nuestra tarea consiste en obtener grupos tales que las probabilidades de mora sean lo más homogéneas posibles.

Objetivos

Metodología

La mayoría de los modelos de calificación crediticia (Credit Scoring) se desarrollan utilizando modelos o metodologías de clasificación comprobados que utilizan “predictores” de campos de datos específicos para estimar la probabilidad de incumplimiento. Dos técnicas comunes son: la regresión logística y los árboles de clasificación o decisión; en el presente trabajo utilizaremos la primera técnica.

Modelo de Regresión Logística

El modelo de regresión logística estima la probabilidad de una variable respuesta binaria/dicotómica basada en una o más variables o características predictoras. Más formalmente, un modelo logístico o logit es aquel en el que las probabilidades logarítmicas de la probabilidad de un evento son una combinación lineal de variables predictoras independientes.

Consideremos una colección de \(n\) variables independientes (explicativas) denotadas por el vector \(X = (X_1,X_2,...,X_n)\). Sea \(P(Y=1 | X) = \pi (x)\) la probabilidad condicional de que el resultado esté presente; el modelo de regresión logística o logit múltiple viene dado por la ecuación:

\[ g(X) = \beta_0 + \beta_1 X_1 + \cdots + \beta_n X_n \] en cuyo caso el modelo de regresión viene dada por la función logística: \[ \pi (X) = \frac{\exp \left\{g(X) \right\} }{1+ \exp \left \{ g(X) \right \} }. \] donde \(\pi\) es la función que toma valores estrictamente entre cero y uno, para todos los números reales.

Si algunas de las variables independientes son cualitativass, variables de escala nominal como etnia, sexo, grupo de tratamiento, etc., no es apropiado incluirlas en el modelo como si fueran variables cuantitativas. Los números utilizados para representar los distintos niveles de estas variables de escala nominal son simplemente identificadores y no tienen significado numérico.

En esta situación, el método de elección es utilizar una colección de variables de diseño o variables ficticias (dummy variables). Supongamos, por ejemplo, que una de las variables independientes es la etnia, que se ha codificado como “mestizo”, “indígena” y “otro”; en este caso, son necesarias dos variables dummy. Una posible estrategia de codificación es que cuando el encuestado es “mestizo”, las dos variables de diseño, \(Z_1\) y \(Z_2\), se fijarían ambas en cero; cuando el encuestado es 'indígena’, \(Z_1\) sería igual a 1 mientras que \(Z_2\) aún sería igual a 0; cuando la raza del encuestado es “otra”, usaríamos \(Z_1 = 0\) y \(Z_2 = 1\).

En general, si una variable cualitativa tiene \(m\) valores posibles, entonces se necesitarán \(m-1\) variables dummy. Esto es cierto ya que, a menos que se indique lo contrario, todos nuestros modelos tienen un término constante. Para ilustrar la notación utilizada para las variables dummy en este trabajo, supongamos que la \(j\_\)ésima variable independiente \(X_j\) tiene \(m_j\) niveles. Las variables dummy \(m_j -1\) se denotarán como \(Z_{jk}\) con \(k \in \{1,2,...,m_j-1\}\) y los coeficientes para estas variables de diseño se denotarán como \(\beta_{jk}\); por lo tanto, el modelo logit con \(n\) variables y la \(j\_\)ésima variable cualitativa con \(m\) niveles sería:

\[ g(X) = \beta_0 + \beta_1 X_1 + \cdots +\sum_{k=1}^{m_j-1} \beta_{jk}X_{jk} + \cdots + \beta_n X_n \]

Al discutir el modelo de regresión logística múltiple, en general, suprimiremos la suma y el doble subíndice necesarios para indicar cuándo se están utilizando las variables dummy. La excepción a esto será la discusión de las estrategias de modelado cuando necesitemos usar el valor específico de los coeficientes para cualquier variable dummy en el modelo.

Software

Para analizar y tratar los datos de la base dada, trabajaremos con el software estadístico R, ya que nos proporciona una gran variedad de paquetes y comandos para desarrollar el trabajo.

Desarrollo

Cargamos todas las librerías que vamos a utilizar en el transcurso del trabajo y, la base de datos:

library(openxlsx)
library(ggplot2)
library(tidyverse)
library(fastDummies)

base_banco <- as.data.frame(read.xlsx("data/bank-additional-full.xlsx"))
c('filas' = nrow(base_banco), 'columnas' =  ncol(base_banco))
##    filas columnas 
##    41188       21

En primer lugar, notamos que la base de datos consta de 41188 objetos, organizada en 21 variables; así, veamos si los nombres de los atributos están en un formato adecuado:

names(base_banco)
##  [1] "age"            "job"            "marital"        "education"     
##  [5] "default"        "housing"        "loan"           "contact"       
##  [9] "month"          "day_of_week"    "duration"       "campaign"      
## [13] "pdays"          "previous"       "poutcome"       "emp.var.rate"  
## [17] "cons.price.idx" "cons.conf.idx"  "euribor3m"      "nr.employed"   
## [21] "y"

Vemos que, en efecto, los nombres de los atributos son adecuados; por tanto, ahora veamos los tipos de datos de cada atributo para poder describir a los mismos y, si es que se lo amerita, cambiar al tipo conveniente.

str(base_banco) 
## 'data.frame':    41188 obs. of  21 variables:
##  $ age           : num  56 57 37 40 56 45 59 41 24 25 ...
##  $ job           : chr  "housemaid" "services" "services" "admin." ...
##  $ marital       : chr  "married" "married" "married" "married" ...
##  $ education     : chr  "basic.4y" "high.school" "high.school" "basic.6y" ...
##  $ default       : chr  "no" "unknown" "no" "no" ...
##  $ housing       : chr  "no" "no" "yes" "no" ...
##  $ loan          : chr  "no" "no" "no" "no" ...
##  $ contact       : chr  "telephone" "telephone" "telephone" "telephone" ...
##  $ month         : chr  "may" "may" "may" "may" ...
##  $ day_of_week   : chr  "mon" "mon" "mon" "mon" ...
##  $ duration      : num  261 149 226 151 307 198 139 217 380 50 ...
##  $ campaign      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ pdays         : num  999 999 999 999 999 999 999 999 999 999 ...
##  $ previous      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ poutcome      : chr  "nonexistent" "nonexistent" "nonexistent" "nonexistent" ...
##  $ emp.var.rate  : chr  "1.1" "1.1" "1.1" "1.1" ...
##  $ cons.price.idx: chr  "93994" "93994" "93994" "93994" ...
##  $ cons.conf.idx : chr  "-36.4" "-36.4" "-36.4" "-36.4" ...
##  $ euribor3m     : chr  "4857" "4857" "4857" "4857" ...
##  $ nr.employed   : chr  "5191" "5191" "5191" "5191" ...
##  $ y             : chr  "no" "no" "no" "no" ...

Vemos que es necesario cambiar los tipos de datos desde la variable emp.var.rate hasta nr.employed:

base_banco[,16:20] <- lapply(base_banco[,16:20], as.numeric)

División de la data

Antes de realizar el análisis exploratorio y hacer algún cambio necesario en las variables, dividamos a la base de datos, tomando el 80\(\%\) de los clientes de la base para el modelamiento y, con el otro 20\(\%\) de infromación realizar la validez del modelo.

library(caret)
set.seed(1984)
training <- createDataPartition(base_banco$y, p = 0.8, list=FALSE)

banco <- base_banco[training,]
banco_test <- base_banco[-training,]

Es decir, ahora, nuestra base de datos tiene la siguiente dimensión:

dim(banco)
## [1] 32951    21

Valores perdidos

Antes de ir con la descripción completa de cada variable, veamos si en la base de datos se presentan valores perdidos:

sum(is.na(banco))
## [1] 0

Como podemos notar, dentro de toda la base de datos, toda información está registrada. Entonces, en la siguiente sección se analizará a cada variable, de donde también, veremos los valores atípicos de aquellas variables que pueden presentar.

Tratamiento de las variables

A continuación, se describe a cada una de las variables existentes donde, además, se mostrará la frecuencia de cada dato que aparece en cada variable. Por otro lado, se estudiarán los valores atípicos de las variables numéricas.

En primer lugar, veamos qué ocurre con la variable respuesta para el modelo.

Variable dependiente (y)

  • Descripción: si el cliente ha suscrito un depósito o no a plazo.

  • Total de valores distintos: se registran dos datos distinto, yes y no.

  • Tabla de frecuencias: mostramos la frecuencia con que aparece esta información

##   Var1  Freq  Porcent
## 1   no 29239 88.73479
## 2  yes  3712 11.26521
  • Gráfico: presentamos el diagrama de frecuencias.

Así, podemos ver que, hay varios clientes a los que no les han dado un depósito; esto puede deberse a la infromación restante que se tiene de los mismos.

Datos bancarios del cliente

1. age

  • Descripción: Edad del cliente, en años.

  • Total de valores distintos:

length(unique(banco[,1]))
## [1] 78
  • Tabla de frecuencias: se muestra la tabla para los primeros 10 valores, ordenados de mayor a menor.
a <- table(banco[,1])
b <- prop.table(a)*100
b <- as.data.frame(b)
names(b)[2] <- 'Porcent'

full_join(as.data.frame(a),b) %>% arrange(desc(Freq)) %>% head(10)
##    Var1 Freq  Porcent
## 1    31 1538 4.667537
## 2    32 1476 4.479378
## 3    33 1473 4.470274
## 4    36 1418 4.303360
## 5    34 1413 4.288185
## 6    35 1398 4.242663
## 7    30 1381 4.191072
## 8    37 1186 3.599284
## 9    29 1157 3.511274
## 10   38 1129 3.426300
  • Gráfico: presentamos el diagrama de frecuencias.
banco %>% ggplot(aes(age)) + geom_bar(col = 'blue', fill='blue') + labs(x = 'Edad', y = 'frecuencia') + scale_x_continuous(breaks = seq(0,95, by = 15)) + theme_light()

  • Datos Atípicos: se presenta el diagrama de caja para detectar los valores atípicos:
banco %>% ggplot(aes(x=0, y=age)) + geom_boxplot() + theme_light() +
labs(y = 'Edad', title = 'Boxplot Edad')

  • Cabe notar que, como se tienen edades altas con poca frecuencia, tratar a esta variable como una numérica podría darnos malas estimaciones, debido a que existen varios valores atípicos; entonces, agrupemos las edades para formar 5 grupos.
grupos <- kmeans(x = banco[, c('age')], centers = 5, nstart = 1)

banco <- banco %>% mutate(gr.edad = grupos$cluster)
banco <- banco %>% mutate(gr.edad = as.factor(gr.edad))

Entonces, veamos el mínimo y máximo de cada grupo formado:

A <- banco %>% filter(gr.edad == 1)
B <- banco %>% filter(gr.edad == 2)
C <- banco %>% filter(gr.edad == 3)
D <- banco %>% filter(gr.edad == 4)
E <- banco %>% filter(gr.edad == 5)

c('minA' = min(A$age), 'maxA' = max(A$age),
  'minB' = min(B$age), 'maxB' = max(B$age),
  'minC' = min(C$age), 'maxC' = max(C$age),
  'minD' = min(D$age), 'maxD' = max(D$age),
  'minE' = min(E$age), 'maxE' = max(E$age))
## minA maxA minB maxB minC maxC minD maxD minE maxE 
##   42   51   28   33   17   27   52   98   34   41

Entonces, vamos a renombrar a dichos grupos de una mejor manera:

library(plyr)
banco[,22] <- revalue(banco[,22], c('1' = 'D', '2' = 'B', '3' = 'A', 
                              '4' = 'E', '5' = 'C'))
detach(package:plyr)

banco[,1] <- banco[,22]
banco <- banco[,-22]

donde se tienen los grupos:

  • A: Menores a 27 años;

  • B: De 28 a 33 años;

  • C: De 34 a 41 años;

  • D: De 42 a 51 años;

  • E: De 52 años en adelante.

banco %>% ggplot(aes(age)) + geom_bar(col = 'blue', fill='blue') + labs(x = 'Rangos de edad', y = 'frecuencia') + theme_light()

2. job

  • Descripción: tipo de trabajo.

  • Total de valores distintos: esta variable consta de 12 categorías

## [1] 12
  • Tabla de frecuencias: se muestra la tabla para los 12 tipos de trabajo, ordenados de mayor a menor.
##             Var1 Freq    Porcent
## 1         admin. 8347 25.3315529
## 2    blue-collar 7431 22.5516676
## 3     technician 5377 16.3181694
## 4       services 3183  9.6597979
## 5     management 2317  7.0316531
## 6        retired 1389  4.2153501
## 7   entrepreneur 1139  3.4566477
## 8  self-employed 1125  3.4141604
## 9      housemaid  845  2.5644138
## 10    unemployed  804  2.4399866
## 11       student  718  2.1789931
## 12       unknown  276  0.8376074
  • Gráfico: presentamos el diagrama de frecuencias.

  • Como podemos notar, esta variable tiene el incoveniente de tener muchas categorías; por ello, necesitamos reducir la cardinalidad de esta. Así, se realiza un análisis de clúster considerando la tasa de morocidad cuando la variable dependiente \(y\) toma valores de
no <- banco %>% filter(y=='no') %>% group_by(job, y) %>% dplyr::summarise(n=n())

yes <- banco %>% group_by(job) %>% dplyr::summarise(n=n())
names(yes)[2] <- 'Total'

cl <- full_join(no, yes, by='job')
cl$NO <- cl$n/cl$Total 
cl$YES <- (cl$Total - cl$n)/cl$Total 
cl <- cl[,-c(2,3,4)]

clust <-data.frame('NO' = cl$NO, 'YES' = cl$YES,
                       row.names = cl$job)

cat <- scale(clust)
d1 <- dist(cat, method = "euclidean")   # Matriz de distancias
hc1 <- hclust(d1, method = "centroid" )     # Cluster Jeráriquico
plot(hc1, cex = 0.6, hang = -1)             # Dendograma y grupos
rect.hclust(hc1, k =4, border = 2:10)

Entonces, se tienen 4 grupos para recategorizar a los tipos de empleos y reducir la cardinalidad de esta variable; denominaremos a los grupos como:

  • T1: housemaid, self-employed, technician, management, unknown, admin y unemployed.

  • T2: blue-collar, entrepreneur y services.

  • T3: student.

  • T4: retired.

library(plyr)
banco[,2] <- revalue(banco[,2], c( 'student' = 'T3', 'retired' = 'T4', 
                                   'housemaid' = 'T1', 'unknown' = 'T1',
                                   'management' = 'T1' , 'entrepreneur' = 'T2',
                                   'blue-collar' = 'T2', 'unemployed' = 'T1', 
                                   'admin.' = 'T1', 'services' = 'T2',
                                   'self-employed' = 'T1', 'technician' = 'T1'))
detach(package:plyr)

3. marital

  • Descripción: estado civil.

  • Total de valores distintos: se tienen 4 categorías

## [1] 4
  • Tabla de frecuencias: se muestra la tabla para los 4 tipos de estado civil.
##       Var1  Freq   Porcent
## 1  married 19970 60.605141
## 2   single  9213 27.959698
## 3 divorced  3705 11.243968
## 4  unknown    63  0.191193
  • Gráfico: presentamos el diagrama de frecuencias.

  • Aunque esta variable no presenta varias categorías, podemos de igual manera agruparlos; utilizando la tasa de morosidad de cada grupo:

4. education

  • Descripción: nivel de educación.

  • Total de valores distintos: se tienen 8 categorías

## [1] 8
  • Tabla de frecuencias: se muestra la tabla para los 8 niveles de educación.
##                  Var1 Freq     Porcent
## 1   university.degree 9760 29.61973840
## 2         high.school 7574 22.98564535
## 3            basic.9y 4860 14.74917301
## 4 professional.course 4192 12.72192043
## 5            basic.4y 3332 10.11198446
## 6            basic.6y 1840  5.58404904
## 7             unknown 1380  4.18803678
## 8          illiterate   13  0.03945252
  • Gráfico: presentamos el diagrama de frecuencias.

  • De igual manera, a esta variable se la puede considerar como una de alta cardinalidad, por ello, vamos a reducir su dimensión.

A estos grupos, los denotaremos como:

  • E1: basic.4y, high.school, professional.course.

  • E2: basic.6y, basic.9y.

  • E3: university,degree, illiterate, unknown.

5. default

  • Descripción: si el cliente tiene o no crédito en incumplimiento.

  • Total de valores distintos: tiene tres categorías

## [1] 3
  • Tabla de frecuencias: se muestra la tabla para los 3 datos.
##      Var1  Freq      Porcent
## 1      no 26063 79.096233802
## 2 unknown  6886 20.897696580
## 3     yes     2  0.006069619
  • Gráfico: presentamos el diagrama de frecuencias.

  • Entonces, esta varible no vamos a recategorizar ya que, como la categoría yes tienen muy poca frecuencia, al unirse con alguna de las otras dos, puede resultar estimaciones inadecuadas para el coeficiente que represente si el cliente tiene crédito de incumpliento o no.

6. housing

  • Descripción: si el cliente tiene o no préstamo de vivienda.

  • Total de valores distintos:

## [1] 3
  • Tabla de frecuencias: se muestra la tabla para los 3 datos.
##      Var1  Freq   Porcent
## 1     yes 17259 52.377773
## 2      no 14888 45.182240
## 3 unknown   804  2.439987
  • Gráfico: presentamos el diagrama de frecuencias.

7. loan

  • Descripción: si el cliente tiene o no préstamo préstamos.

  • Total de valores distintos:

## [1] 3
  • Tabla de frecuencias: se muestra la tabla para los 3 datos.
##      Var1  Freq   Porcent
## 1      no 27149 82.392037
## 2     yes  4998 15.167977
## 3 unknown   804  2.439987
  • Gráfico: presentamos el diagrama de frecuencias.

Último contacto de la campaña actual

8. contact

  • Descripción: tipo de comunicación de contacto.

  • Total de valores distintos: se tienen dos categorías.

## [1] 2
  • Tabla de frecuencias: se muestra la tabla para los 2 datos.
##        Var1  Freq  Porcent
## 1  cellular 20836 63.23329
## 2 telephone 12115 36.76671
  • Gráfico: presentamos el diagrama de frecuencias.

  • Para esta variable no se tiene ningún inconveniente, en cuanto a recategorización.

9. month

  • Descripción: último mes de contacto del año.

  • Total de valores distintos: se tienen 10 meses registrados.

## [1] 10
  • Tabla de frecuencias: se muestra la tabla para los 10 datos.
##    Var1  Freq    Porcent
## 1   may 10995 33.3677278
## 2   jul  5747 17.4410488
## 3   aug  4897 14.8614610
## 4   jun  4282 12.9950533
## 5   nov  3261  9.8965130
## 6   apr  2099  6.3700646
## 7   oct   599  1.8178507
## 8   sep   471  1.4293952
## 9   mar   455  1.3808382
## 10  dec   145  0.4400473
  • Gráfico: presentamos el diagrama de frecuencias.

  • Vemos que en el mes de mayo es donde más se le ha contactado al cliente; siendo diciembre el mes que menos se han comunicado.

10. day_of_week

  • Descripción: último día de contacto de la semana.

  • Total de valores distintos: se tienen registrados 5 días.

## [1] 5
  • Tabla de frecuencias: se muestra la tabla para los 5 días.
##   Var1 Freq  Porcent
## 1  thu 6839 20.75506
## 2  mon 6827 20.71864
## 3  wed 6499 19.72323
## 4  tue 6454 19.58666
## 5  fri 6332 19.21641
  • Gráfico: presentamos el diagrama de frecuencias.

  • Notamos que hay muy poca diferencia entre las frecuencias de los días en que se han contactado.

11. duration

  • Descripción: duración del último contacto, en segundos.

  • Total de valores distintos: el número de datos distintos es de:

## [1] 1466
  • Tabla de frecuencias: se muestra la tabla para los 10 primeros valores, de mayor a menor.
##    Var1 Freq   Porcent
## 1    85  141 0.4279081
## 2    73  135 0.4096993
## 3   135  133 0.4036296
## 4    72  132 0.4005948
## 5    87  132 0.4005948
## 6   136  132 0.4005948
## 7    92  131 0.3975600
## 8   124  129 0.3914904
## 9    93  128 0.3884556
## 10   89  127 0.3854208
  • Gráfico: presentamos el diagrama de frecuencias.

  • Notamos que hay más frecuencia en los tiempos de, aproximandamente, 125 y 250 segundos y, hay muy pocas llamdas que ha durado un tiempo superior a 1000 segundos.

  • Por otro lado, veamos si existen valores atípicos en esta variable, observando el diagrama de cajas:

banco %>% ggplot(aes(x=0, y=duration)) + geom_boxplot() + theme_light() +
labs(y = 'Duración', title = 'Boxplot duration')

En esta variable, tenemos un total de

boxplot.stats(banco$duration)$out %>% length()
## [1] 2374

valores atípicos; entonces, un un criterio para sustituir estos valores es el siguiente:

  • Si los datos se encuentran por debajo del umbral inferior \((Q_1 − 1.5IQR)\), los reemplezamos por la media del conjunto de los datos.

  • Si los datos se encuentra por encima del umbral superior \((Q_3 + 1.5IQR)\), los reemplazamos por la mediana del conjunto de los datos.

reemplazo <- function(x, BI, BS){
  x[x < BI] <- mean(x)
  x[x > BS] <- median(x)
  x
}

BI <- boxplot.stats(banco$duration)$stats[1]
BS <- boxplot.stats(banco$duration)$stats[5]
banco$duration <- reemplazo(banco$duration, BI, BS)

c('mínimo' = min(banco$duration), 'Q1'=quantile(banco$duration,0.25),
  'mediana Q2'=median(banco$duration), 'media'=mean(banco$duration), 'Q3'=quantile(banco$duration,0.75),
  'maximo'=max(banco$duration), 'rango'=max(banco$duration)-min(banco$duration), 'desv.est'=sd(banco$duration))
##     mínimo     Q1.25% mediana Q2      media     Q3.75%     maximo      rango 
##     0.0000   102.0000   180.0000   201.3386   264.5000   642.0000   642.0000 
##   desv.est 
##   135.2699

Otros atributos

12. campaign

  • Descripción: número de contactos realizados durante esa campaña y para ese cliente.

  • Total de valores distintos: el número de datos distintos es de:

## [1] 40
  • Tabla de frecuencias: se muestra la tabla para los 10 primeros valores, de mayor a menor.
##    Var1  Freq    Porcent
## 1     1 14073 42.7088707
## 2     2  8473 25.7139389
## 3     3  4286 13.0071925
## 4     4  2118  6.4277260
## 5     5  1269  3.8511730
## 6     6   791  2.4005341
## 7     7   512  1.5538223
## 8     8   327  0.9923826
## 9     9   224  0.6797973
## 10   10   178  0.5401960
  • Gráfico: presentamos el diagrama de frecuencias.

  • Notamos que hay mayor repetición entre 1 y 2 del total de contactos; así, podemos intuir que, habrán valores atípicos puesto que, sí se han realiado contactos más de 10 veces, siendo estos de poca frecuencia. Veamos el diagrama de caja:
banco %>% ggplot(aes(x=0, y=campaign)) + geom_boxplot() + theme_light() +
labs(y = '# Contactos', title = 'Boxplot campaign')

En esta variable, tenemos un total de

boxplot.stats(banco$campaign)$out %>% length()
## [1] 1941

valores atípicos; entonces los sustituimos utilizando el criterio anterior.

BI <- boxplot.stats(banco$campaign)$stats[1]
BS <- boxplot.stats(banco$campaign)$stats[5]
banco$campaign <- reemplazo(banco$campaign, BI, BS)

c('mínimo' = min(banco$campaign), 'Q1'=quantile(banco$campaign,0.25),
  'mediana Q2'=median(banco$campaign), 'media'=mean(banco$campaign), 'Q3'=quantile(banco$campaign,0.75),
  'maximo'=max(banco$campaign), 'rango'=max(banco$campaign)-min(banco$campaign), 'desv.est'=sd(banco$campaign))
##     mínimo     Q1.25% mediana Q2      media     Q3.75%     maximo      rango 
##   1.000000   1.000000   2.000000   2.043094   3.000000   6.000000   5.000000 
##   desv.est 
##   1.242236

13. pdays

  • Descripción: número de días que pasaron después de que el cliente fue contactado por última vez desde una campaña anterior.

999 significa que el cliente no fue contactado previamente

  • Total de valores distintos: el número de datos distintos es de:
## [1] 26
  • Tabla de frecuencias:
##    Var1  Freq      Porcent
## 1   999 31749 96.352159267
## 2     3   350  1.062183242
## 3     6   336  1.019695912
## 4     4    85  0.257958787
## 5     7    50  0.151740463
## 6     9    49  0.148705654
## 7    12    46  0.139601226
## 8     2    45  0.136566417
## 9     5    37  0.112287943
## 10   10    36  0.109253133
## 11   13    30  0.091044278
## 12   11    24  0.072835422
## 13    1    22  0.066765804
## 14   14    19  0.057661376
## 15   15    17  0.051591757
## 16    8    14  0.042487330
## 17    0    11  0.033382902
## 18   16    10  0.030348093
## 19   17     7  0.021243665
## 20   18     5  0.015174046
## 21   19     2  0.006069619
## 22   21     2  0.006069619
## 23   22     2  0.006069619
## 24   20     1  0.003034809
## 25   25     1  0.003034809
## 26   26     1  0.003034809
  • Gráfico: presentamos el diagrama de frecuencias.

  • Como podemos observar, hay mucha concentración en el valor de 999, por tanto, si consideramos como numérica a esta variable, puede producirnos sesgo en nuestras estimaciones; por ello, a esta variable la vamos a considerar como dicotómica, es decir, 1 si el número de días es 999 y 0 si es otro valor.
banco$pdays <- if_else(banco$pdays == 999, 1, 0)

15. poutcome

  • Descripción: resultado de la campaña de marketing anterior.

  • Total de valores distintos: se tienen 3 categorías:

## [1] 3
  • Tabla de frecuencias:
##          Var1  Freq   Porcent
## 1 nonexistent 28491 86.464751
## 2     failure  3373 10.236412
## 3     success  1087  3.298838
  • Gráfico: presentamos el diagrama de frecuencias.

  • Vemos que en la campaña de marketing anterior, hay más resultados no existentes que los conocidos como éxitos o fracasos.

Contexto social y económico

16. emp.var.rate

  • Descripción: tasa de variación del empleo - indicador trimestral.

  • Total de valores distintos: se tienen 10 datos distintos:

## [1] 10
  • Tabla de frecuencias:
##    Var1  Freq     Porcent
## 1   1.4 12985 39.40699827
## 2  -1.8  7326 22.23301266
## 3   1.1  6223 18.88561804
## 4  -0.1  2937  8.91323480
## 5  -2.9  1316  3.99380899
## 6  -3.4   867  2.63117963
## 7  -1.7   625  1.89675579
## 8  -1.1   527  1.59934448
## 9    -3   137  0.41576887
## 10 -0.2     8  0.02427847
  • Gráfico: presentamos el diagrama de frecuencias.

  • Veamos si existen valores atípicos en esta variable:
banco %>% ggplot(aes(x=0, y=emp.var.rate)) + geom_boxplot() + theme_light() +
labs(x = 'Tasa Variación de empleo', title = 'Boxplot emp.var.rate')

En esta variable, no hay presencia de valores atípicos.

boxplot.stats(banco$emp.var.rate)$out %>% length()
## [1] 0

17. cons.price.idx

  • Descripción: índice de precios al consumidor - indicador mensual.

Mide los cambios en el tiempo del nivel general de los precios, correspondientes al consumo final de bienes y servicios de los hogares de estratos de ingreso

  • Total de valores distintos: el total de datos distintos es:
## [1] 26

Sin embargo, al obtener la tabla de frecuencias, se percató que hay un dato el cual fue mal interpretado al momento de leer la data, pues se tiene un valor de 93.2, de donde, este valor debe ser 93200; entonces, corregimos ese error:

for(i in 1:length(banco$cons.price.idx)){
  if(banco$cons.price.idx[i] == 93.2) {banco$cons.price.idx[i] <- 93200}
}
  • Tabla de frecuencias: mostramos los 10 primeros valores, de mayor a menor:
##     Var1 Freq   Porcent
## 1  93994 6223 18.885618
## 2  93918 5360 16.266578
## 3  92893 4595 13.944949
## 4  93444 4111 12.476101
## 5  94465 3514 10.664320
## 6  93200 2884  8.752390
## 7  93075 1965  5.963400
## 8  92201  601  1.823920
## 9  92963  579  1.757155
## 10 92431  373  1.131984
  • Gráfico: presentamos el diagrama de frecuencias.

  • Así, corroboremos si hay presencia de datos atípicos en esta variable, por medio del diagrama de caja:
banco %>% ggplot(aes(x=0, y=cons.price.idx)) + geom_boxplot() + theme_light() +
labs(x = 'índice precios al consumidor', title = 'Boxplot cons.price.idx')

En consecuencia, podemos ver, que no hay valores atípicos.

18. cons.conf.idx

  • Descripción: índice de confianza del consumidor - indicador mensual.

Está diseñado para evaluar la confianza en general, la salud financiera relativa y el poder adquisitivo del consumidor promedio.

  • Total de valores distintos: el total de datos distintos es:
## [1] 26
  • Tabla de frecuencias: mostramos los 10 primeros valores, de mayor a menor:
##     Var1 Freq   Porcent
## 1  -36.4 6223 18.885618
## 2  -42.7 5360 16.266578
## 3  -46.2 4595 13.944949
## 4  -36.1 4111 12.476101
## 5  -41.8 3514 10.664320
## 6    -42 2884  8.752390
## 7  -47.1 1965  5.963400
## 8  -31.4  601  1.823920
## 9  -40.8  579  1.757155
## 10 -26.9  373  1.131984
  • Gráfico: presentamos el diagrama de frecuencias.

  • Así, corroboremos si hay presencia de datos atípicos en esta variable, por medio del diagrama de caja:
banco %>% ggplot(aes(x=0, y=cons.conf.idx)) + geom_boxplot() + theme_light() +
labs(x = 'índice confianza del consumidor', title = 'Boxplot cons.price.idx')

Del diagrama anterior, vemos que sí existen valores atípicos, para ser exactos, hay:

boxplot.stats(banco$cons.conf.idx)$out %>% length()
## [1] 373
boxplot.stats(banco$cons.conf.idx)$out[1:10]
##  [1] -26.9 -26.9 -26.9 -26.9 -26.9 -26.9 -26.9 -26.9 -26.9 -26.9

de donde podemos ver que es el valor de -26.9; por tanto, reemplazamos estos valores por el de la mediana:

for(i in 1:length(banco$cons.conf.idx)){
  if(banco$cons.conf.idx[i] == -26.9) {banco$cons.conf.idx[i] <- median(banco$cons.conf.idx)}
}

19. euribor3m

  • Descripción: tipo euribor a 3 meses - indicador diario.

euribor es un índice de referencia publicado diariamente que indica el tipo de interés promedio al que un gran número de bancos europeos dicen concederse préstamos a corto plazo

  • Total de valores distintos: el total de datos distintos es:
## [1] 313

Sin embargo, se percató que, al momento de leer esta variable en \(R\), a los dígitos con 3 decimales y cuya parte entera es distinta de cero, los leyó como un número entero de 4 dígitos; entonces, corregimos dichos valores:

for(i in 1:length(banco$euribor3m)){
  if(banco$euribor3m[i] >= 1000) {banco$euribor3m[i] <- banco$euribor3m[i]/1000}
}
  • Tabla de frecuencias: mostramos los 10 primeros valores, de mayor a menor:
##     Var1 Freq  Porcent
## 1  4.857 2274 6.901156
## 2  4.962 2096 6.360960
## 3  4.963 1965 5.963400
## 4  4.961 1517 4.603806
## 5  4.856  986 2.992322
## 6  1.405  942 2.858790
## 7  4.964  942 2.858790
## 8  4.965  839 2.546205
## 9  4.864  835 2.534066
## 10  4.96  826 2.506752
  • Gráfico: presentamos el diagrama de frecuencias.

  • Así, visualicemos el diagrama de caja, para ver si existen valores atípicos:
banco %>% ggplot(aes(x=0, y=euribor3m)) + geom_boxplot() + theme_light() +
labs(x = 'Euribor', title = 'Boxplot euribor')

Entonces, no tenemos muchos inconvenientes con esta variable.

20. nr.employed

  • Descripción: número de empleados - indicador trimestral.

  • Total de valores distintos: el total de datos distintos es:

## [1] 11
  • Tabla de frecuencias: mostramos la freceuncia de los valores, de mayor a menor:
##      Var1  Freq    Porcent
## 1  5228.1 12985 39.4069983
## 2  5099.1  6800 20.6367030
## 3    5191  6223 18.8856180
## 4  5195.8  2937  8.9132348
## 5  5076.2  1316  3.9938090
## 6  5017.5   867  2.6311796
## 7  4991.6   625  1.8967558
## 8  4963.6   527  1.5993445
## 9  5008.7   526  1.5963097
## 10 5023.5   137  0.4157689
  • Gráfico: presentamos el diagrama de frecuencias.

  • Así, visualicemos el diagrama de caja, para ver si existen valores atípicos:
banco %>% ggplot(aes(x=0, y=nr.employed)) + geom_boxplot() + theme_light() +
labs(x = 'Empleados', title = 'Boxplot nr.employed')

Entonces, vemos que no existen valores atícos en esta variable.

boxplot.stats(banco$nr.employed)$out %>% length()
## [1] 0


MODELOS CREDITICIOS

Una vez analizadas y tratadas las variables que más inconvenientes nos hubiesen producido en el modelamiento, realizaremos los modelos de regresión logística en base a las dos filosofías, para así dar resultados respecto a la puntuación de calidad crediticia.


FILOSOFÍA PIT

El sistema de clasificación PIT (Punto en el tiempo) es la clasificación de riesgo de los clientes cambia que cambia de acuerdo a épocas de expansión y recesión de la economía, y de cambios en la política de crédito.

Así, para desarrollar el modelo, tomaremos, en principio, todas las variables presentes en la base de datos, que son las idiosincrácticas, sistémicas y macroeconómicas.

Así, como se había visto, a las variable categóricas las convertimos en dummies; de donde, en \(R\), este proceso es rápido ya que solo utilizamos un comando:

categoricas <- c('age', 'job', 'marital', 'education', 'housing', 'loan', 'poutcome','contact' )

banco_training <- banco %>% dummy_cols(select_columns = categoricas) %>% select(-categoricas, -month, -day_of_week, -default)

Así entonces, corremos el modelo de regresión logística, especificando que la variable dependiente es \(y\) que, para esto, necesitamos cambiar el valor de ‘no’ por 0 y el valor de ‘yes’ por 1:

banco_training$y <- if_else(banco_training$y == 'yes', 1, 0)
banco_training$y <- as.numeric(banco_training$y)

#write.xlsx(banco_training, '00/training.xlsx')

logit_banco <-glm(y~., data = banco_training, family='binomial')

Así, veamos qué variables son explicativas para el modelo:

summary(logit_banco)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = banco_training)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3124  -0.4175  -0.2850  -0.2184   2.9419  
## 
## Coefficients: (10 not defined because of singularities)
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -5.940e+01  1.496e+01  -3.971 7.16e-05 ***
## duration              4.302e-03  1.321e-04  32.563  < 2e-16 ***
## campaign              6.201e-03  1.723e-02   0.360 0.718962    
## pdays                -1.062e+00  2.107e-01  -5.040 4.65e-07 ***
## previous             -2.416e-01  2.225e-01  -1.086 0.277374    
## emp.var.rate         -7.387e-01  6.873e-02 -10.747  < 2e-16 ***
## cons.price.idx        8.671e-04  1.012e-04   8.566  < 2e-16 ***
## cons.conf.idx         3.543e-02  5.605e-03   6.322 2.59e-10 ***
## euribor3m             2.324e-01  7.223e-02   3.218 0.001293 ** 
## nr.employed          -4.512e-03  1.298e-03  -3.477 0.000507 ***
## age_A                 1.158e-01  9.485e-02   1.221 0.222100    
## age_B                 4.499e-02  7.322e-02   0.614 0.538908    
## age_C                -1.822e-01  7.007e-02  -2.600 0.009322 ** 
## age_D                -1.785e-01  7.277e-02  -2.453 0.014160 *  
## age_E                        NA         NA      NA       NA    
## job_T1               -3.089e-01  9.394e-02  -3.288 0.001009 ** 
## job_T2               -5.314e-01  1.004e-01  -5.294 1.20e-07 ***
## job_T3               -1.581e-01  1.428e-01  -1.107 0.268170    
## job_T4                       NA         NA      NA       NA    
## marital_divorced     -1.618e-03  4.249e-01  -0.004 0.996962    
## marital_married       4.292e-02  4.210e-01   0.102 0.918804    
## marital_single        8.348e-02  4.219e-01   0.198 0.843140    
## marital_unknown              NA         NA      NA       NA    
## education_E1         -1.273e-01  4.645e-02  -2.741 0.006127 ** 
## education_E2         -1.389e-01  6.631e-02  -2.095 0.036174 *  
## education_E3                 NA         NA      NA       NA    
## housing_no            2.827e-03  4.101e-02   0.069 0.945050    
## housing_unknown       1.317e-01  1.368e-01   0.963 0.335679    
## housing_yes                  NA         NA      NA       NA    
## loan_no               3.369e-02  5.695e-02   0.591 0.554204    
## loan_unknown                 NA         NA      NA       NA    
## loan_yes                     NA         NA      NA       NA    
## poutcome_failure     -7.591e-01  2.134e-01  -3.557 0.000375 ***
## poutcome_nonexistent         NA         NA      NA       NA    
## poutcome_success             NA         NA      NA       NA    
## contact_cellular      7.912e-01  6.171e-02  12.821  < 2e-16 ***
## contact_telephone            NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23199  on 32950  degrees of freedom
## Residual deviance: 17492  on 32924  degrees of freedom
## AIC: 17546
## 
## Number of Fisher Scoring iterations: 6

Quitamos las variables que no aportan en nada al modelo y luego, vamos retirando una a una, las variables que tienen \(p\_valores\) mayores al \(5\%\); pues, el \(p\_valor\) de cada coeficiente evalúa la hipótesis nula \((H_0)\) de que el coeficiente es igual a cero (no hay efecto), por lo que si un \(p\_valor\) es bajo (usualmente \(< 0,05\)) indica que se puede rechazar la hipótesis nula.

En otras palabras, es probable que un predictor que tenga un \(p\_valor\) bajo sea una adición significativa al modelo porque los cambios en el valor del predictor se relacionan con cambios en la variable de respuesta; por el contrario, un \(p\_valor\) más grande (insignificante) sugiere que los cambios en el predictor no están asociados con cambios en la respuesta.

Así, se obtuvo el siguiente modelo::

banco_training <- banco_training %>% select(-age_E, -job_T4, -marital_unknown, -education_E3, -housing_yes, -loan_unknown, -loan_yes, -poutcome_nonexistent, -poutcome_success, -contact_telephone, -marital_divorced , -housing_no, -campaign, -loan_no, -housing_unknown, -marital_married, -previous, -job_T3, -marital_single, -age_B, -age_A,  -euribor3m, -nr.employed)

logit_banco <-glm(y~., data = banco_training, family='binomial')
summary(logit_banco)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = banco_training)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3128  -0.4172  -0.2839  -0.2196   2.9881  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -1.025e+02  4.836e+00 -21.191  < 2e-16 ***
## duration          4.287e-03  1.318e-04  32.521  < 2e-16 ***
## pdays            -1.285e+00  7.484e-02 -17.175  < 2e-16 ***
## emp.var.rate     -7.512e-01  1.989e-02 -37.764  < 2e-16 ***
## cons.price.idx    1.094e-03  5.167e-05  21.166  < 2e-16 ***
## cons.conf.idx     4.846e-02  4.215e-03  11.496  < 2e-16 ***
## age_C            -2.348e-01  4.966e-02  -4.728 2.27e-06 ***
## age_D            -2.355e-01  5.510e-02  -4.273 1.93e-05 ***
## job_T1           -2.417e-01  6.887e-02  -3.510 0.000449 ***
## job_T2           -4.744e-01  7.675e-02  -6.181 6.37e-10 ***
## education_E1     -1.291e-01  4.615e-02  -2.798 0.005141 ** 
## education_E2     -1.434e-01  6.596e-02  -2.174 0.029716 *  
## poutcome_failure -5.361e-01  6.127e-02  -8.750  < 2e-16 ***
## contact_cellular  8.089e-01  5.898e-02  13.714  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23199  on 32950  degrees of freedom
## Residual deviance: 17511  on 32937  degrees of freedom
## AIC: 17539
## 
## Number of Fisher Scoring iterations: 6


Validez del modelo

Así, una vez generado nuestro modelo crediticio, vamos a utilizar nuestra base de entrenamiento para comparar los resultados y validar adecuamente nuestro modelo. Entonces, tratemos a la base de entrenamiento para poder predecir sus valores:

  ## AGE
for(i in 1:length(banco_test[,1])){
  if(banco_test[i,1] <= 27){banco_test[i,1] = 'A'}
  else if(banco_test[i,1] > 27 & banco_test[i,1] <= 33){banco_test[i,1] = 'B'}
  else if(banco_test[i,1] > 33 & banco_test[i,1] <= 42){banco_test[i,1] = 'C'}
  else if(banco_test[i,1] > 41 & banco_test[i,1] <= 51){banco_test[i,1] = 'D'}
  else if(banco_test[i,1] > 51){banco_test[i,1] = 'E'}
}

library(plyr)
  ## JOB
banco_test[,2] <- revalue(banco_test[,2], c( 'student' = 'T3', 'retired' = 'T4', 
                                   'housemaid' = 'T1', 'unknown' = 'T1',
                                   'management' = 'T1' , 'entrepreneur' = 'T2',
                                   'blue-collar' = 'T2', 'unemployed' = 'T1', 
                                   'admin.' = 'T1', 'services' = 'T2',
                                   'self-employed' = 'T1', 'technician' = 'T1'))

  ## EDUCATION
banco_test[,4] <- revalue(banco_test[,4], c('university.degree' = 'E3', 'illiterate' = 'E3',
                                  'unknown' = 'E3', 'basic.6y' = 'E2', 
                                  'basic.9y' = 'E2', 'basic.4y' = 'E1',
                                  'high.school' = 'E1', 'professional.course' = 'E1'))
detach(package:plyr)

  ## PDAYS
banco_test$pdays <- if_else(banco_test$pdays == 999, 1, 0)

  ## PREVIOUS
banco_test$previous <- if_else(banco_test$previous == 0, 1, 0)

Una vez organizada nuestra base de datos, utilizamos las mismas variables expliciativas del modelo obtenido, para la predicción de las probabilidades de cada individuo:

categoricas <- c('age', 'job', 'marital', 'education', 'housing', 'loan', 'poutcome','contact')

banco_test <- banco_test %>% dummy_cols(select_columns = categoricas) %>% select(-categoricas, -month, -day_of_week, -default)

banco_test <- banco_test %>% select(-age_E, -job_T4, -marital_unknown, -education_E3, -housing_yes, -loan_unknown, -loan_yes, -poutcome_nonexistent, -poutcome_success, -contact_telephone, -marital_divorced , -housing_no, -campaign, -loan_no, -housing_unknown, -marital_married, -previous, -job_T3, -marital_single,  -age_B, -age_A, -euribor3m, -nr.employed)

banco_test$y <- if_else(banco_test$y == 'yes', 1, 0)
banco_test$y <- as.numeric(banco_test$y)


Concepto de Odds o razón de probabilidad

En la regresión lineal simple, se modela el valor de la variable dependiente \(Y\) en función del valor de la variable independiente \(X\). Sin embargo, en la regresión logística, tal como se ha descrito, se modela la probabilidad de que la variable respuesta \(Y\) pertenezca al nivel de referencia 1 en función del valor que adquieran los predictores, mediante el uso del logaritmo de odds.

La transformación de probabilidades a ODDs es monótnoa, si la probabilidad aumenta también lo hacen los ODDs, y viceversa. El rango de valores que pueden tomar los ODDs es de \([0,1]\), pues se recurre a la trasformación logit explicada anteriormente.

odds <- predict(logit_banco, banco_test, type = 'response')
odds[1:20]
##          1          2          3          4          5          6          7 
## 0.05472909 0.03433924 0.01887779 0.06438036 0.03404255 0.02254449 0.02624662 
##          8          9         10         11         12         13         14 
## 0.06237745 0.07263059 0.01873508 0.04990452 0.03076028 0.03394286 0.06704541 
##         15         16         17         18         19         20 
## 0.03663850 0.02954750 0.03047631 0.03860607 0.02666302 0.03067402

Así, una vez predichas las probabilidades de cada individuo de la base de testeo, realizaremos algunas pruebas estadísticas para validar nuestro modelo:


Matriz de Confusión

Este método de clasificación nos ayuda a medir cómo se comporta nuestro modelo cuando se lo aplicamos a base de testeo. La matriz tiene la siguiente estructura:


donde:

  • VP: es la cantidad de positivos que fueron clasificados correctamente como positivos por el modelo.

  • VN: es la cantidad de negativos que fueron clasificados correctamente como negativos por el modelo.

  • FN: es la cantidad de positivos que fueron clasificados incorrectamente como negativos.

  • FP: es la cantidad de negativos que fueron clasificados incorrectamente como positivos.

Ahora, presentamos la matriz de confusión para la base test.Se tiene un total de 501 falsos positivos y 247 falsos negativos. Esto representa una precisión del 90%. Entonces el modelo es adecuado.

test.pred <- ifelse(odds > .5, "1", "0")

matriz<-confusionMatrix(as.factor(banco_test$y), as.factor(test.pred))
matriz$table
##           Reference
## Prediction    0    1
##          0 7062  247
##          1  501  427
matriz$overall[1]
##  Accuracy 
## 0.9091902

A continuación, se presentan los gráficos de las curvas ROC para la base de datos de entrenamiento y la base test. También se presenta el área bajo las curvas, que es el coeficiente AUC y el estadístico de Kolmogorov-Smirnov.

require(ROCR)

train.prob <- predict(logit_banco, type = "response")
train.pred <- ifelse(train.prob > .5, "1", "0")

train.roc <- prediction(train.prob, banco_training$y)
plot(performance(train.roc, "tpr", "fpr"), col = "red", main = "Curva ROC para la base de entrenamiento")
abline(0, 1, lty = 8, col = "blue")

# AUC
train.auc <- performance(train.roc, "auc")
c(AUC=slot(train.auc, "y.values"))
## $AUC
## [1] 0.81948
# KS
ks.train <- performance(train.roc, "tpr", "fpr")
train.ks <- max(attr(ks.train, "y.values")[[1]] - (attr(ks.train, "x.values")[[1]]))
c(KS = train.ks)
##       KS 
## 0.513674
test.prob <- predict(logit_banco, newdata = banco_test, type = "response")
test.pred <- ifelse(test.prob > .5, "1", "0")

test.roc <- prediction(test.prob, banco_test$y)
plot(performance(test.roc, "tpr", "fpr"), col = "red", main = "Curva ROC para la base test")
abline(0, 1, lty = 8, col = "blue")

# AUC
test.auc <- performance(test.roc, "auc")
c(AUC=slot(test.auc, "y.values"))
## $AUC
## [1] 0.8980304
# KS
ks.test <- performance(test.roc, "tpr", "fpr")
test.ks <- max(attr(ks.test, "y.values")[[1]] - (attr(ks.test, "x.values")[[1]]))
c(KS=test.ks)
##        KS 
## 0.7050205

Notemos que el valor AUC de la base de entrenamiento es 0,82 y de la base test es 0,89; por lo tanto, el modelo se considera adecuado. Notemos que el estadístico KS de la base de entrenamiento es 0,51 y el de la base de test es 0,7. En general, el modelo se considera válido.

A continuación, prcedemos a unir los valores predichos en la base original y, vamos a predecir los valores para la base de entrenamiento para, de igual manera, juntarlos en una sola base. Además, vamos a unir la información de la variable indicadora para posteriormente realizar el ciclo económico de cada uno de los grupos.

indicadora <- as.data.frame(read.xlsx("data/bank-additional-full-1.xlsx"))
ind_training <- indicadora[training,]  ## training, índices de la muestra
ind_test <- indicadora[-training,]
  ## TESTING
Testing <- banco_test
Testing$probs <- odds

Testing$month <- ind_test$month1

  ## TRAINING
Training <- banco_training
odds <- predict(logit_banco, banco_training, type = 'response')
Training$probs <- odds

Training$month <- ind_training$month1

banco <- rbind(Training, Testing)


Grupos Homogéneos

Para la formación de los grupos homogéneos vamos a agrupar las probabilidades, utilizando el método de \(k-\)medias; de donde, se formarán 4 grupos:

grupos <- kmeans(x = banco[, c('probs')], centers = 4, nstart = 1)
  ## GH = grupos homogeneos
banco <- banco %>% mutate(GH = grupos$cluster)
banco <- banco %>% mutate(GH = as.factor(GH))

Entonces, veamos los valores máximos y mínimos de cada grupo formado para organizarlos:

C1 <- banco %>% filter(GH == 1)
C2 <- banco %>% filter(GH == 2)
C3 <- banco %>% filter(GH == 3)
C4 <- banco %>% filter(GH == 4)

c('minC1' = min(C1$probs), 'maxC1' = max(C1$probs),
  'minC2' = min(C2$probs), 'maxC2' = max(C2$probs),
  'minC3' = min(C3$probs), 'maxC3' = max(C3$probs),
  'minC4' = min(C4$probs), 'maxC4' = max(C4$probs))
##        minC1        maxC1        minC2        maxC2        minC3        maxC3 
## 2.832617e-01 5.671461e-01 2.220446e-16 1.057109e-01 1.057462e-01 2.832166e-01 
##        minC4        maxC4 
## 5.674784e-01 9.999724e-01

Entonces, notamos que el primer grupo es el \(C_2\), el segundo grupo es \(C_3\), el tercer grupo es \(C_1\) y último grupo, los que tiene altas probabilidades de mora, es \(C_4\).

Para evitarnos confusiones, renombramos estos grupos respecto a su orden. Así, Una vez especificados los grupos, veamos las pruebas de diferencias de medias para cada uno; además, para asegurarnos que los grupos no son homogéneos entre sí. Entonces, para aplicar una de las pruebas estadísticas que contrastan la hipótesis de igualdad de medias, veamos si la varianza de los grupos, en cuanto a las probabilidades, son iguales o no:

Para esto, haremos uso de la prueba de Barlet, la cual contrasta la hipótesis nula de si todas las varianzas de los grupos son iguales:

G1 <- C2; G2 <- C3; G3 <- C1; G4 <- C4

G1$GH <- as.character(G1$GH)
G2$GH <- as.character(G2$GH)
G1$GH <- as.character(G1$GH)
G2$GH <- as.character(G2$GH)

variables <- c(G1$probs, G2$probs)
categorias <- c(G1$GH, G2$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 7815.2, df = 1, p-value < 2.2e-16

Podemos ver que, como el \(p\_valor\) es menor a un \(\alpha = 0.05\), se rechaza la hipótesis nula y, en consecuencia, las varianzas entre el grupo 1 y grupo 2 son diferentes. De manera análoga, se contrasta para los otros grupos, dándonos para todos los casos que el \(p\_valor\) es menor al nivel de significancia \(\alpha\).

variables <- c(G1$probs, G3$probs)
categorias <- c(G1$GH, G3$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 13892, df = 1, p-value < 2.2e-16
variables <- c(G1$probs, G4$probs)
categorias <- c(G1$GH, G4$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 18325, df = 1, p-value < 2.2e-16
variables <- c(G2$probs, G3$probs)
categorias <- c(G2$GH, G3$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 1001.9, df = 1, p-value < 2.2e-16
variables <- c(G2$probs, G4$probs)
categorias <- c(G2$GH, G4$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 2241.9, df = 1, p-value < 2.2e-16
variables <- c(G3$probs, G4$probs)
categorias <- c(G3$GH, G4$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 228.06, df = 1, p-value < 2.2e-16

Con esto, realizamos la prueba de Games y Howell, para la comparación de medias entre los grupos:

library(userfriendlyscience)
banco$GH <- as.factor(banco$GH)
GH <- oneway(banco$GH, y = banco$probs, posthoc = 'games-howell')
GH
## ### Oneway Anova for y=probs and x=GH (groups: 1, 2, 3, 4)
## 
## Omega squared: 95% CI = [.94; .94], point estimate = .94
## Eta Squared: 95% CI = [.94; .94], point estimate = .94
## 
##                                      SS    Df     MS         F     p
## Between groups (error + effect) 1085.22     3 361.74 207730.75 <.001
## Within groups (error only)        71.72 41184      0                
## 
## 
## ### Post hoc test: games-howell
## 
##      diff ci.lo ci.hi      t      df     p
## 2-1 -0.36 -0.36 -0.35 239.87 2910.65 <.001
## 3-1 -0.23 -0.23 -0.22 143.81 3805.82 <.001
## 4-1  0.34  0.33  0.35 114.24 2897.03 <.001
## 3-2  0.13  0.13  0.13 209.97 7752.35 <.001
## 4-2  0.69  0.69  0.70 268.48 1766.89 <.001
## 4-3  0.57  0.56  0.57 214.31 1942.44 <.001

En consecuencia, como los \(p\_valores\) de la comparación de cada dos grupos son menores al \(5\%\), concluimos que los grupos tienen medias distintas; es decir, los 4 grupos formados son homogéneos pero heterogéneos entre sí.

Ahora entonces, por medio del uso de la variable indicadora month, tenemos 8 tiempos para realizar el ciclo económico; de donde, por cada grupo y por cada mes, se calculará el procentaje de clientes malos \(y = no = 0\) en ese grupo:

banco$GH <- as.character(banco$GH)
library(plyr)
banco$GH <- revalue(banco$GH, c('1' = 'G3', '2' = 'G1', '3' = 'G2', 
                              '4' = 'G4'))
detach(package:plyr)

for(i in 1:length(banco$month)){
  if(banco$month[i] == 6 | banco$month[i] == 7 | banco$month[i] == 8){
    banco$month[i] <- 5
  }else if(banco$month[i] == 9){
    banco$month[i] <- 6
  }else if(banco$month[i] == 10){
    banco$month[i] <- 7
  }else if(banco$month[i] == 11 | banco$month[i] == 12 | banco$month[i] == 13 |
           banco$month[i] == 14 | banco$month[i] == 15 | banco$month[i] == 16 |
           banco$month[i] == 17 | banco$month[i] == 18 | banco$month[i] == 19 |
           banco$month[i] == 20 | banco$month[i] == 21 | banco$month[i] == 22 |
           banco$month[i] == 23 | banco$month[i] == 24 | banco$month[i] == 25 |
           banco$month[i] == 26){
    banco$month[i] <- 8
  }
}

no <- banco %>% filter(y == 0) %>% group_by(GH, month) %>% dplyr::summarise(tmalos = n())

yes <- banco %>% filter(y == 1) %>% group_by(GH, month) %>% dplyr::summarise(tbuenos = n())

mes <- full_join(no, yes, by = c('GH', 'month'))
mes$tmalos <- mes$tmalos/(mes$tmalos + mes$tbuenos)

Así también, asignamos la PD de cada grupo:

PD <- mes %>% group_by(GH) %>% dplyr::summarise(PD = mean(tmalos, na.rm = TRUE))
PD
## # A tibble: 4 x 2
##   GH       PD
##   <chr> <dbl>
## 1 G1    0.963
## 2 G2    0.826
## 3 G3    0.589
## 4 G4    0.386
ggplot(data=mes) + geom_line(aes(x = month, y = tmalos, color = GH), size=1) +
  labs(x = "Periodo", y = "PD",
       title = "Ciclo Económico") +
  theme_bw()+
  theme(legend.background = element_rect()) +
  scale_color_hue(name = 'GH')

mes <- full_join(no, yes, by = c('GH', 'month'))
mes$total <- mes$tmalos + mes$tbuenos
bdd <- mes %>% group_by(month) %>% mutate(porc = round(100*(total/sum(total,na.rm=T)),2))

ggplot(data=bdd,
       aes(x = factor(month, level= c("1","2","3","4","5","6","7","8")),
           y = porc, fill = GH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
  xlab("Meses") + ylab("Porcentaje") +
  labs(title = "Porcentaje de clientes por grupos en cada mes",
       fill= "Grupos") +
  scale_fill_brewer(palette="Spectral") 
## Warning: Removed 1 rows containing missing values (position_stack).

De los gráficos podemos ver que, en el grupo G1, no existe información para el mes 5.


FILOSOFÍA TTC

El sistema de clasificación TTC (A través del ciclo) es un sistema donde la clasificación de los clientes en su mayoría no se ve afectado por las condiciones macroeconómicas del país

Para realizar este modelo utilizando la filsofía TTC, vamos tomar las variables idiosincráticas, es decir, tomaremos las variables: age, job, marital, education, housing, loan, contact y duration

Previamente se había guardado la base de entrenamiento y testeo con los datos ya tratados, para desarrollar el modelo de credit scoring usando la filosofía TTC. Así, en primer lugar, formulamos el modelo utilizando las variables antes mencionadas siguiendo el mismo procedimiento para el modelo anterior.

trainingTTC <- trainingTTC[,c(1,2,3,4,6,7,8,11,21)]
categoricas <- c('age', 'job', 'marital', 'education', 'housing', 'loan', 'contact')

banco_training <- trainingTTC %>% dummy_cols(select_columns = categoricas) %>% select(-categoricas)

banco_training$y <- if_else(banco_training$y == 'yes', 1, 0)
banco_training$y <- as.numeric(banco_training$y)

logit_banco <-glm(y~., data = banco_training, family='binomial')

Así, veamos qué variables son explicativas para el modelo:

summary(logit_banco)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = banco_training)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7944  -0.5000  -0.3773  -0.2520   2.8271  
## 
## Coefficients: (8 not defined because of singularities)
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -2.7655230  0.4100562  -6.744 1.54e-11 ***
## duration           0.0042704  0.0001186  35.993  < 2e-16 ***
## age_A              0.1251227  0.0848045   1.475  0.14010    
## age_B             -0.1833562  0.0654346  -2.802  0.00508 ** 
## age_C             -0.3620118  0.0626473  -5.779 7.53e-09 ***
## age_D             -0.4407489  0.0653480  -6.745 1.53e-11 ***
## age_E                     NA         NA      NA       NA    
## job_T1            -0.8016449  0.0814685  -9.840  < 2e-16 ***
## job_T2            -1.1667824  0.0870962 -13.396  < 2e-16 ***
## job_T3             0.0390040  0.1257665   0.310  0.75646    
## job_T4                    NA         NA      NA       NA    
## marital_divorced  -0.0126858  0.4021515  -0.032  0.97483    
## marital_married    0.0832394  0.3987470   0.209  0.83464    
## marital_single     0.2094749  0.3995897   0.524  0.60012    
## marital_unknown           NA         NA      NA       NA    
## education_E1      -0.2083156  0.0419486  -4.966 6.84e-07 ***
## education_E2      -0.3001179  0.0614558  -4.883 1.04e-06 ***
## education_E3              NA         NA      NA       NA    
## housing_no        -0.0148573  0.0375335  -0.396  0.69222    
## housing_unknown    0.1812860  0.1264645   1.433  0.15172    
## housing_yes               NA         NA      NA       NA    
## loan_no            0.0552194  0.0519925   1.062  0.28821    
## loan_unknown              NA         NA      NA       NA    
## loan_yes                  NA         NA      NA       NA    
## contact_cellular   1.0333090  0.0467471  22.104  < 2e-16 ***
## contact_telephone         NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23199  on 32950  degrees of freedom
## Residual deviance: 20476  on 32933  degrees of freedom
## AIC: 20512
## 
## Number of Fisher Scoring iterations: 5

De igual manera, vamos retirando una a una las variables que no aportan al modelo:

banco_training <- banco_training %>% select(-age_E, -job_T4, -marital_unknown, -education_E3, -housing_yes, -loan_unknown, -loan_yes, -contact_telephone, -marital_divorced, -job_T3, -housing_no, -loan_no, -housing_unknown, -marital_married, -age_A)

logit_banco <-glm(y~., data = banco_training, family='binomial')
summary(logit_banco)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = banco_training)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7478  -0.4996  -0.3782  -0.2521   2.8363  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.6220007  0.0787946 -33.276  < 2e-16 ***
## duration          0.0042695  0.0001186  36.008  < 2e-16 ***
## age_B            -0.2329333  0.0531352  -4.384 1.17e-05 ***
## age_C            -0.4052630  0.0529558  -7.653 1.97e-14 ***
## age_D            -0.4799909  0.0578124  -8.303  < 2e-16 ***
## job_T1           -0.8058047  0.0632386 -12.742  < 2e-16 ***
## job_T2           -1.1639380  0.0701931 -16.582  < 2e-16 ***
## marital_single    0.1807164  0.0408463   4.424 9.68e-06 ***
## education_E1     -0.2071977  0.0418668  -4.949 7.46e-07 ***
## education_E2     -0.2911222  0.0613017  -4.749 2.04e-06 ***
## contact_cellular  1.0342564  0.0466346  22.178  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 23199  on 32950  degrees of freedom
## Residual deviance: 20484  on 32940  degrees of freedom
## AIC: 20506
## 
## Number of Fisher Scoring iterations: 5

Validez del Modelo

Una vez obtenido el modelo con la base de entranimiento, ahora, vamos a validar el mismo, utilizando la base de testeo guardad previamente.

testingTTC <- testingTTC[,c(1,2,3,4,6,7,8,11,21)]

categoricas <- c('age', 'job', 'marital', 'education', 'housing', 'loan', 'contact')

banco_test <- testingTTC %>% dummy_cols(select_columns = categoricas) %>% select(-categoricas)

banco_test <- banco_test %>% select(-age_E, -job_T4, -marital_unknown, -education_E3, -housing_yes, -loan_unknown, -loan_yes, -contact_telephone, -marital_divorced, -job_T3, -housing_no, -loan_no, -housing_unknown, -marital_married, -age_A)

banco_test$y <- if_else(banco_test$y == 'yes', 1, 0)
banco_test$y <- as.numeric(banco_test$y)
odds <- predict(logit_banco, banco_test, type = 'response')
odds[1:20]
##          1          2          3          4          5          6          7 
## 0.06401966 0.04558393 0.02662430 0.08439756 0.04265887 0.02399327 0.03197234 
##          8          9         10         11         12         13         14 
## 0.10355166 0.08150926 0.01881222 0.06580304 0.02668990 0.03840487 0.07619754 
##         15         16         17         18         19         20 
## 0.04546208 0.02743823 0.03436322 0.04682742 0.03831409 0.02848246

Así, una vez predichas las probabilidades de cada individuo de la base de testeo, realizaremos algunas pruebas estadísticas para validar nuestro modelo:


Matriz de Confusión

Corremos la matriz de confusión y vemos la precisión que tiene el modelo respecto a la base de datos.

test.pred <- ifelse(odds > .5, "1", "0")

matriz<-confusionMatrix(as.factor(banco_test$y), as.factor(test.pred))
matriz$table
##           Reference
## Prediction    0    1
##          0 7089  220
##          1  653  275
matriz$overall[1]
##  Accuracy 
## 0.8940148

Vemos que se tiene un total de 653 falsos positivos y 220 falsos negativos. Además, se tiene que representa una precisión del 90%. Entonces el modelo es adecuado.

require(ROCR)

train.prob <- predict(logit_banco, banco_training, type = "response")
train.pred <- ifelse(train.prob > .5, "1", "0")

train.roc <- prediction(train.prob, banco_training$y)
plot(performance(train.roc, "tpr", "fpr"), col = "red", main = "Curva ROC para la data de entrenamiento")
abline(0, 1, lty = 8, col = "blue")

# AUC
train.auc <- performance(train.roc, "auc")
c(AUC=slot(train.auc, "y.values"))
## $AUC
## [1] 0.7506053
# KS
ks.train <- performance(train.roc, "tpr", "fpr")
train.ks <- max(attr(ks.train, "y.values")[[1]] - (attr(ks.train, "x.values")[[1]]))
c(KS = train.ks)
##        KS 
## 0.3966375

De los valores AUC y KS obtenidos para la base de entrenamiento, podemos concluir que el modelo es adecuado; como también, si se usa la base de testeo:

test.prob <- predict(logit_banco, newdata = banco_test, type = "response")
test.pred <- ifelse(test.prob > .5, "1", "0")

test.roc <- prediction(test.prob, banco_test$y)
plot(performance(test.roc, "tpr", "fpr"), col = "red", main = "Curva ROC para la data de testeo")
abline(0, 1, lty = 8, col = "blue")

# AUC
test.auc <- performance(test.roc, "auc")
c(AUC=slot(test.auc, "y.values"))
## $AUC
## [1] 0.8626456
# KS
ks.test <- performance(test.roc, "tpr", "fpr")
test.ks <- max(attr(ks.test, "y.values")[[1]] - (attr(ks.test, "x.values")[[1]]))
c(KS=test.ks)
##        KS 
## 0.5815846

A continuación, procedemos a unir los valores predichos en la base original y, vamos a predecir los valores para la base de entrenamiento para, de igual manera, juntarlos en una sola base. Además, vamos a unir la información de la variable indicadora para posteriormente realizar el ciclo económico de cada uno de los grupos.

  ## TESTING
Testing <- banco_test
Testing$probs <- odds

Testing$month <- ind_test$month1

  ## TRAINING
Training <- banco_training
odds <- predict(logit_banco, banco_training, type = 'response')
Training$probs <- odds

Training$month <- ind_training$month1

banco <- rbind(Training, Testing)
str(banco)
## 'data.frame':    41188 obs. of  13 variables:
##  $ duration        : num  261 149 226 151 198 217 380 55 222 137 ...
##  $ y               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ age_B           : int  0 0 0 0 0 0 0 0 0 1 ...
##  $ age_C           : int  0 0 1 1 0 1 0 1 0 0 ...
##  $ age_D           : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ job_T1          : int  1 0 0 1 0 0 1 0 0 0 ...
##  $ job_T2          : int  0 1 1 0 1 1 0 1 1 1 ...
##  $ marital_single  : int  0 0 0 0 0 0 1 0 1 1 ...
##  $ education_E1    : int  1 1 1 0 0 0 1 0 1 1 ...
##  $ education_E2    : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ contact_cellular: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ probs           : num  0.0744 0.0337 0.0313 0.0299 0.0239 ...
##  $ month           : num  1 1 1 1 1 1 1 1 1 1 ...


Grupos Homogéneos

Para la formación de los grupos homogéneos vamos a agrupar las probabilidades, utilizando el método de \(k-\)medias; de donde, se formarán 4 grupos:

grupos <- kmeans(x = banco[, c('probs')], centers = 4, nstart = 1)
  ## GH = grupos homogeneos
banco <- banco %>% mutate(GH = grupos$cluster)
banco <- banco %>% mutate(GH = as.factor(GH))

Entonces, veamos los valores máximos y mínimos de cada grupo formado para organizarlos:

C1 <- banco %>% filter(GH == 1)
C2 <- banco %>% filter(GH == 2)
C3 <- banco %>% filter(GH == 3)
C4 <- banco %>% filter(GH == 4)

c('minC1' = min(C1$probs), 'maxC1' = max(C1$probs),
  'minC2' = min(C2$probs), 'maxC2' = max(C2$probs),
  'minC3' = min(C3$probs), 'maxC3' = max(C3$probs),
  'minC4' = min(C4$probs), 'maxC4' = max(C4$probs))
##     minC1     maxC1     minC2     maxC2     minC3     maxC3     minC4     maxC4 
## 0.1029139 0.2424824 0.0104721 0.1029029 0.5076631 0.9999902 0.2425322 0.5057944

Entonces, notamos que el primer grupo es el \(C_2\), el segundo grupo es \(C_1\), el tercer grupo es \(C_4\) y último grupo, los que tiene altas probabilidades de mora, es \(C_3\).

Así, veamos si cumplen las pruebas estadísticas mencionadas anteriormente:

G1 <- C2; G2 <- C1; G3 <- C4; G4 <- C3

G1$GH <- as.character(G1$GH)
G2$GH <- as.character(G2$GH)
G1$GH <- as.character(G1$GH)
G2$GH <- as.character(G2$GH)

variables <- c(G1$probs, G2$probs)
categorias <- c(G1$GH, G2$GH)
bartlett.test(variables, categorias)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  variables and categorias
## Bartlett's K-squared = 2533.5, df = 1, p-value < 2.2e-16

Podemos ver que, como el \(p\_valor\) es menor a un \(\alpha = 0.05\), se rechaza la hipótesis nula y, en consecuencia, las varianzas entre el grupo 1 y grupo 2 son diferentes. De manera igual manera se corroboró para las otras combinaciones de grupos.

Con esto, realizamos la prueba de Games y Howell, para la comparación de medias entre los grupos:

banco$GH <- as.factor(banco$GH)
GH <- oneway(banco$GH, y = banco$probs, posthoc = 'games-howell')
GH
## ### Oneway Anova for y=probs and x=GH (groups: 1, 2, 3, 4)
## 
## Omega squared: 95% CI = [.89; .89], point estimate = .89
## Eta Squared: 95% CI = [.89; .89], point estimate = .89
## 
##                                     SS    Df     MS         F     p
## Between groups (error + effect) 528.74     3 176.25 108103.54 <.001
## Within groups (error only)       67.14 41184      0                
## 
## 
## ### Post hoc test: games-howell
## 
##      diff ci.lo ci.hi      t       df     p
## 2-1 -0.10 -0.10 -0.09 250.12 17543.34 <.001
## 3-1  0.53  0.52  0.54 101.50   763.60 <.001
## 4-1  0.18  0.18  0.19 157.10  4593.70 <.001
## 3-2  0.62  0.61  0.64 119.92   758.49 <.001
## 4-2  0.28  0.28  0.28 246.48  4015.39 <.001
## 4-3 -0.35 -0.36 -0.33  64.85   828.51 <.001

En consecuencia, vemos que los grupos homogéneos son heterogéneos entre sí; entonces, por medio del uso de la variable indicadora month, tenemos 8 tiempos para realizar el ciclo económico; de donde, por cada grupo y por cada mes, se calculará el procentaje de clientes malos \(y = no = 0\) en ese grupo:

banco$GH <- as.character(banco$GH)
library(plyr)
banco$GH <- revalue(banco$GH, c('1' = 'G2', '2' = 'G1', '3' = 'G4', 
                              '4' = 'G3'))
detach(package:plyr)

for(i in 1:length(banco$month)){
  if(banco$month[i] == 6 | banco$month[i] == 7 | banco$month[i] == 8){
    banco$month[i] <- 5
  }else if(banco$month[i] == 9){
    banco$month[i] <- 6
  }else if(banco$month[i] == 10){
    banco$month[i] <- 7
  }else if(banco$month[i] == 11 | banco$month[i] == 12 | banco$month[i] == 13 |
           banco$month[i] == 14 | banco$month[i] == 15 | banco$month[i] == 16 |
           banco$month[i] == 17 | banco$month[i] == 18 | banco$month[i] == 19 |
           banco$month[i] == 20 | banco$month[i] == 21 | banco$month[i] == 22 |
           banco$month[i] == 23 | banco$month[i] == 24 | banco$month[i] == 25 |
           banco$month[i] == 26){
    banco$month[i] <- 8
  }
}

no <- banco %>% filter(y == 0) %>% group_by(GH, month) %>% dplyr::summarise(tmalos = n())

yes <- banco %>% filter(y == 1) %>% group_by(GH, month) %>% dplyr::summarise(tbuenos = n())

mes <- full_join(no, yes, by = c('GH', 'month'))
mes$tmalos <- mes$tmalos/(mes$tmalos + mes$tbuenos)

Así también, asignamos la PD de cada grupo:

PD <- mes %>% group_by(GH) %>% dplyr::summarise(PD = mean(tmalos, na.rm = TRUE))
PD
## # A tibble: 4 x 2
##   GH       PD
##   <chr> <dbl>
## 1 G1    0.938
## 2 G2    0.856
## 3 G3    0.716
## 4 G4    0.495
ggplot(data=mes) + geom_line(aes(x = month, y = tmalos, color = GH), size=1) +
  labs(x = "Periodo", y = "PD",
       title = "Ciclo Económico") +
  theme_bw()+
  theme(legend.background = element_rect()) +
  scale_color_hue(name = 'GH')

mes <- full_join(no, yes, by = c('GH', 'month'))
mes$total <- mes$tmalos + mes$tbuenos
bdd <- mes %>% group_by(month) %>% mutate(porc = round(100*(total/sum(total,na.rm=T)),2))

ggplot(data=bdd,
       aes(x = factor(month, level= c("1","2","3","4","5","6","7","8")),
           y = porc, fill = GH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
  xlab("Meses") + ylab("Porcentaje") +
  labs(title = "Porcentaje de clientes por grupos en cada mes",
       fill= "Grupos") +
  scale_fill_brewer(palette="Spectral")

Conclusiones