libname fic 'C:\Documents and Settings\hraissi\Bureau\Sens\exemple';


/******LA PARTIE PARAMETRIQUE - REGRESSION LOGISTIQUE ****/
data fic.bdf;

infile 'C:\Documents and Settings\hraissi\Bureau\Sens\Ks.txt' firstobs=2;
input y x1-x20;
x2 = x2/100; 
x5 = x5/10000; /*on ramene les variables a la meme echelle */
run;

data fic.bdf2;
set fic.bdf;

d1a=0; d1b=0; d1c=0;  /* x1 balance of current account */
if x1 = 2 then d1a = 1;
if x1 = 3 then d1b = 1;
if x1 = 4 then d1c = 1;
drop x1;

/* variable x2 - duration in month */

d3a=0; d3b=0; d3c=0; d3d=0;  /* x3 payment of previous credits */
if x3 = 1   then d3a = 1;
if x3 = 2 then d3b = 1;
if x3 = 3  then d3c = 1;
if x3 = 4 then d3d = 1;
drop x3;

d4a=0; d4b=0; d4c=0; d4d=0;  /* x4 purpose of credit */
if x4 = 1 or x4 = 2  then d4a = 1;
if x4 = 3 then d4b = 1;
if x4 = 4 or x4 = 5 or x4 = 6 then d4c = 1;
if x4 = 9 then d4d = 1;
drop x4;

/* variable x5 - amount of credit */

if  x6 <= 2 then x6 = 0;
if x6 > 2 then x6 = 1; /* x6 value of savings  */


d7a=0; d7b=0; d7c=0; d7d=0;  /* x7 employement */
if x7 = 2 then d7a = 1;
if x7 = 3 then d7b = 1;
if x7 = 4 then d7c = 1;
if x7 = 5 then d7d = 1;
drop x7;

d8a=0; d8b=0; d8c=0;  /* x8 installement in % of available income */
if x8 = 2 then d8a = 1;
if x8 = 3 then d8b = 1;
if x8 = 4 then d8c = 1;
drop x8;

d9a=0; d9b=0; d9c=0;  /* x9 marital status */
if x9 = 2 then d9a = 1;
if x9 = 3 then d9b = 1;
if x9 = 4 then d9c = 1;
drop x9;


d10a=0; d10b=0;   /* x10 further debtors/ guarantors */
if x10 = 2 then d10a = 1;
if x10 = 3 then d10b = 1;
drop x10;


d12a=0; d12b=0; d12c=0;  /* x12 most valuable available assets */
if x12 = 2 then d12a = 1;
if x12 = 3 then d12b = 1;
if x12 = 4 then d12c = 1;
drop x12;

d14a=0; d14b=0;   /* x14 further running credits */
if x14 = 2 then d14a = 1;
if x14 = 3 then d14b = 1;
drop x14;

d15a=0; d15b=0;   /* x15 is the type of appartment */
if x15 = 2 then d15a = 1;
if x15 = 3 then d15b = 1;
drop x15;

if x16 >1 then x16 = 0 ; /* x16 is the number of previous credits at this banks*/
							/*including the running one. If x16 =0 then there are */
							/* at least two credits */
d17a=0; d17b=0; d17c=0;
if x17 = 2 then d17a = 1;
if x17 = 3 then d17b = 1;
if x17 = 4 then d17c = 1;
drop x17;
run;


/* POUR ILLUSTRATION, UN NOMBRE REDUIT DES VARIABLES EST SELECTIONE, 
LES MEMES QUE DANS LE LIVRE DE REFERENCE */

data fic.bdf3;
set fic.bdf2 ;
keep y d1a d1b d1c x2 x5 d7a d7b d7c d7d d8a d8b d8c;
run;


proc freq data = fic.bdf3; tables y; run;


data fic.est fic.test;
set fic.bdf3;
aux = uniform(1);
if aux > 0.15 then output fic.est;
else output fic.test;
run;

proc freq data = fic.est;
tables y; run;

