Logistinen regressio ja luottoriskien estimointi R-ohjelmalla

Tässä tekstissä käsittelen logistista regressiota ja luottoriskien estimointia. Muodostan mallin, jolla pankin asiakkaiden maksukyvyttömyyden (default) todennäköisyyttä voidaan arvioida perustuen erilaisiin muuttujiin. Muuttujina ovat mm. asiakkaan luottohistoria, haetun luoton määrä, työllisyystilanne, säästöt, asumismuoto ja ikä. Data (german.data) ja lisätietoa datasta löytyy täältä. Menetelmät, joita käytän perustuvat D, Sharman oppaaseen, joka löytyy täältä. Datasta Default-muuttuja toimii selitettävänä muuttujana ja muut 20 muuttujaa selittävinä muuttujina.
Default-muuttuja on binäärimuuttuja, joka saa arvon 1 jos asiakas on maksukyvytön ja arvon 0 jos asiakas ei ole maksukyvytön, eli pystyi maksamaan luottonsa takaisin. Tarvitsemme siis luotonmyöntöön vaikuttavien tietojen lisäksi tiedon siitä, kuinka asiakkaille kävi, eli pystyivätkö he maksamaan luottonsa vai eivät. Tämän tiedon avulla voimme muodostaa mallin ja mahdollisesti tehdä havaintoja, kuten esim. että työssäkäyvät selviytyvät luotoistaan työttömiä useammin. Jos nyt kaksi muuten täysin identtistä asiakasta tulisi hakemaan pankista lainaa, niin mallimme mukaan työllinen selviytyisi luotostaan todennäköisemmin. Se, että toinen maksaa lainansa takaisin todennäköisemmin ei tarkoita välttämättä, että toinen ei saisi lainaa ollenkaan. Lainanantaja voi nimittäin korkoa kasvattamalla pitää tuotto-odotuksensa vakiona, eli vaikka maksukyvyttömyyttä esiintyisi enemmän tietyllä ryhmällä, niin suuremmat korot voivat kompensoida nämä tappiot. Alla on esimerkki yksinkertaisesta yhtälöstä, joka kuvaa riskillisen tuoton ja riskittömän tuoton yhteyttä. Riskitön tuotto "r(g)" voi kuvata esim. EKP:n talletuskorkoa, eli tuottoa, jonka pankki saa rahalleen ilman riskiä. Yhtälössä "(1-p)" on "ei maksukyvyttömyyden" todennäköisyys, "u" on raha, joka saadaan takaisin maksukyvyttömyyden sattuessa (oletetaan, että kaikki menetetään, jos asiakas on maksukyvytön, eli u=0), "z" riskipreemio ja "r(f)" riskillinen tuotto . Riskipreemio kuvaa tuottoa, jonka pankki haluaa lisäksi, kun riskillisen tuoton odotusarvo vastaa riskitöntä tuottoa. Yhtälöstä (1) voidaan ratkaista ns. spreadi, joka kuvaa riskillisen- ja riskittömän tuoton erotusta (S=r(f)-r(g)). Tämä spreadi on esitetty yhtälössä (2). (Huom. ratkaisussa on oletettu, että p ja r(f) ovat pieniä, jotta p*r(f)~0

1.  (1-p)(1+r(f)) + p*u = 1+r(g)+z
2.  S=p+z

Spreadiin vaikuttaa siis sekä maksukyvyttömyyden todennäköisyys että riskipreemio. Vielä yksinkertaistan, kun oletetaan, että riskipreemio on nolla, voimme kirjoittaa spreadin S=p eli r(f)-r(g)=p

Jos siis pankki ei pyydä riskipreemiota ja haluaa odotusarvoisesti saada saman tuoton sekä työlliseltä että työttömältä, niin sen pitää tietää p mahdollisimman tarkasti. Jos työllisen todennäköisyys maksukyvyttömyyteen olisi 1%:n ja työttömän vastaava luku 3%:a, niin pankki pyytäisi esimerkin tilanteessa työlliseltä korkoa ¨~1%+riskitön korko ja vastaavasti työttömältä  ~3%+riskitön korko. Tällä tavoin pankin odotusarvo kummallekin saatavalla olisi sama.


Nyt takaisin malliin ja sen luomiseen. Ensin lataan R:ään datan ja nimeän sarakkeet, koska niillä ei ole valmiina nimiä.
R-koodi:

rm(lists=ls())    #poistaa aiemmat tallennetut muuttujat. Hyvä käyttää, jos attach-komentoa käytetään.
data <- read.delim(file.choose(), header=FALSE,sep="")

names(data) <- c("check_acc", "duration", "history", "purpose", "amount", "sav_acc", "employment", "install_rate", "ptatus", "other_debtor", "resid", "property", "age", "other_install", "housing", "other_credits", "job", "num_depend", "teleph", "foreign", "default")
attach(data)
str(data)                              #str-komennolla nähdään ovatko muuttujat kategoriallisia vai jatkuvia

Seuraavaksi voimme katsoa dataa kuvien avulla. Jatkuvia muuttujia, kuten age voidaan katsoa esim. plot-ja hist-komennoilla. Muuttujia, jotka ovat faktoreina (esim.history) voidaan katsoa lattice-pakkauksen histogram-komennolla. Katsotaan seuraavaksi age-muuttujaa:
R-koodi:

library(lattice)

histogram(age,type="percent",breaks=seq(10,80,by=10))   #breaks-komennolla jaetaan eri ikäryhmiin

Näemme, että aineistossa suurin osa on 20-40-vuotiaita. Voimme tarkastella eri ikäluokkien osuuksia default-luokissa vcd-pakkauksen mosaic-komennolla. Jaetaan age-muuttuja ensin kolmeen ikäluokkaan ja muutetaan default-muuttuja faktoriksi, joka saa arvot 0 ja 1. Sen arvo 1 tarkoittaa maksukyvyttömyyttä. Sitten voimmekin tehdä Mosaic-kuvan.
R-koodi:

agec <- as.factor(ifelse(age <= 30, "0-30",ifelse(age <=40,"31-40","41-80")))
data$default <- as.factor(ifelse(data$default==1,"0","1")) #päivitetään myös data.frameen
attach(data)
library(vcd)
mosaic(~default+age,highlighting="age",main="age & default")

Kuvasta nähdään, että 0-30-vuotiaiden osuus on maksukyvyttömistä suurin. Seuraavaksi tutkin, mitkä ovat maksukyvyttömien osuudet eri ikäluokissa. Tämän voi tehdä table-komennolla seuraavasti: 
R-koodi:

defaultagetable <- table(default,age) # Defaultin luokat menevät table:n riveiksi ja age:n luokat sarakkeiksi
prop.table(defaultagetable,margin=2) # Lasketaan maksukyvyttömien ja maksukykyisten osuudet
a <- barplot(prop.table(defaultagetable,margin=2),beside=T,legend.text=c("Maksukykyisten osuus",
"Maksukyvyttömien osuus"),yaxt="n",ylab="%",xlab="age") # Tehdään pylväskuvio
axis(2,at=seq(0,1,by=0.1))
text(a, 0, round(fdf,2),cex=1.1,pos=3)

 Testaan nyt Khiin neliö -testillä (Chi-square test), ovatko nämä erot ikäluokkien välillä tilastollisesti merkitseviä. "0-30" -luokka eroaa muistaa luokista merkitsevästi, kun taas "31-40" ja "41-80" -luokat eivät eroa keskenään merkitsevästi.

Seuraavaksi lisään "age & default" -kuvaan kolmannen muuttujan, joka on other_debtor, eli muuttuja, joka on kuvattu aineistokuvauksessa seuraavasti:

Other debtors / guarantors 
A101 : none 
A102 : co-applicant 
A103 : guarantor 


Muuttuja siis kertoo, onko hakijalla lainalle takaaja ja hakeeko hakija lainaa yksin. Katsotaan seuraavaksi kuvaa.
R-koodi:

mosaic(default~age+other_debtor,main="age, default & other_debtor")

Kuvasta nähdään, että jokaisessa ikäryhmässä "co-applicant"-hakijat olivat suhteellisesti useimmin maksukyvyttömiä. Hakijat, joilla oli takaaja, pystyivät joka ikäluokassa suhteellisesti useimmin maksamaan lainansa takaisin. 

Seuraavaksi muutan osan muuttujista faktori-muotoon, koska ne eivät ole numeerisia muuttujia eli ne saavat arvoja esim. välillä 1-4 ja kuvaavat luokittelua. 
R-koodi:

data$install_rate <- as.factor(data$install_rate)
data$resid <- as.factor(data$resid)
data$other_credits <- as.factor(data$other_credits)

data$num_depend <- as.factor(data$num_depend)

Nyt  jaan datan "trainingdataan" ja "testdataan", eli jaan datan kahteen osaan, joista toisesta (trainingdata) muodostan mallin, jolla ennustan toisen arvoja (testdata). Tosiasiassa mallin toimivuus tulisi testata myös cross-validationilla ja/tai jakamalla data kolmee osaan, joista kolmatta validationsettiä ei käytetä mallia luodessa, vaan sitä kokeillaan vain kerran, kun malli on valmis. Kun mallia verrataan testdataan, niin muuttujia lisäillessä ja mallia "tuunatessa" voi syntyä ns. overfitting, jolloin malli toimii testdataan paremmin kuin mallin ulkopuoliseen validationsettiin. Tavoitteena on siis muodostaa malli, johon syöttämällä testdatan muuttujien arvot voidaan saada mahdollisimman hyviä ennusteita maksukykyisyydestä. Testdatastahan voimme tarkistaa, kuinka hyvin malli on onnistunut, koska testdatassa on default-muuttuja, joka kertoo maksoiko asiakas lainansa takaisin vai ei (tosiasiassa 1-arvo voi tarkoittaa myös, että asiakas on maksanut lainansa myöhässä tai lainanmaksussa on ollut muita vaikeuksia). Otan trainingdataan sattumanvaraisesti 60%:a arvoista, eli testdataan jää 40%:a arvoista. Huom. Datassa on tällä hetkellä mukana myös jatkuvia muuttujia, kuten age ja amount, joten teen aluksi mallin normaalilla datalla. Kokeilen tämän jälkeen mallia, jossa on pelkästään kategoriallisia muuttujia, eli myös muuttujat, kuten age ja amount on muutettu kategorioihin. Sharma suosittelee käytettäväksi mallissa vain kategorisia muuttujia, mutta toisaalta informaatiokadon vuoksi myös jatkuvien muuttujien malleilla on kannattajansa. Jatkuvan muuttujan muuttaminen kategorialliseksi voi olla hyödyllistä, jos muuttujan yhteys selitettävään muuttujaan ei ole lineaarinen. Epälineaarisia yhteyksiä voi toisaalta mallintaa myös esim. neliöjuuri-muunnoksella, jolloin mallissa käytetään muuttujaa sekä normaalissa että neliöjuuri-muodossa. Epälineaarisuuksia tarkastellessa on aina tarkasteltava muuttujaa myös eri kategorioissa, koska epälineaarisuus saattaa johtua  muuttujan erilaisesta käyttäytymisestä jossain kategoriassa (yhteyden ollessa lineaarinen muissa kategorioissa). Jos muuttuja käyttäytyy eri tavoin eri kategorioissa, on interarction-termien käyttö suositeltavaa. Suosittelen, että muuttujia tarkastellaan aina ennen mallin tekemistä huolella ja mahdollinen kategorisointi on pystyttävä perustelemaan teorialla. Lähtökohtaisesti pyrkisin pitämään jatkuvan muuttujan jatkuvana. Joskus selitettävän ja selittävän -muuttujan välillä löytyy yhteys vasta esim. logaritmi-muunnoksen jälkeen, joten eri muunnokset on hyvä käydä läpi ennen kuin lähtee kategorisoimaan muuttujaa.
R-koodi:

e <- sample(nrow(data), 0.6*nrow(data))  #ottaa sattumanvaraisesti 60%:a rivinumeroista
trainingdata <- data[e,]                              # valitaan saadut rivinumerot trainingdataan ja testdataan

testdata <- data[-e,]

Nyt trainingdatasta voidaan muodostaa logistinen regressiomalli. 
R-koodi: 

model1 <- glm(default~check_acc+duration+history+purpose+
amount+
sav_acc+
employment+
install_rate+
ptatus+
other_debtor+
resid+
property+
age+
other_install+
housing+
other_credits+
job+
num_depend+
teleph+

foreign,family=binomial(link=logit),data=trainingdata)
summary(model1)

Mallissa 1 on useampia tilastollisesti merkitseviä muuttujia. Alla tuloksia. (Standard error suluissa)

Results of the model1

	Dependent variable:
	default
check_accA12	-0.550*
	(0.282)
check_accA13	-0.933**
	(0.461)
check_accA14	-1.944***
	(0.321)
duration	0.036***
	(0.012)
historyA31	0.849
	(0.727)
historyA32	-0.059
	(0.565)
historyA33	-0.240
	(0.613)
historyA34	-0.527
	(0.557)
purposeA41	-1.859***
	(0.496)
purposeA410	-2.696**
	(1.134)
purposeA42	-0.553
	(0.346)
purposeA43	-0.809**
	(0.319)
purposeA44	-2.336
	(1.480)
purposeA45	-0.166
	(0.650)
purposeA46	0.092
	(0.563)
purposeA48	-0.696
	(1.465)
purposeA49	-0.297
	(0.423)
amount	0.0001
	(0.0001)
sav_accA62	-0.541
	(0.377)
sav_accA63	-0.532
	(0.502)
sav_accA64	-1.364*
	(0.734)
sav_accA65	-1.044***
	(0.365)
employmentA72	0.759
	(0.637)
employmentA73	0.454
	(0.617)
employmentA74	0.070
	(0.654)
employmentA75	0.517
	(0.614)
install_rate2	0.150
	(0.389)
install_rate3	0.332
	(0.439)
install_rate4	0.731*
	(0.387)
ptatusA92	-0.084
	(0.490)
ptatusA93	-0.630
	(0.477)
ptatusA94	0.296
	(0.579)
other_debtorA102	0.331
	(0.637)
other_debtorA103	-0.844*
	(0.505)
resid2	1.186***
	(0.400)
resid3	0.872**
	(0.429)
resid4	0.452
	(0.403)
propertyA122	0.446
	(0.325)
propertyA123	0.139
	(0.312)
propertyA124	0.950*
	(0.554)
age	-0.019
	(0.012)
other_installA142	-0.682
	(0.651)
other_installA143	-0.520
	(0.318)
housingA152	-0.391
	(0.314)
housingA153	-0.221
	(0.636)
other_credits2	0.455
	(0.312)
other_credits3	0.073
	(0.848)
other_credits4	0.851
	(1.384)
jobA172	-0.169
	(0.807)
jobA173	-0.198
	(0.779)
jobA174	0.257
	(0.807)
num_depend2	0.261
	(0.328)
telephA192	-0.372
	(0.265)
foreignA202	-2.272***
	(0.854)
Constant	-0.051
	(1.309)
Observations	600
Log Likelihood	-270.011
Akaike Inf. Crit.	650.022
Note:	*p<0.1; **p<0.05; ***p<0.01

>

Nyt voimme arvioida mallin ennustustarkkuutta laskemalla ROC-käyrän.
R-koodi:

testdata$model1p <- predict(model1, testdata, type = "response") # predict-komento laskee mallilla ennusteet testdatan arvoilla ja muodostaa testdataan uuden "model1p"-sarakkeen
library(ROCR)
model1points <- prediction(testdata$model1p, testdata$default) #vertaa ennusteita oikeisiin arvoihin
performancemodel1 <- performance(model1points,"tpr","fpr")
plot(performancemodel1, col = "blue")     #kuvaa ROC-käyrän
abline(0,1,lty=5)                                         #tekee 45 asteen suoran viivan

ROC-käyrän muotoja ja mallin onnistumistarkkuutta havainnollistaa alla oleva kuva.

Lähde: http://gim.unmc.edu/dxtests/roccomp.jpg

ROC-käyrässä True positive rate (TPR) tarkoittaa lukua, joka saadaan, kun jaetaan onnistuneet luotonhylkäykset, kun luotonsaaja on maksukyvytön, kaikilla maksukyvyttömien määrällä. False positive rate (FPR) tarkoittaa lukua, joka saadaan, kun jaetaan luoton hylkäykset, kun asiakas olisi ollut maksukykyinen, kaikkien maksukykyisten määrällä. Käyrän arvot saadaan, kun käytetään erilaisia ns. cut-off -arvoja ja lasketaan TPR ja FPR näillä arvoilla. Esim. kun cut-off-arvo on 0.1, niin luotot hylätään, jos asiakkaan maksukyvyttömyyden todennäköisyys on yli tämän. Hylätyistä luotoista lasketaan sitten oikein ja väärin luokitellut, jolloin saadaan ROC-käyrän piste tälle cut-off-arvolle.  ROC-käyrästä näemme, että cut-off arvo alkaa luvusta 1 ja oikealle mentäessä se laskee aina nollaan saakka, jolloin käyrät kohtaavat. Lasken vielä AUC-arvon, joka kertoo, kuinka paljon aluetta ROC-käyrän alle jää. Tämä auttaa ROC-käyrien vertailussa, kun ero ei visuaalisesti ole aivan selvä. ROC-käyrille saa laskettua myös luottamusvälit. 
R-koodi:

model1_auc <- performance(model1points, "auc")

Mallin AUC-arvo on 0.7898

Lasken myös Kolmogorov-Smirnov (KS) -arvon
R-koodi:


ks1 <- max(attr(performancemodel1, "y.values")[[1]] - (attr(performancemodel1 "x.values")[[1]]))



KS-arvo kertoo FPR:n ja TPR:n kumulatiivisen erotuksen maksimikohdan. KS-arvo 1 tarkoittaa, että malli on onnistunut jakamaan maksukyvyttömät ja maksukykyiset selkeästi kahteen eri luokkaan. KS-arvo 0 tarkoittaisi, ettei malli ole onnistunut luokkiin jakamisessa. KS-arvo 1 voitaisiin saada esim. jos kaikki maksukykyiset saisivat arvon <0.5 ja kaikki maksukyvyttömät arvon >= 0.5. Tällainen KS-arvo tarkoittaisi, että voisimme laittaa cut-off arvoksi 0.5, jolloin hylkäisimme kaikki maksukyvyttömien lainat ja myöntäisimme kaikki maksukykyisten lainat. 

KS-arvosta löytyy lisää tietoa mm. täältä. 

Kokeilin myös mallia, jossa jatkuvat muuttujat (age, duration ja amount) oli muutettu kategorisiksi, mutta mallin tulokset eivät eronneet merkittävästi. Mallia sai paranneltua, kun lisäsi siihen interaction-muuttujia. Potentiaalisia interaction-muuttujia löytää Machine learningiin (Koneoppiminen) kuuluvan Random Forestin Variable importance -listasta. Parhaan mallin löysin käyttämällä Koneoppimisen menetelmiä, joilla sain AUC-arvon, joka oli suurempi kuin D. Sharman paperissaan saama suurin arvo. Koneoppiminen on erittäin hyödyllinen menetelmä mallintamisessa, koska sen keinoin saadaan useista malleista luotua ennusteita, jotka peittoavat yksittäisten mallien tulokset. Näistä menetelmistä kirjoittelen jatkossa tarkemmin.