// multiplier.cpp
//
// This program implements the model described in "The Microfoundations of the
// Keynesian Multiplier Process," by Peter Howitt, to be published in the
// inaugural issue of the Journal of Economic Interaction and Coordination.
// It is written for the Borland C++Builder 5 environment.
//
// January 4, 2005

#pragma hdrstop
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <conio.h>

#define n 50                                   // Number of goods
#define bsize 3                                // Number of each type of person
#define K 80                                   // Number of store locations
#define xMax n                                 // Maximal target of a new entrant
#define lambda 1.0   		    	       // Autonomous search intensity
#define alpha 0.75			       // Adjustment speed
#define theta 0.03      		       // Weekly failure rate of unprofitable stores
#define C 2.0				       // Setup cost in flow terms
#define f 6.0                                  // The overhead cost
#define numruns 10000                          // Number of runs for each configuration
#define gamma .12                              // Size of shock
#define T 250                                  // Number of weeks in each run
#define RAND 7                                 // Initial seed for each new configuration
#define tol 0.2*f                              // How much overhead the owner is willing to cover

#define N (bsize*n*(n-1))		       // Number of people - bsize of each type
#define NC (bsize*(n-1)*(n-1))                 // Number of non-money consumers
#define F(x,y,z) ((y-z>0)?((y-z)/x):0)         // Pricing formula
#define rnd01 ((1+random(1000))/1000.0)
#define capacity (N-(n-1)*f)

int cu[2], willing_cu[2];
int s[N+1], d[N+1];			       // s[r] d[r] are r's supply and demand goods
int store[N+1];				       // Trader r buys from store[r]
int employer[N+1];                             // Trader r sells to employer[r]
int NS;					       // Number of stores
int active[K+1];			       // 0 if inactive shop location
int g[K+1];     			       // Good traded at shop k (along with commodity 1)
int numprod[n+1];
int numcons[n+1];
int producer[n+1][N+1];                        // producer[i][j] is the jth producer of i
int consumer[n+1][N+1];                        // consumer[i][j] is the jth consumer of i
int owner[K+1];                                // The individual that owns shop k
int familyshop[N+1];                           // The shop owned by r
int h, i, j, k, a, b;                          // Indices
int r, fr;                                     // Used to index individuals and friends
int t;                                         // Used to index week
int run;                                       // Used to index the run
int line[N+1];                                 // Random lineup before each week's search
double W, P;                                   // Trial prices in Research()
double Tq, Tm;                                 // Targets established in Research()
double trm[K+1], trq[K+1];		       // Target values for m and q
double U;               		       // utility indices
double part;				       // number of participants
double w[K+1];				       // wholesale prices
double p[K+1];                                 // inverse-retail prices
double q[K+1];				       // q[k] = supply of input to shop k
double m[K+1];                                 // m[k] = money income of shop k
double Csurp, Psurp;                           // Consumer and producer surpluses
double AvGDP[T+1];                             // Average GDP after t weeks
int numfails[numruns+1];                       // Number of shop failures in each run
double minGDP[numruns+1];                      // Minimum GDP gap in each run
double varshock[numruns+1];                    // Variance of shock
double assym[numruns+1];                       // Assymetry of shock
time_t begin, finish, firstbegin;
FILE *stream;

void Init(void);
void InitRun(void);
void lineup(void);
int Research(void);
int comrade(void);
int soulmate(void);
double Utility(void);
int Calc1(void);
void Calc2(void);
void shock(void);
double var(void);
double skew(void);