proc logistic data = fic.est descending;
model y = d1a d1b d1c x2 x5 d7a d7b d7c d7d d8a d8b d8c ;
output out=fic.logit_estun pred=pun;
run;

/*proc logistic data = fic.est descending;
model y = x5 d7a d7b d7c d7d d8a d8b d8c ;
output out=fic.logit_estdeux pred=pdeux;
run;*/
/*proc print data=fic.test;
run;*/
PROC IML;
use fic.test;
read all var{d1a d1b d1c x2 x5 d7a d7b d7c d7d d8a d8b d8c} into logit;

coef_logit = { 
0.3970 ,     
0.9136 ,     
1.9611 ,     
-3.2764,      
-0.5294,      
0.0922 ,     
0.3298 ,     
0.8678 ,     
0.6812 ,     
-0.3013,      
-0.6471,      
-0.8799};     

n=nrow(logit);
xbeta=J(n,1,0);
inter=J(n,1,1.1626);/*1.1626 est l'intercept estimé*/
xbeta=logit*coef_logit+inter;

prob=J(n,1,0);
do i=1 to n;
prob[i]=exp(xbeta[i])/(1+ exp(xbeta[i]));
end;

CREATE result
            var{p};
CLOSE result;

/* on a créé une table pour les résultats */

pre=J(n,1,0);
pre=prob;

print pre;
EDIT result ;  Append from pre; CLOSE result ;
run;
quit;

proc print data=result;
run;

data fic.score; 
merge result fic.test; 
keep p y ;
run;
proc print data=fic.score;
run;

proc sort data=fic.score out=fic.sco;
by p;
run;


/*proc sort data = fic.logit_estun out = fic.prev_un;
by pun; run;
quit;

proc sort data = fic.logit_estdeux out = fic.prev_deux;
by pdeux; run;
quit;

data fic.score; 
merge fic.prev_un fic.prev_deux;
by score_logit; 
keep pun pdeux y ;
run;

proc print data= fic.score;
run;*/

/****************************************************************************************
** Macro de traçage des courbes de performance, sélection et discrimination d'un score **
*****************************************************************************************
** Etienne MAROT, Groupe de Recherche Opérationnelle, Crédit Lyonnais (2001)		   **
*****************************************************************************************
** 																					   **
** PARAMETRES :																		   **
*****************************************************************************************
** 																					   **
** DATA			La table à utiliser													   **
** SCORE		La variable de score (ou les deux variables de score à comparer)	   **
** CRITERE		Le critère binaire modélisé	(1 pour les bons, 0 pour les mauvais)	   **
** TAUX			Le taux d'échantillonnage pour le traçage des courbes (défaut: 1)	   **
** DEVICE		Pilote pour l'impression des courbes (défaut: win)					   **
** COURBES		Série de 4 booleens pour resp. la courbe de perf, sel et disc, et	   **
**				pour le graphe comparant la densité sur la pop. des défaillants	 et	   **
**				celle des non défaillants approchées par la méthode du noyau		   **
**				(par défaut: 1 1 1 1)												   **
**																					   **
*****************************************************************************************
**																					   **
** EXEMPLE D'APPEL :																   **
** %courbes(data=a,score=score,critere=y,courbes=1 1 1 1,taux=1);					   **
**																					   **
****************************************************************************************/

%macro courbes(data=,score=,critere=,taux=1,courbes=1 1 1 1,device=win) /des="Courbes caractéristiques d'un score";

options nomprint nonotes nomlogic nosymbolgen; /* Moi j'aime pas qu'on regarde mon code */

%let crbperf=%scan(&courbes,1); /* Récupération des choix d'affichage des courbes */
%let crbsel=%scan(&courbes,2);
%let crbdisc=%scan(&courbes,3);
%let crbdens=%scan(&courbes,4);

