//Scilab Code for the Microcredit Projet
//@date: June 2010
//@authors: Leo AUGE- Aurore LEBRUN - Anaïs PIOZIN


//////////////////////////////////
//          Model datas          //
//////////////////////////////////
clf();
somme = 1000; //loaned amount
remboursement = 22; // amount to pay back each week
n = 50; // number of pay back
e=52; // sample (week, month,year,...)
p = 0.8; //pay back probability



////////////////////////////////////////////////////////////
//          Calculation of the polynom solutions          //
////////////////////////////////////////////////////////////
funcprot(0);

// Calculate all the polynom solutions with the Yunus model
function [sols]=solutions()
  q=poly(0,"q");
  [sols]=roots(remboursement*q^(n+1)-(somme+remboursement)*q+somme); //
endfunction

//Test for this function
sols=solutions();


// Function that select the right solution among all the solutions (it is the real one between 0 and 1)
// sols is the solutions vector
function[solution]=solutionChoisie(sols)
  solution=[];
  imgSols = imag(sols); //imaginary part of the solutions
  for i=1:length(imgSols),
    if ( imgSols(i)==0 & real(sols(i))>0 & real(sols(i))<1-0.00000001 ) then //pick the real solution between 0 and 1
          solution=[solution,sols(i)];
    end,
  end,
  solution = real(solution);
endfunction,

//Test for this function
sol=solutionChoisie(sols);



///////////////////////////////////////
//          Rate calculation         //
///////////////////////////////////////

// Calculation of the rate (case where there is no lateness)
function [r]=taux()
  sols = solutions();
  solution = solutionChoisie(sols);
  r = -e*log(solution);
endfunction

//Test for this function
r = taux();

//disp() print the message or the element between the parentesis
disp("*****   Rate calculation example  *****");
disp("Applied rate (case where there is no lateness)");
disp(r);



///////////////////////////////////////////////////////
//          Effective average rate calculation       //
///////////////////////////////////////////////////////

// Calculate the effective average rate in two different ways
// The paramater p is the payback probability
function[r,r2]=taux_effectif_moyen(p)
  sols = solutions();
  y=solutionChoisie(sols);
  r=e*log((p+y*(1-p))/y);
  r2=e*p*(1/y-1);
endfunction

//Plotting of the effective average rate according to the payback probability
p1=[0:0.02:1];
[taux,r2]=taux_effectif_moyen(p1);
xset("window",0) //window in which the graph will be plot
xtitle("Effective average rate calculation depending on the payback probability"); //title of the window
xlabel("Payback probability"); // legend of the x axis
ylabel("Effective average rate"); //legend of the y axis
plot(p1,taux,'r'); // function which plot the calculated taux effectif moyen 
plot(p1,r2,"+"); // function which plot the theoric taux effectif moyen 
legend("Calculated effective average rate","Theoritical effective average rate",2); //legend of the graph



//////////////////////////////////////////////////////////////
//          Simulation of a geometric distribution          //
//////////////////////////////////////////////////////////////

// simulate a value according to the geometric law
function x = randGeometrique()
  y = rand(); // rand()=function which simulate equally a number between 0 and 1
  if (y<=p) then x = 1;
  else 
    k = log(1-y)/log(1-p);
    x = ceil(k);
  end,
endfunction,


// simulate a vector of value chosen according to the geometric law
// the parameter nb is the length of the wanted vector
function x = vecteurGeom(nb)
  x=[];
  for i=1:nb
    x(i) = randGeometrique(),
  end,
endfunction,



/////////////////////////////////////////////////////////////////////
//          Rate calculation with a geometric distribution         //
/////////////////////////////////////////////////////////////////////

// function which calculate the rate for a geometric distribution
// the parameter x is a vecteur sistributed according to the geometric law
// (use the function vecteurGeom to get an x)
function [r]=tauxSim(x)
  q=poly(0,"q"); //declaration of q as a variable
  y=q^x(length(x));
  for i=length(x)-1:-1:1
    y=q^x(i)*(1+y);
  end
  y=remboursement*y-somme;
  sols = roots(y);
  solution = solutionChoisie(sols);
  r = -e*log(solution);
endfunction

// Test of the previous function
disp("");
disp("*****   Distribution example   *****");
vect = vecteurGeom(n)';
s=sum(vect);
disp("Number of weeks to pay back");
disp(s);
rSim = tauxSim(vect);
disp("Interest rate for this simulation : ");
disp(rSim);



////////////////////////////////////////////////////////////////////
//          Histogram plotting for geometric simulations          //
////////////////////////////////////////////////////////////////////

//calculate lots of geometric distribution and the associated rate
t=[];
for i=1:1000
  x=vecteurGeom(n);
  r=tauxSim(x);
  t=[t,r];
end

// Plotting of these rates with an histogram
// (Figure n°1)
disp("");
disp("*****   Lots of distribution simulation   *****");
xset("window",1);
clf();
classe=min(t):0.001:max(t);
rb=taux_effectif_moyen(p)
disp("Effective average rate  : ");
disp(rb);
effectiveAverageRate=rb*ones(n,1);
y=0:1:n-1;
meanAverage = mean(t)*ones(n,1);
disp("Average of the calculated interest rate  : ");
disp(mean(t));

