// Aug 28, 2006.  The program generating results on wage-flexibility in "Macroeconomics
// with Autonomous, Intelligent Agents," to be presented by Peter Howitt at the Festschrift
// for Axel Leijonhufvud at UCLA, August 30, 2006.

#pragma hdrstop
#pragma warning(disable:4996)
#pragma warning(disable:4101)
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <conio.h>

const int numruns = 30;							// Number of runs for each configuration
const int RAND = 1;							    // Initial seed for each new configuration
const int T = 5000;								// Maximal number of weeks in each run
const int RptPer = 1;							// Reporting Period (in years)
const int n = 50;      							// Number of goods
const int m = n*(n-2);							// Number of people - 1 of each type
const int K =2*n;								// Number of store locations
const double F = 2.0;							// Size of fixed cost
const double capacity = m - n*F;				// Capacity GDP
const int xMax = 100;							// Maximal Salestarget of a new entrant
const int I = 2;								// Expected Frequency of innovation opportunities
const double beta = 0.01;						// Fraction of matches that breakup each week

const double kappa=1.0;							// desired inventory-sales ratio
const double ksi = 0.05;						// speed of adjustment for inventories
const double alf = 0.2; 						// Adjustment speed for shops, households and policy makers
const double theta = 0.01; 						// Failure rate of unprofitable stores
const double eps = 20.0;						// demand parameter (alpha / (1-alpha))
const double rho = 2/eps;						// Average markup of entrants
const double lambda = .5;						// Search propensity
const double v = 0.005;							// Propensity to spend out of wealth
const double liq = 0.25;						// Liquidity value of permanent income
const double kk = 8.0;							// Upper threshold on money holdings

const double tauMax = 0.4;						// Maximal tax rate
const double bstar = 5;							// target debt/weekly GDP ratio
const double bp = 1.0;							// slope of Taylor Rule
const double pistar = 0.0;						// the (weekly) inflation target

const double moneysupply=100.0;					// Initial Money supply
const double bondsupply=0.0;			  		// Initial Bond supply
const double taxrate=v*bstar;					// Initial taxrate
const double debtratio=bstar;					// Initial bond/income ratio
const double bondprice=1/(v+pistar);			// Initial bond price

#define rnd01 ((1+random(1000))/1000.0)
int s[m+1], d[m+1];								// s[r] is r's labor and d[r], 1+d[r] his demand goods
int employer[m+1];								// Trader r sells to employer[r]
int store[2][m+1];								// Trader r buys d[r] from store[0,r]
int NS;											// Number of stores
int active[K+1];								// 0 if inactive shop location
int g[K+1];  									// Non-money good traded at shop k
int produces[n+1][n-1];							// Names of those endowed with good i
int consumes[n+1][n-1];							// Names of those with i as primary consumption good
int numprod[n+1];								// Number of those endowed with i
int numcons[n+1];								// Number of those liking i
int owner[K+1];									// The individual that owns shop k
int familyshop[m+1];							// The shop owned by r
int h, i, j, q;
int k;											// Used to index firms
int r, fr;										// Used to index individuals and friends
int t;											// Used for week
int run=0;										// Used to index the run
int a, b[2];									// Used to index firms in search and Utility
int line[m+1];									// Random lineup before each week's search
int numfail;									// Number of failures 
double W[K+1];									// Wholesale prices
double Pm[K+1];									// Inverse Retail prices.
double Sales[K+1], Input[K+1];					// Actual Sales and Input
double trSales[K+1], trInput[K+1];				// Target Sales and Input
double Orders[K+1], virtInput[K+1];				// Virtual Sales (Orders) and Input
double mu[K+1];									// Store k's markup
double mkup;									// Provisional markup
double Pmeffect[2][m+1], Weffect[m+1];			// effective values of money and endowments
double S[K+1];									// Shop k's stocks of g[k]
double fi, X, Y, Z;
double GDP, P, w;								// Real GDP, Price level and real wage
double BS;										// Stock of bonds
double MS;										// Stock of fiat money
double Debt;									// Bonds plus money
double Deficit;									// Change in Debt
double PercentDeficit;							// Like it sounds
double tau;										// Sales tax rate
double by;										// bond to income ratio
double Wagg;									// Average wage rate
double Employment;
double speed;									// Speed of wage adjustment
   
time_t begin, finish, firstbegin;
FILE *stream;

double M[m+1];									// Money holdings
double B[m+1];									// Bond holdings
double D;										// Bond purchases
double Mstar;									// Target money holdings

double Yp[m+1];									// Nominal (weekly) permanent incomes
double E[m+1];									// Planned nominal expenditures
double Pb;										// Bond price				
double Pstar;									// Price-level target (moving)

