/////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////                                                                                                 //////
///// Processus/Diffusions à champ moyen : Modèle de Stabilité interbancaire J.P. Fouque  et L.H. Sun  //////
/////                                                                                                  //////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////

clf();
stacksize("max");


T=1;
n=100; //Pas de temps T/n
N=10;  //Nombre de banques (N T/n<1/10)

D=-0.7; //barrière de défaut

R=1000; // Répétition de la simulation (<10^4 pour éviter trop d'attente)

//Volatilité

sigma=1;

ND_X=[];
ND_Y=[];


// Taux de prêt interbancaire autorisé (ajusté par les autorités interbancaires pour augmenter la stabilité 
// du système)

alpha=10; // alpha=0, 1, 10, 100




for r=1:R


dW=grand(n,N,"nor",0,sqrt(T/n));
time_mesh=(T/n)*[0:n]';



X=zeros(1,N);
Y=X;


// Taux de réserve monétaire des N banques (placements, investissements, etc)

delta_min=0;
delta_max=.1;
delta=grand(1,N,"unf",delta_min,delta_max);  

//delta=zeros(1,N)




for t=1:n

//// Evolution (prêts/emprunts) des réserves (log-)monétaires X(t,i) de i=1:N banques 
////  à chaque instant t,au taux normalisé alpha/N

X=[X;X(t,:)+delta*(T/n)+alpha*(mean(X(t,:))-X(t,:))*(T/n)+sigma*dW(t,:)];

//// Evolution sans corrélations  

Y=[Y;Y(t,:)+delta*(T/n)+sigma*dW(t,:)];
end




if r==1 then


clf(0);xset("window",0);show_window();

 plot2d(time_mesh,D*ones(1:n+1),5)

for NN=1:N
 plot2d(time_mesh,X(:,NN),NN)
end


clf(1);xset("window",1);show_window();

 plot2d(time_mesh,D*ones(1:n+1),5)
for NN=1:N
 plot2d(time_mesh,Y(:,NN),NN)
end


end

//Nombre de défauts observés

Defauts_X=double(X<D);
Defauts_Y=double(Y<D);

Nb_Defauts_X=sum(Defauts_X,'r');
Nb_Defauts_Y=sum(Defauts_Y,'r');



ND_X=[ND_X,sum((Nb_Defauts_X>0))];
ND_Y=[ND_Y,sum((Nb_Defauts_Y>0))];


end



//Estimation des probabilités d'avoir i défauts 
Proba_X=zeros(1:N);
Proba_Y=zeros(1:N);


for i=1:N
p=double((ND_X==i)); 
q=double((ND_Y==i));    
Proba_X(i)=sum(p)/R;
Proba_Y(i)=sum(q)/R;
end

clf(2);xset("window",2);show_window();
subplot(121);
plot(Proba_X);
subplot(122);
plot(Proba_Y);