main( ) {
  Init();
  for(run=1;run<=numruns;run++)  {
    InitRun();
    for(t=0;t<=T;t++) {

                // Entry
   for(b=1;b<=n/5;b++) {
      do {
        r=random(numcons[1])+1; r=consumer[1][r];
        if( NS<K && familyshop[r]==0) {
          j=Research();
          if(j>0) {
            NS++; k=1; while(active[k]==1) k++;
            active[k]=1; g[k]=s[r];
            trq[k]=Tq; trm[k]=Tm;
            p[k]=F(trm[k],trq[k],f);
            w[k]=F(trq[k],trm[k],C);
            employer[r]=k; familyshop[r]=k; owner[k]=r;
            for(h=0;h<=1;h++) if(willing_cu[h]==1) store[cu[h]]=k;
          }
        }
      }
      while (NS==0);
   }

                 // Matching
      for(i=1;i<=N;i++)  {
        r=line[i];
        U=Utility();
        if(rnd01<( U>0 ? lambda : lambda ) && familyshop[r]==0) {
          k=random(K)+1;
          if(s[r]!=1) {
            fr=comrade();
            if(w[employer[fr]]>w[employer[r]]) employer[r]=employer[fr];
            if(g[k]==s[r] && w[k]>w[employer[r]]) employer[r]=k;
          }
          if(d[r]!=1) {
            fr=soulmate();
            if(p[store[fr]]>p[store[r]]) store[r]=store[fr];
            if(g[k]==d[r] && p[k] > p[store[r]]) store[r]=k;
          }
        }
      }

                 // Trade
      for(k=0;k<=K;k++) if(active[k]>0) {q[k]=0; m[k]=0;}
      for(r=1;r<=N;r++) {
        if(w[employer[r]]>0 && (p[store[r]]>0 || d[r]==1)) q[employer[r]]++;
        if(p[store[r]]>0) {
          if(employer[r]>0) m[store[r]]+=w[employer[r]];
          if(s[r]==1) m[store[r]]++;
        }
      }

                // Exit
      for(k=1;k<=K;k++) if(active[k])
      if( ( q[k]-p[k]*m[k]-f < -tol ||
            m[k]-w[k]*q[k] < 0 ) && rnd01<=theta ) {
        active[k]=0; NS--; familyshop[owner[k]]=0; owner[k]=0;
        q[k]=0; m[k]=0; g[k]=0; trq[k]=0; trm[k]=0; p[k]=0; w[k]=0;
        for(r=1;r<=N;r++) {
          if(store[r]==k) store[r]=0;
          if(employer[r]==k) employer[r]=0;
        }
        numfails[run]++;
      }

                // Update targets and prices
      for(k=1;k<=K;k++) if(active[k]) {
        trq[k]+=alpha*(q[k]-trq[k]);
        trm[k]+=alpha*(m[k]-trm[k]);
        w[k]=F(trq[k],trm[k],C);
        p[k]=F(trm[k],trq[k],f);
      }

                   // Periodic check during run
      Calc1();
      Calc2();
     // printf("%6.0f ", part );
     // printf("%4d %6.0f %6d\n", NS, Psurp + Csurp, t);
      AvGDP[t]+=((Psurp + Csurp)/capacity - AvGDP[t])/run;
      minGDP[run] = min((Psurp + Csurp)/capacity, minGDP[run]);
      if(t==0) {
        shock();
        varshock[run]=var();
        assym[run]=skew();
      }

            // End of the week loop
    }
    finish = clock();
    printf("Run %4d  seconds=%.2f ", run, (finish - begin)/CLK_TCK);
    printf(" var=%8.5f skew=%9.5f ", varshock[run], assym[run]);
    printf(" fails=%4d  MinGDP=%7.5f\n", numfails[run], minGDP[run]);

            // End of the run loop
  }

  /*stream=fopen("skew.txt", "w+");
  fprintf(stream, "n   bsiz K  xMax lambda alpha  theta C    f   Nrun gamma   T     RND  tol\n");
  fprintf(stream, "%2d %2d   %3d %3d %6.3f %6.3f %6.3f %3.1f  %2.1f %2d    %2.2f  %5d %3d   %2.1f\n", n, bsize, K, xMax,
           lambda, alpha, theta, C, f, numruns, gamma,  T,  RAND, tol);
  fprintf(stream, "\n\nWeek      Fail        MinGDP     varshock     assymetry");
  for(run=1;run<=numruns;run++) {
    fprintf(stream, "\n %3d   %7d       %5.5f     %5.5f  %11.5f", run, numfails[run], minGDP[run], varshock[run],
             assym[run]);
  }
  fclose(stream);
  */

  time(&finish);
  stream=fopen("AvGDP.txt", "w+");
  fprintf(stream, "\n Started : %s", ctime(&firstbegin));
  fprintf(stream, " Finished: %s\n\n", ctime(&finish));
  fprintf(stream, " week    AvGDP\n");
  for(t=0;t<=T;t++) fprintf(stream, "%5d    %6.5f\n", t, AvGDP[t]);
  fclose(stream);
  printf("\n Started : %s", ctime(&firstbegin));
  printf(" Finished: %s\n\n", ctime(&finish));
  printf("Press any key to end the program...\n\n");
  getch();
  return(0);
}

void Init(void) {
  firstbegin=clock();
  time(&firstbegin);
  familyshop[0]=0; d[0]=0; s[0]=0; store[0]=0; employer[0]=0;
  for(t=0;t<=T;t++) AvGDP[t]=0;
}