double AvGDP[T/50+1];							// Cross-run average GDP at the end of year t/50
double AvFail[T/50+1];							// Cross-run average number of failures in year t/50
double MeanGDP[numruns+1];						// Time-average GDP
double MeanFail[numruns+1];						// Time-average number of failures

void Init(void);
void InitRun(void);
void lineup(void);
int comrade(void);
int soulmate(void);
int Calc1(void);
void setprices(void);
double Dem(double, double);
int random(int);
int Research(double);
double TaylorRule(double, double);
double FiscalRule(double, double, double);
void SetupShop(int, int, double, double);

int main(void) {
	Init();
	for(speed=0.5;speed<=3.0;speed+=0.5) {
		for(run=1;run<=numruns;run++) {
			InitRun();
			for(t=1;t<=T;t++) {		    
				// Entry
				do {
					for(r=1;r<=m;r++) {
						if(random(m)<I && NS<K && familyshop[r]==0 ) {
							X=random(xMax)+1;
							mkup=2*rho*rnd01;
							if(Wagg*(mkup*X-F)>Yp[r] && (Research(mkup)) ) {
								NS++;
								k=1;
								while(active[k]==1) k++;
								SetupShop(k,r,X,mkup);								
							}
						}
					}
				}
				while (NS==0);
				
				// Matching
				for(i=1;i<=m;i++) {
					r=line[i];
					if(rnd01<=(Weffect[r] ? lambda : 1) && familyshop[r]==0) {
						for(h=0;h<=1;h++) {
							fr=comrade(); a=employer[fr];
							if( virtInput[a]<=trInput[a] ) {
								if(Weffect[fr] > Weffect[r] && familyshop[fr]==0) {
									employer[r] = a;
									Weffect[r] = Weffect[fr];
								}
								if(Weffect[r]<W[a] && owner[a]==fr) {
									employer[r] = a;
									Weffect[r] = W[a];
								}
							}
						}
					}
					if(rnd01<=((Pmeffect[0][r] * Pmeffect[1][r]) ? lambda : 1))
						for(h=0;h<=1;h++) {
							fr=soulmate();
							for(q=0;q<=1;q++) if(Pmeffect[q][r] < Pmeffect[q][fr]) {
								store[q][r] = store[q][fr];
								Pmeffect[q][r] = Pmeffect[q][fr];
							}
						}
					k=random(K)+1;
					if( g[k]==s[r] && familyshop[r]==0 && Weffect[r] <= W[k]) {
						if(virtInput[k]<=trInput[k]) {
							employer[r] = k;
							Weffect[r] = W[k];
						}
					}
					for(q=0;q<=1;q++) {
						j=1+((d[r]-1+q)%n);
						if( (g[k]==j) &&  (Pmeffect[q][r]<=Pm[k]) ) {
							store[q][r] = k;
							Pmeffect[q][r] = Pm[k];
						}
					}
				}
	     
				// Policy decisions
				GDP=0; for(k=1;k<=K;k++) if (active[k]) GDP+=(Input[k]>F?Input[k]-F:0);
				X=0;  for(k=1;k<=K;k++) if (Pm[k] && Input[k]>F) X+=(Input[k]-F)/Pm[k];
				if(GDP*X>0) P=X/GDP;
				Pstar = Pstar * exp(pistar);
				Pb = TaylorRule(X,Pb);
				X=0; for(k=1;k<=K;k++) if (active[k]) X+=(Sales[k]/Pm[k]);
				Debt=MS+Pb*BS;
				BS=0; MS=0; for(i=1;i<=m;i++) { BS += B[i]; MS += M[i]; }
				Deficit=MS+Pb*BS-Debt;
				if(Debt>0) PercentDeficit = Deficit / Debt;
				if(X>0) by = Pb*BS/X;
				tau = FiscalRule(by, PercentDeficit, tau);
	     
				// Bond market
				for(i=1;i<=m;i++) M[i] += B[i];
				for(i=1;i<=m;i++) {
					if(familyshop[i]==0) {
						Yp[i] += alf * (Weffect[i] - Yp[i]); 
						E[i] =  v * Pb * ( B[i] + liq*Yp[i] );
						E[i] = __max( E[i], 0.0);
						Mstar = (kk/2)*E[i];
					}
					else {
						k=familyshop[i];
						Yp[i] =  ( (1+mu[k]) * trSales[k]  -  trInput[k] + 1  ) * W[k];
						E[i] =  v * Pb * ( B[i] + liq*Yp[i] );
						E[i] = __max( E[i], 0.0);
						Mstar = (kk/2)*(E[i] +  W[k]*(trInput[k]-1));
					}
					if(M[i]<(2/kk)*Mstar || M[i] > 2*Mstar) {
						D = (M[i]-Mstar)/Pb;
						D = __max( D, -B[i] );
						D = __min( D, M[i] / Pb );
						M[i] -= Pb * D;
						M[i] = __max(M[i],0.0);
						B[i] += D;
						B[i] = __max(B[i],0.0);
					}
				}
	       
				// Trade and Adjust stocks
				for(k=1;k<=K;k++) if(active[k]>0) {
					Sales[k]=0; Input[k]=0; Orders[k]=0; virtInput[k]=0;
				}
				for(i=1;i<=m;i++) {
					r=line[i];
					a=employer[r]; for(h=0;h<=1;h++) b[h]=store[h][r];
					if(a>0) virtInput[a]+=1;
					if(owner[a]==r) {
						Input[a]++; S[a]++;
					}
					else {
						Weffect[r] = __min( W[a], M[owner[a]]);
						if(Weffect[r]>0) {
							Input[a]++; M[r] += Weffect[r]; M[owner[a]] -= Weffect[r]; S[a]++;
						}
					}
					E[r] = __min( E[r], M[r] );
					for(h=0;h<=1;h++) if(b[h]>0) {
						X=Dem( Pm[b[h]],Pm[b[1-h]] ) * E[r];
						Orders[b[h]]+=X*Pm[b[h]];
						fi=__min(X*Pm[b[h]],S[b[h]]);
						Pmeffect[h][r]=(X>0 ? fi/X : Pm[b[h]]);
						if(fi>0){
							Sales[b[h]] += fi;	M[owner[b[h]]] += (1 - tau) * fi / Pm[b[h]];
							S[b[h]] -= fi;		M[r] -= fi / Pm[b[h]];
						}
					}
				}
				for(k=1;k<=K;k++) if(active[k]>0) {
					S[k]-=__min(F,S[k]);
				}
				// Exit
				for(k=1;k<=K;k++) if(active[k])
				if( (Yp[owner[k]]<0 && rnd01<=theta) || rnd01 < theta*.01) {
					active[k]=0; g[k]=0; NS-=1; numfail++;
					familyshop[owner[k]]=0; owner[k]=0;
					Pm[k]=0; W[k]=0; trSales[k]=0;trInput[k]=0; Orders[k]=0; virtInput[k]=0;
					Sales[k]=0; Input[k]=0; S[k]=0;
						for(r=1;r<=m;r++) {
						if(employer[r]==k) { employer[r]=0; Weffect[r]=0;}
						for(h=0;h<=1;h++) if(store[h][r]==k) { store[h][r]=0; Pmeffect[h][r]=0;}
					}
				}
			
				// Update targets and prices
				for(k=1;k<=K;k++) if(active[k]) {
					trSales[k]= trSales[k] + alf*(Orders[k]-trSales[k]);
					trInput[k]=F+trSales[k]+ksi*(kappa*trSales[k]-S[k]);
					trInput[k]=__max(trInput[k],0);
					setprices();
				}
				// Random breakup of matches
				for (i=1;i<=m;i++) if(familyshop[i]==0 && rnd01 <= beta) {
					employer[i]=0; Weffect[i]=0;
					for(h=0;h<=1;h++) { store[h][i]=0; Pmeffect[h][i]=0;}
				}
				// Periodic (yearly) calculations during run
				if(t%(50*RptPer)==0) {
					lineup();
					Calc1();
					printf("%4d    %4.0f   %5.3f  %5.4f  %5.1f  %5.1f  %5.1f  %5.3f %4d %4d\n", NS, GDP, P/Pstar, w, Pb, MS/P, by, tau, numfail, t/50);
					AvGDP[t/50] +=(GDP-AvGDP[t/50])/run;
					AvFail[t/50] +=(numfail-AvFail[t/50])/run;
					MeanGDP[run] +=(GDP-MeanGDP[run])*50/t;
					MeanFail[run] += (numfail-MeanFail[run])*50/t;
					numfail=0;
				}
			}
			printf("End of run number %2d.  Adjustment speed %2.2f.\n\n", run, speed);
			stream=fopen("Results.fil", "a");
			fprintf(stream, "            %2.2f      %2d    %6.1f         %4.1f \n", speed, run, MeanGDP[run], MeanFail[run]);
			fclose(stream);
		}
		
	}
	printf("Press any key to end the program...\n\n");
	getch();
	return(0);
}

