/*GAUSSI Programme : RNDWALK = Random Walk Process*/
/*Monte-Carlo Study for Random Walk Process: March 28, 1995   revised September 2,1998*/
/*Reexamination of William Feller's Model*/
/*Systemic Risk Simulation */

/********************* Defining Key Parameters *************************/
sim= 100;	 			/* number of trials */
k= 100;					/* number of banks or players */
v= 151;						/* initial average capital endowment */
n= 100;                    		/* how many round of play */

/********************** Defining Variable Size  *************************/
W = ones(n+1,k)*v;		/* equal initial capital endowment */

/********************* Creating Random Variables ************************/
t=1;
do while t <= n;

X = zeros((sim+1),k);       /* make these matrices back to zero every time t counts on */
Y = zeros(sim,k);
Z = zeros(sim,n);

  j=1;
  do while j <= k;

  seed1=j;
  Z=rndus(sim,n,seed1)-0.5;   /* uniform random variable (-0.5,0.5)*/
/*format /m1 /ldn 5,3;*/
/*print "rnd-matrix =" Z;*/

    i=1;
    do while i <= sim;

    if Z[i,t]>=0;
       Y[i,j]=1;
    elseif Z[i,t]<0;
       Y[i,j]=-1;
    endif;

    X[i+1,j]=X[i,j]+Y[i,j];	/* Random Walk Process */

    i=i+1;
    endo;

  j=j+1;
  endo;

/*format /m1 /ldn 5,3;*/
/*print X;*/

/**************  Exclude players whose previous position is negative **************/
a=0;
c=0;

i=1;
do while i <= k;

    if W[t,i] < 0;
       X[sim+1,i]=0;
       c = c+1;
    elseif W[t,i] >= 0;
       a = a+X[sim+1,i];
    endif;

i=i+1;
endo;

/************* Start the zero-sum operation among survivors *******************/
i=1;
do while i <= k;

    if W[t,i] >= 0;
       X[sim+1,i]=X[sim+1,i] - a/(k-c);
    endif;

i=i+1;
endo;

/** Winner and Loser Specified **/
win=0;
loss=0;

i=1;
do while i <= k;

if W[t,i] < 0;
   X[sim+1,i] = X[sim+1,i];

elseif W[t,i] >= 0;
   if X[sim+1,i] + W[t,i] >= 0;

     if X[sim+1,i] < 0;
        loss = loss - X[sim+1,i];
     elseif X[sim+1,i] >= 0;
        win = win + X[sim+1,i];
     endif;

   elseif X[sim+1,i] + W[t,i] < 0;
     loss = loss + W[t,i];
   endif;
endif;
i=i+1;
endo;

/*print "# of winner =" win;
print "# of loser =" loss;*/

/** Payoff Recharged  **/
i=1;
do while i <= k;

if X[sim+1,i] + W[t,i] >= 0;
     if X[sim+1,i] < 0;
        W[t+1,i] = W[t,i] + X[sim+1,i];
     elseif X[sim+1,i] >= 0; 
        W[t+1,i] = W[t,i] + loss * X[sim+1,i] / win;
     endif;
elseif X[sim+1,i] + W[t,i] < 0;
     W[t+1,i] = -100;
endif;

i=i+1;
endo;

t=t+1;
endo;

/************ count the number of default *******************************/
def = 0;

i=1;
do while i <= k;
    if W[n+1,i] < 0;
       def = def+1;
    endif;
i = i+1;
endo;

/*format /m1 /ldn 3,1;*/
/*print "Asset-matrix =" W;*/

print "# of failure =" def;

/****** Graphics 1 : Change in each bank's asset position ***********/
M = zeros(n+1,k);  /* make up a matrix to represent x-axis */

t=1;
do while t <= n+1;
  j=1;
  do while j <= k;
      M[t,j] = t;
  j=j+1;
  endo;
t=t+1;
endo;

output file = a:M.fmt;
output on;
M;
output off;

output file = a:W.fmt;       /* save W in the file */
output on;
W;
output off;   /* note W and M must be equal size */

library pgraph;
graphset;

load X[n+1,k]=a:M.fmt;
load Y[n+1,k]=a:W.fmt;

 xlabel("Number of Rounds");
 ylabel("Asset Level");
 cmdstr = "-p" $+ "-pr=3" $+ "_pline={2 6 0 0 0 0 0 9 10}";
 graphprt(cmdstr);
 XY(X,Y);

/***** Graphics 2 : How many come backs to the level of original endowment ? **/
L = zeros(k,1);   /* matrix to record times of return to the original endowment */

i=1;
do while i <= k;
a=0;
   t=1;
   do while t <= n;
       if (W[t,i]-v)*(W[t+1,i]-v) <= 0;
           a=a+1;
      elseif (W[t,i]-v)*(W[t+1,i]-v) > 0;
           a=a;
      endif;
   t=t+1;
   endo;
 L[i,1] = a;
i=i+1;
endo;

output file = a:L.fmt;
output on;
L;
output off;

library pgraph;
graphset;

load X[k,1]=a:L.fmt;
{C,M,F}=hist(X,50);

/** Graphics 2-2 : Distribution of Asset Level (Normal Dist.)**/
/*L = zeros(k,1);   /* matrix to record times of return to the original endowment */

t=1;
do while t <= k;
    L[t,1] = W[n+1,t];
t=t+1;
endo;

output file = a:L2.fmt;
output on;
L;
output off;

library pgraph;
graphset;

load X[k,1]=a:L2.fmt;
{C,M,F}=hist(X,50);*/

/***** Graphics 3 : How long can leads or lags be sustained? ******/
/*P = zeros(k,1);   /* matrix to record times of return to the original endowment */
i=1;
do while i <= k;
   a=0;
   t=0;
   do while t <= n-1;
       if (W[n-t,i]-v)*(W[n-1-t,i]-v) > 0;
           a=a+1;
      elseif (W[n-t,i]-v)*(W[n-1-t,i]-v) <= 0;
           goto fip;
      endif;
   t=t+1;
   endo;
fip:
P[i,1] = a;
i=i+1;
endo;

output file = a:P.fmt;
output on;
P;
output off;

library pgraph;
graphset;

load X[k,1]=a:P.fmt;
{C,M,F}=hist(X,50);*/

/***** Graphics 4 : How many leads ? arc-sine graph******/
/*P = zeros(k,1);   /* matrix to record times of leads */
i=1;
do while i <= k;
   a=0;
   t=2;
   do while t <= n+1;
       if W[t,i] > v;
           a = a+1;
      elseif W[t,i] <= v;
           a = a;
      endif;
   t=t+1;
   endo;
P[i,1] = a/n;
i=i+1;
endo;

output file = a:P2.fmt;
output on;
P;
output off;

library pgraph;
graphset;

load X[k,1]=a:P2.fmt;
{C,M,F}=hist(X,50);*/