/* Controle de la validité des paramètres */
%if (&crbsel ne 1 and &crbsel ne 0) or (&crbperf ne 1 and &crbperf ne 0) or 
	(&crbdisc ne 1 and &crbdisc ne 0) or (&crbdens ne 1 and &crbdens ne 0)%then %do;
	%put ;
	%put ############ Mauvais paramétrage de l%str(%')affichage des courbes.;
	%goto fin;
%end;
%if &taux<=0 or &taux>1 %then %do;
	%put ############ Mauvais taux d%str(%')échantillonnage : il doit être dans %nrstr(]0;1]).;
	%goto fin;
%end;

/* Si on a plusieurs scores à comparer : on prend les deux premiers */
%let score2=%scan(&score,2);
%let score=%scan(&score,1);
%if &score2 ne %then %let s=s; %else %let s=; /* pour les titres des graphiques */

/* A partir de là tout va bien, donc on peut y aller */
goptions device=&device keymap=winansi devmap=winansi; /* Pour afficher à l'écran et avoir les accents dans les titres */

proc datasets nolist; /* On efface le catalogue de sortie au cas où il existe : c'est pour l'unicité des noms de graphiques */
delete __crb (memtype=catalog);
quit;

proc sort data=&data(keep=&score &score2 &critere where=(ranuni(0)<=&taux)) out=__tmp; /* Tri décroissant selon le score, avec échantillonnage */
by descending &score;
run;

%if &score2 ne %then %do; /* s'il y a un second score, on fait un tri sur les mêmes individus */
	proc sort data=__tmp out=__tmp2(rename=(&critere=&critere.2));
	by descending &score2;
	run;

	data __tmp; /* on juxtapose les deux tables */
	merge __tmp __tmp2;
	run;
%end;


%if &SYSERR>4 %then %do; /* s'il y a une erreur lors de cette étape, on ne peut rien faire après => on arrête tout */
	%put ############ Erreur d%str(%')exécution. Mauvais paramétrage.;
	%goto fin;
%end;

data __tmp; /* Lecture de la table par note, et comptage des défaillants et des individus */
set __tmp;
nbIND=_N_;
retain nbDEF1 nbDEF2 0;
nbDEF1=nbDEF1+(1-&critere);
%if &score2 ne %then %do;
	nbDEF2=nbDEF2+(1-&critere.2);
%end;
run;

/* Récupération du nombre d'individus et du nombre de défaillants */

proc sql noprint;
select count(&critere), sum(1-&critere) into: nbINDT separated by ' ', :nbDEFT separated by ' ' from __tmp;
quit;

data __tmp; /* Calcul de la performance, de la sélection et de la discrimination */
set __tmp;
nbIND=nbIND/&nbINDT; /* Valeur de l'abscisse */
%if &crbsel=1 %then %do;
	nbSEL1=nbDEF1/&nbDEFT; /* Courbe de sélection */
	refSEL=max(1+(nbIND-1)/(&nbDEFT/&nbINDT),0); /* Courbe de sélection du score optimal */
	if nbIND ne 1 then do;
		Gini1=nbIND-nbSEL1; /* Calcul des différences d'aires pour l'indice de Gini */
		Gini0=nbIND-refSEL;
	end; 
	%if &score2 ne %then %do;
		nbSEL2=nbDEF2/&nbDEFT;
		if nbIND ne 1 then do;
			Gini2=nbIND-nbSEL2;
		end;
	%end;
%end;
%if &crbperf=1 %then %do;
	nbPERF1=nbDEF1/&nbDEFT/nbIND; /* Courbe de performance */
	%if &score2 ne %then %do;
		nbPERF2=nbDEF2/&nbDEFT/nbIND;
	%end;
%end;
%if &crbdisc=1 %then %do;
	nbPASBON1=nbDEF1/&nbDEFT; /* Abscisse de la courbe de discrimination */
	nbBON1=(nbIND*&nbINDT-nbDEF1)/(&nbINDT-&nbDEFT); /* Ordonnée de la courbe de discrimination */
	%if &score2 ne %then %do;
		nbPASBON2=nbDEF2/&nbDEFT;
		nbBON2=(nbIND*&nbINDT-nbDEF2)/(&nbINDT-&nbDEFT);
	%end;