void Init(void) {
	r=0; d[0]=0; s[0]=0; Pm[0]=0.0; W[0]=0.0;
	for(i=0;i<=n;i++) {
		numprod[i]=0; numcons[i]=0;
		for(j=0;j<=n-2;j++) {
			produces[i][j]=0;
			consumes[i][j]=0;
		}
	}for(i=1;i<=n;i++) for(j=1;j<=n;j++) for(k = 1; k <= ( (i==j||i==1+(j%n)) ?  0 : 1); k++) {
		r++; s[r]=i; d[r]=j;
		numprod[i]++; produces[i][numprod[i]]=r;
		numcons[j]++; consumes[j][numcons[j]]=r;
	}
	for(i=0;i<=numruns;i++) {MeanGDP[i]=0; MeanFail[i]=0;}
	if(run==0) {
		stream=fopen("Results.fil", "w");
		fprintf(stream, "\Adjustment speed     run  Mean GDP   Mean # of failures.\n");
		fclose(stream);
	}
}

void InitRun(void) {
	printf("Num                   Real  \n");
	printf("Shops    GDP   P/P*   Wage    Pb      MS/P       by   tau    fail  year\n");
	numfail=0;
	srand(RAND+run-1);
	begin = clock();
	t=0;	NS=0; D = 0;
	Pstar=exp(-pistar);
	P=Pstar;
	MS=moneysupply;
	BS=bondsupply;
	tau=taxrate;
	Wagg=P*(1-tau)/(1+rho);
	by=debtratio;
	Pb=bondprice;
	for(r=0;r<=m;r++) {
		Yp[r]=0;
		E[r]=0;
		employer[r]=0;
		familyshop[r]=0; Weffect[r]=0;
		if(r>0) { M[r]=MS/m; B[r]=BS/m; }
		else { M[r]=0.0; B[r]=0.0; }
		for(h=0;h<=1;h++) {store[h][r]=0; Pmeffect[h][r]=0; }
	}
	for(k=0;k<=K;k++) {
		active[k] = 0; g[k]=0; owner[k]=0;
		Pm[k]=0; W[k]=0; trSales[k]=0; trInput[k]=0; Orders[k]=0; virtInput[k]=0;
		Sales[k]=0;Input[k]=0; S[k]=0; mu[k]=rho;
	}
	lineup();
}

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

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

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