void InitRun(void) {
  // printf("Number  \n");
  // printf("Active   NS    GDP   Week   \n");
  srand(RAND+run-1);
  begin = clock();
  for(k=0;k<=K;k++) {
    active[k]=0; owner[k]=0;
    p[k]=0; w[k]=0; trm[k]=0; trq[k]=0; g[k]=0;
  }
  for(i=0;i<=n;i++) {
    numprod[i]=0; numcons[i]=0;
    for(j=0;j<=N;j++) {producer[i][j]=0; consumer[i][j]=0;}
  }
  r=0;
  for(i=1;i<=n;i++) for(j=1;j<=n;j++) if(i != j) {
    for(k=1;k<=bsize;k++) {
    r++; s[r]=i; d[r]=j; familyshop[r]=0;
    numprod[i]++; producer[i][numprod[i]]=r;
    numcons[j]++; consumer[j][numcons[j]]=r;
    }
  }
  NS=n-1;
  for(i=1;i<=n-1;i++) {
    active[i]=1; owner[i]=producer[i+1][1];
    g[i]=i+1; familyshop[owner[i]]=i;
    w[i]=1-C/bsize;
    trm[i]=bsize*(n-1)-C*(n-2);
    trq[i]=(n-1)*bsize;
    p[i]=F(trm[i],trq[i],f);
    for(h=1;h<=(n-1)*bsize;h++) {
      r=producer[i+1][h];  employer[r]=i;
      r=consumer[i+1][h];  store[r]=i;
    }
  }
  lineup();
  numfails[run]=0;
  minGDP[run]=1.0;
}

void lineup(void) {
  int i,j,k;
  int start[N+1];
  line[0]=0;
  for(i=0;i<=N;i++) start[i]=i;
  for(j=1;j<=N;j++) {
    k=random(N-j+1)+1;
    line[j]=start[k];
    for(i=k;i<=N-j;i++) start[i]=start[i+1];
  }
}

int Research(void) {
  for(h=0;h<=1;h++) {
    cu[h]=0; willing_cu[h]=0;
  }
  Tm=random(xMax)+1; Tq=random(xMax)+1;
  W = F(Tq,Tm,C);
  if(W<=w[employer[r]]) return(0);
  P = F(Tm,Tq,f);
  for(h=0;h<=1;h++) {
    k=random(numcons[s[r]])+1;
    cu[h]=consumer[s[r]][k];
    if(P>p[store[cu[h]]]) willing_cu[h]=1;
  }
  if(willing_cu[0]+willing_cu[1]==0) return(0);
  return(1);
}

double Utility(void) {
  double X=0;
  if(d[r]==1) X=w[employer[r]];
  if(s[r]==1) X=p[store[r]];
  if(d[r]!=1 && s[r]!=1) X=w[employer[r]]*p[store[r]];
  return(X);
}

int comrade(void) {
  int k=1+random(numprod[s[r]]-1);
  int fr= producer[s[r]][k];
  if(fr>=r) fr=producer[s[r]][k+1];
  return(fr);
}

int soulmate(void) {
  int k=1+random(numcons[d[r]]-1);
  int fr= consumer[d[r]][k];
  if(fr>=r) fr=consumer[d[r]][k+1];
  return(fr);
}

int Calc1(void) {    // The calculations of number of active participants
  part=0;
  for(r=1;r<=N;r++) {
    U=Utility();
    if(U>0 || familyshop[r]>0) part++;
  }
  return (1);
}

void Calc2(void) {   // The calculations made for GDP.
  Csurp=0;
  for(r=1;r<=N; r++) Csurp += Utility();
  Psurp=0;
  for(k=1;k<=K;k++) if(active[k]) {
    Psurp += m[k]-w[k]*q[k];
  }
}

void shock(void) {
  int i,j,k;
  int start[NC+1];
  for (k=0;k<=NC;k++)start[0]=0;
  k=1;
  for(i=2;i<=n;i++) for(j=1;j<=numcons[i];j++) {
    start[k]=consumer[i][j];
    k++;
  }
  for(k=1;k<=gamma*NC;k++) {
    i=random(NC-k+1)+1;
    j=start[i];
    store[j]=0;
    for(j=i;j<=NC-k;j++) start[j]=start[j+1];
  }
}

double var(void) {
  int i,j;
  double average=0, variance=0;
  double x[n];
  for(j=0;j<=n;j++) x[j]=0;
  for(i=1;i<=N;i++) x[d[i]-1]+=(store[i]==0);
  for(j=1;j<=n-1;j++) average+=x[j];
  average=average/(n-1);
  for(j=1;j<=n-1;j++) variance+=(x[j]-average)*(x[j]-average);
  variance=variance/(n-1);
  return(variance);
}

double skew(void) {
  int i,j;
  double average=0, skewness=0;
  double x[n];
  for(j=0;j<=n;j++) x[j]=0;
  for(i=1;i<=N;i++) x[d[i]-1]+=(store[i]==0);
  for(j=1;j<=n-1;j++) average+=x[j];
  average=average/(n-1);
  for(j=1;j<=n-1;j++) skewness+=(x[j]-average)*(x[j]-average)*(x[j]-average);
  skewness=skewness/(n-1);
  return(skewness);
}