%end;
run;

%if &crbsel=1 %then %do; /* Calcul de l'indice de Gini par sommation de la variable créée plus haut */

	proc sql noprint; /* on somme les aires individuelles */
	%if &score2 ne %then %do;
		select sum(Gini1), sum(Gini2), sum(Gini0) into: Gini1 separated by ' ', :Gini2 separated by ' ', : Gini0 separated by ' ' from _LAST_;
	%end;
	%else %do;
		select sum(Gini1), sum(Gini0) into: Gini1 separated by ' ', : Gini0 separated by ' ' from _LAST_;
	%end;
	quit;
	%let Gini1=%sysevalf(&Gini1/&Gini0); /* et le ratio vaut l'indice de Gini */
	%if &score2 ne %then %let Gini2=%sysevalf(&Gini2/&Gini0); 
%end;

%if %eval(&crbsel+&crbperf+&crbdisc+&crbdens)=4 %then %do; /* Si les 4 courbes sont demandées, on les affiche ensemble */
	goptions nodisplay; /* Et du coup on n'affiche pas les autres */
%end;

proc gplot data=__tmp gout=__crb; /* Et on se lance pour faire des superbes graphiques */
%if &crbsel=1 %then %do; /* Représentation de la courbe de sélection */
	title "Selection curves (Lorenz curves)";
	axis1 label=(color=black justify=center 
		'Proportion of selected applicants (from higher to lowest score)'); /*Part des individus sélectionnés');*/
	axis2 label=(angle=90 color=black justify=center 'Proportion of defaults (Y=0)') order=(0 to 1 by 0.1);
	%if &score2= %then %do;
		plot (nbSEL1 nbIND refSEL Gini1)*nbIND /noframe grid haxis=axis1 vaxis=axis2 name='crbsel' overlay legend=legend1;
		symbol1 i=join v=none c=red width=2;
		symbol2 i=join v=none c=black width=2 line=2;
		symbol3 i=join v=none c=green width=2 line=4;
		symbol4 i=none v=none c=white width=1;
		legend across=1 cborder=black position=(middle inside left) offset=(10,9) mode=share 
			value=(tick=1 'Single-index Score' 
					tick=2 'Random Score function' 
					tick=3 'Optimal Score function' 
					tick=4 "Gini : &Gini1") shape=symbol(7,3) label=none;
	%end;
	%else %do;
		plot (nbSEL1 nbSEL2 nbIND refSEL Gini1 Gini2)*nbIND /noframe grid haxis=axis1 vaxis=axis2 name='crbsel' overlay legend=legend1;
		symbol1 i=join v=none c=red width=2 line=1;
		symbol2 i=join v=none c=orange width=2 line=1;
		symbol3 i=join v=none c=black width=2 line=2;
		symbol4 i=join v=none c=green width=2 line=4;
		symbol5 i=none v=none c=white width=1;
		symbol6 i=none v=none c=white width=1;
		legend across=1 cborder=black position=(middle inside left) offset=(10,7) mode=share 
			value=(tick=1 'Single-index Score' tick=2 'Logit Score' tick=3 'Random Score' 
tick=4 'Optimal Score' tick=5 "Gini 1 : &Gini1" tick=6 "Gini 2 : &Gini2") shape=symbol(7,3) label=none;
	%end;
	run;
%end;