int Calc1(void) {
	Wagg=0.0; Employment=0.0;
	for(k=1;k<=K;k++) if (active[k]) {Wagg+=W[k]*(Input[k]-1); Employment+= Input[k]-1;}
	Wagg=Wagg/Employment;
	w=Wagg/P; if(w>=9.999) w=9.999;
	return (0);
}

void setprices(void) {
	if(W[k] < 0.8 * Wagg && Yp[owner[k]] > 1.2*Wagg) W[k]=Wagg;
	//if(abs(trInput[k]/virtInput[k] - 1)>0.05)
		W[k] += __max( speed*0.25*alf*(trInput[k]/virtInput[k] - 1 ), -0.05 ) * W[k];
		W[k] = __min(W[k], 1000000.0);
	Pm[k] =  (1-tau)/((1+mu[k])*W[k]);
}

double Dem(double pm1, double pm2) {
	double X1;
	if(pm1==0) return(0.0);
	if(pm2==0) return(1.0);
	X1 = pow(pm1, eps)/(pow(pm1,eps) + pow(pm2,eps));
	if(X1>0) return(X1);
	else return(0.0);
}

int random(int upper) {
	return int(upper*rand()/(RAND_MAX+1.0));
}

int Research(double mkup) {
	int z1 = comrade();
	int z2 = 1+random(numcons[s[r]]);
	int z3 = consumes[s[r]][z2];
	if( Weffect[z1]<Wagg && Pmeffect[0][z3]<(1-tau)/((1+mkup)*Wagg) )
		return(1);
	else
		return(0);
}

double TaylorRule(double X, double Pb) {
	if(X==0) return(Pb);
	double X1 = pow( Pstar * capacity / X, bp );
	if(X1>0) return(__max(alf* X1/(v+pistar) + (1-alf)*Pb, 0.00001));
	else return(Pb);
}

double FiscalRule(double by, double PD, double tau) {
	double X = tau + 100*alf*(PD - exp(pistar) + 1) + (alf/100)*(by-bstar);
	X =	 __min( X, tauMax);
	return(X);
}
void SetupShop(int k, int r, double X, double mkup) {
	active[k]=1;
	g[k] = s[r];
	mu[k]=mkup;
	familyshop[r]=k;
	owner[k]=r;
	S[k]=0;
	employer[r]=k;
	Weffect[r]=1;
	trSales[k] = X;
	trInput[k]=F+trSales[k]+ksi*kappa*trSales[k];
	W[k] = Wagg;
	virtInput[k]=trInput[k];
	Pm[k] = (1-tau)/((1+mu[k])*W[k]);
}