xtitle("Interest rate distribution")
xlabel("Interest rate");
ylabel("number of occurance");
histplot(classe,t, style=2);
plot(effectiveAverageRate,y,'r');
plot(meanAverage,y,"black");
ecartT = sqrt(variance(t));
normal = 1/(ecartT*sqrt(2*%pi))*exp((-1/2)*((classe-mean(t))/ecartT)^2);

plot(classe, normal,color=13);
legend("Interest rate distribution","Effective average rate","Mean of the distribution","Gaussian law",2);



  //////////////////////////////////////////////////////
 //          Testing of normal law properties        //
//////////////////////////////////////////////////////
disp("Third moment (must be equal to 0 to be a normal law) : ");
disp(moment(t,3));
disp("Kurtosis (must be equal to 3 to be a normal law) : ");
disp(moment(t,4)/(moment(t,2)^2));



/////////////////////////////////////////////////////////////////////
//          Calculation of the probability of a d-default          //
/////////////////////////////////////////////////////////////////////

// Calculate the probabilty to have a d-default
// The parameter p is the payback probability
// The parameter d is the "d" of "d-default"
function pi = default(prob,d)
  pi = (1-prob)^d;
endfunction

// Plotting of the probability to have a default according to the payback probability
// (Figure n°2)
pVar=[0:0.001:1];
xset("window",2);
clf();
xtitle("Default in function of the probability");
xlabel("Pay back probability");
ylabel("Probability of a default");
// uncomment the following line to have the 5% level on the graph
//plot(pVar,[default(pVar,2)' default(pVar,3)' default(pVar,5)' default(pVar,9)' 0.05*ones(1,length(pVar))']);
plot(pVar,[default(pVar,2)' default(pVar,3)' default(pVar,5)' default(pVar,9)']);
legend("2-default","3-default","5-default","9-default");
//xs2gif(2,"DefaultProbability")



/////////////////////////////////////////////////////////////////////////////////////
//          Calculation of the pay back probabilty to have a default of 5%         //
/////////////////////////////////////////////////////////////////////////////////////

// Calculate the payback probability to have a probability of def% to have a d-default
// The parameter def is the wanted d-default probability
// The parameter d is the "d" of "d-default"
function prob = payBack(def,d)
  prob = 1-(def)^(1/d)
endfunction,
disp("");
disp("*****   D-default   *****");
prob = payBack(0.05,5)
disp("Pay back probability to have a 2-default of 5% : ");
disp(prob)
rateDefault = taux_effectif_moyen(prob)
disp("Effective average rate for this probability : ");
disp(rateDefault);




/////////////////////////////////
////          Cetelem         //
///////////////////////////////
remboursementC = 103; //amount to pay back every month
sommeC = 1000;  //loaned amount
period = 12; //period to pay back :month
p=0.7 // probability to pay back on time


//function to simulate rate in the case of Cetelem (consumer credit companie)
// When you don't refund on time, they charge a penality rate on what you forget to pay.
// So when you want to refund, you have to pay what you forget to pay plus the penalities.
function [r]=reportSC()
  q=poly(0,"q");
  t = 0.01; //penality for not paying back on times
  x=vecteurGeom(period);
  y=0;
  refundedWeeks = 0;
  i=1;
  
  //Building of the polynomial
  while (refundedWeeks<10), // 10 because we wanted to refund during 10 month no more
    if (x(i)==1) then
      y = y+remboursementC*q^(sum(x(1:i)));  
    else
      rembR = remboursementC;
      for k=1:(x(i)-1)
        rembR = remboursementC + (1+t)*rembR;
      end,
      y = y+rembR*q^(sum(x(1:i)));
    end
    refundedWeeks = refundedWeeks + x(i);
    i = i+1;
  end
  y=y-sommeC;
  
  //find the rate
  sols=roots(y);
  solution = solutionChoisie(sols);
  r = -period*log(solution);
endfunction

//Test for this function
disp("");
disp("*****   Consumer credit   *****");
cetelem=reportSC();
disp("Calculated rate");
disp(cetelem);

////////////////////////////////////////////////////////////////////
//          Histogram plotting for consumer credit model          //
////////////////////////////////////////////////////////////////////

//function calculate several rates in case of Cetelem
consumerRate=[];
for i=1:10000
  consumerRate = [consumerRate,reportSC()];
end,


//Building of a Histogramm in the case of Cetelem
// (Figure n°4)
xset("window",4);
clf();
classe=min(consumerRate):0.001:max(consumerRate);
rb=taux_effectif_moyen(p)
effectiveAverageRate=rb*ones(1.15*n,1);
yc=0:1:1.15*(n-1);
meanAverage= mean(consumerRate)*ones(1.15*n,1);
disp("Average of the calculated interest rates  : ");
disp(mean(consumerRate))

xtitle("Interest rate for a consumer credit model")
xlabel("Interest rate");
ylabel("number of occurance");
histplot(classe,consumerRate, style=2);

plot(meanAverage,yc,"black");
legend("Interest rates","Mean of the distribution",1);


  x=vecteurG