%if &crbperf=1 %then %do; /* Représentation de la courbe de performance */
	title "Performance Curve";
	axis1 label=(color=black justify=center 
'Proportion of selected applicants (from higher to lowest score)');/*
Part des individus sélectionnés');*/
	axis3 label=(angle=90 color=black justify=center 
'Performance (#{Y=0} in selection / #{Y=0})');
	%if &score2= %then %do;
		plot nbPERF1*nbIND /noframe grid haxis=axis1 vaxis=axis3 vref=1 lvref=2 name='crbperf';
		symbol1 i=join v=none c=red width=2;
	%end;
	%else %do;
		plot (nbPERF1 nbPERF2)*nbIND /noframe grid haxis=axis1 vaxis=axis3 vref=1 lvref=2 name='crbperf' overlay;
		symbol1 i=join v=none c=red width=2 line=1;
		symbol2 i=join v=none c=orange width=2 line=1;
	%end;
	run;
%end;


%if &crbdisc=1 %then %do; /* Représentation de la courbe de discrimination */
	title "ROC curve";
	axis1 label=(color=black justify=center 'Proportion of defaults 
				(from higher to lowest score)'); /*Part de défaillants');*/
	axis4 label=(angle=90 color=black justify=center 'Proportion of non-defaults (Y=1)'); /*Part de non-défaillants');*/
	%if &score2= %then %do;
		plot (nbBON1 nbPASBON1)*nbPASBON1 /noframe grid haxis=axis1 vaxis=axis4 vref=1 lvref=2 href=1 lhref=2 name='crbdisc' overlay;
		symbol1 i=join v=none c=red width=2;
		symbol2 i=join v=none c=black width=2 line=3;
		run;
	%end;
	%else %do;
		quit;

		data __tmp2;
		merge __tmp(keep=nbPASBON1 nbBON1 rename=(nbPASBON1=nbPASBON)) __tmp(keep=nbPASBON2 nbBON2 rename=(nbPASBON2=nbPASBON));
		by nbPASBON;
		run;

		proc gplot data=__tmp2 gout=__crb;
		plot (nbBON1 nbBON2 nbPASBON)*nbPASBON /noframe grid haxis=axis1 vaxis=axis4 vref=1 lvref=2 href=1 lhref=2 name='crbdisc' overlay;
		symbol1 i=join v=none c=red width=2 line=1;
		symbol2 i=join v=none c=orange width=2 line=1;
		symbol3 i=join v=none c=black width=2 line=3;
		run;
	%end;
%end;
quit;

%if &crbdens=1 %then %do; /* Calcul et représentation des deux densités lissées (noyau) sur un meme graphique */

	proc sql noprint; /* on récupère le min et le max du score */
	select min(&score), max(&score) into:min1 separated by ' ', :max1 separated by ' ' from __tmp;
	quit;
	
	%if &score2 ne %then %do;
		proc sql noprint; /* on récupère le min et le max du score */
		select min(&score2), max(&score2) into:min2 separated by ' ', :max2 separated by ' ' from __tmp;
		quit;
	%end;

	%let ampli1=%sysevalf((&max1-(&min1))/200); /* Amplitude du score, pour régler l'abscisse du graphique */
	%if &score2 ne %then %let ampli2=%sysevalf((&max2-(&min2))/200); 

	ods exclude all; /* on ne veut pas afficher les résultats des proc KDE */

	proc kde data=__tmp(where=(&critere=0)) out=__f0 bwm=2.5; /* Densité sur les défaillants : méthode du noyau */
	var &score;
	run;

	proc kde data=__tmp(where=(&critere=1)) out=__f1 bwm=2.5; /* Densité sur les non-défaillants : idem */
	var &score;
	run;

	data __tmp1; /* On regroupe pour avoir une abscisse commune */
	merge __f0(rename=(density=dens0)) __f1(rename=(density=dens1));
	by &score;
	run;
	
	%if &score2 ne %then %do;
		proc kde data=__tmp(where=(&critere=0)) out=__f0 bwm=0.95; /* Densité sur les défaillants : méthode du noyau */
		var &score2;
		run;

		proc kde data=__tmp(where=(&critere=1)) out=__f1 bwm=2.5; /* Densité sur les non-défaillants : idem */
		var &score2;
		run;

		data __tmp2; /* On regroupe pour avoir une abscisse commune */
		merge __f0(rename=(density=dens0)) __f1(rename=(density=dens1));
		by &score2;
		run;
	%end;

	ods exclude none; /* on rétablit les sorties comme il faut */

	proc gplot data=__tmp1 gout=__crb;
	title "Overlapping of Default and non-Default applicants densities";
	%if &score2 ne %then %do;
		title2 "Single-index Score";
	%end;
	axis1 label=(color=black justify=center 'Value of the Score');
	axis5 label=(angle=90 color=black justify=center 'Density');
	plot (dens1 dens0)*&score /noframe grid haxis=axis1 vaxis=axis5 name='crbdens1' overlay legend=legend1;
	symbol1 i=join v=none c=green width=2;
	symbol2 i=join v=none c=red width=2;
	legend across=1 cborder=black position=(top inside right) offset=(-10,-5) mode=share 
		value=(tick=1 'Non-Default' tick=2 'Default') shape=symbol(7,3) label=none;
	run;
	quit;
	
	%if &score2 ne %then %do;
		proc gplot data=__tmp2 gout=__crb;
		title "Overlapping of Default and non-Default applicants densities";
		title2 "Logit Score";
		axis1 label=(color=black justify=center 'Value of the Score');
		axis5 label=(angle=90 color=black justify=center 'Density');
		plot (dens1 dens0)*&score2 /noframe grid haxis=axis1 vaxis=axis5 name='crbdens2' overlay legend=legend1;
		symbol1 i=join v=none c=blue width=2;
		symbol2 i=join v=none c=orange width=2;
		legend across=1 cborder=black position=(top inside right) offset=(-10,-5) mode=share 
			value=(tick=1 'Non-Default' tick=2 'Default') shape=symbol(7,3) label=none;
		run;
		quit;
	%end;
