/*GAUSSI Programme : RNDWALK = Random Walk Process*/
/*Monte-Carlo Study for Random Walk Process: March 28, 1995   revised September 8,1998*/
/*Reexamination of William Feller's Model*/
/*Systemic Risk Simulation*/
/*Programme with Capital Injection (back to position 0) Considered*/
/*Seed not change after injection*/

/********************* 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 */
R = zeros(n,k);              /* record amount of capital injection in each round and bank */
R2 = zeros(n+1,k);          /* record of accumulated injection history */

/********************* 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;

/*  if R2[t,j] <= 0;*/
     seed1 = j;
     Z = rndus(sim,n,seed1)-0.5;   /* uniform random variable (-0.5,0.5)*/
/*  elseif R2[t,j] > 0;
     seed1 = j+1000;
     Z = rndus(sim,n,seed1)-0.5;
  endif;*/

    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;

/*********** Start the zero-sum operation *******************/
a=0;
i=1;
do while i <= k;
    a = a + X[sim+1,i];
i=i+1;
endo;

i=1;
do while i <= k;
    X[sim+1,i] = X[sim+1,i] - a/k;
i=i+1;
endo;

/************  Injection of Capital **************/
i=1;
inj=0;
do while i <= k;

    if X[sim+1,i] + W[t,i] < 0;
       inj = -(X[sim+1,i] + W[t,i]);
       X[sim+1,i] = X[sim+1,i] + inj; /*capital inj injected*/

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

    R[t,i] = inj;
    R2[t+1,i] = R2[t,i] + R[t,i];
i=i+1;
endo;

/***************** Payoff Recharged *************************/
i=1;
do while i <= k;
    W[t+1,i] = W[t,i] + X[sim+1,i];
i=i+1;
endo;

t=t+1;
endo;

/******** Graphics *********************/
/*M = zeros(n+1,k); /* matrix representing X-axis in the graph */

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;
output on;
W;
output off;

library pgraph;
graphset;

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

XY(X,Y);*/

/******** How many returns 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,100);*/

/***** How long can leads or lags would be sustained ? *******/
/*L = zeros(k,1);
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:
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);*/

/***** Record of capital injection and identify who needed the capital *****/
/*output file = a:R.fmt;
output on;
R;
output off;*/

a=0;  /* how much capital injected */
  i=1;
  do while i <= k;
      a=a+R2[n+1,i];
  i=i+1;
  endo;

b=0; /* how many banks injection taken place */
  i=1;
  do while i <= k;
      if R2[n+1,i] > 0;
         b=b+1;
      elseif R2[n+1,i] <= 0;
         b=b;
      endif;
  i=i+1;
  endo;

/*print R;*/
print "amout of injection = " a;
print "# of banks received injection = " b;

/** Graphics 2-2 : Distribution of Asset Level **/
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 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);