%end;

goptions cback=white display; /* On remet les options graphiques normalement */

%if %eval(&crbperf+&crbsel+&crbdisc+&crbdens)=4 %then %do; /* Affichage des graphes s'il y en a 4 */
	proc greplay tc=__crb nofs igout=__crb gout=__crb;  /* Création du patron pour l'affichage des 4 graphes */
	tdef crb  des="The curves for a score"
        1/llx=0   lly=0 
          ulx=0   uly=50 
          urx=50  ury=50 
          lrx=50  lry=0 
          color=white 
 
        2/llx=0   lly=50 
          ulx=0   uly=100 
          urx=50  ury=100 
          lrx=50  lry=50 
          color=white 
 
        3/llx=50  lly=50 
          ulx=50  uly=100 
          urx=100 ury=100 
          lrx=100 lry=50 
          color=white 
	
	%if &score2= %then %do; 
        4/llx=50  lly=0 
          ulx=50  uly=50 
          urx=100 ury=50 
          lrx=100 lry=0 
          color=white
	%end;
	%else %do;
        4/llx=50  lly=0 
          ulx=50  uly=50 
          urx=75 ury=50 
          lrx=75 lry=0 
          color=white
        5/llx=75  lly=0 
          ulx=75  uly=50 
          urx=100 ury=50 
          lrx=100 lry=0 
          color=white
	%end;
	;
	template crb; 
	%if &score2= %then %do; 
	treplay 2:crbsel 3:crbperf 1:crbdisc 4:crbdens1; /* On réaffiche les graphes comme on en a envie */
	%end;
	%else %do;
	treplay 2:crbsel 3:crbperf 1:crbdisc 4:crbdens1 5:crbdens2;
	%end;
	quit; 
%end;

proc datasets nolist;
delete __N __tmp __tmp1 __tmp2 __f0 __f1; /* On efface les tables temporaires pour que ce soit propre */
quit;

title;
title2;

%fin:
%put Fin de la macro COURBES.;

options notes; /* On réactive l'affichage des notes dans la log */

%mend courbes; /* Fin de la macro */


/* Exemple d'exécution */



%courbes(data=fic.sco,score= p ,critere=y,courbes=1 1 1 1,taux=1);
quit;

quit;