/* Program: equifin.c
   Written by Pedro Amaral, Erwan Quintin
   Date started : April, 2004
   Purpose : computes equilibria for "Limited Enforcement, Financial Intermediation, and Economic Development " and calibrates parameters*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdflib.h"
#define NRANSI
#include "nrutil.h"
#define NMAX 100
#define GET_PSUM \
					for (j=0;j<ndim;j++) {\
					for (sum=0.0,i=0;i<mpts;i++) sum += p[i][j];\
					psum[j]=sum;}
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}




/*global variables*/
double alpha=0.3, nu=0.85, rho=-0.02,delta=0.88,A=30,theta=0.55;
int gridz=110; grida=100;
/*function find, finds zero of FOC in profit maximization problem*/
double find(double z, double k, double w, double nstar);

/*function out, calculate output given inputs*/
double outp2(double z,double k, double n);

/*function cumnorm, cumulative normal distribution*/

void cumnor(double* arg, double* result, double* ccum);

/*function amoeba from numerical recipes: downhill simplex*/

void amoeba(double p[3][2], double y[3], int ndim, double ftol,
	double (*funk)(double lam[2]), int *nfunk);

/*function amotry from numerical recipes: used by amoeba*/
double amotry(double p[3][2], double y[3], double psum[2], int ndim,
	double (*funk)(double lambda[2]), int ihi, double fac);


/* function ir : returns value of IR constraint given a and z */

double ir(double anow, /*assets*/
	   double dnow, /* loan*/
		double wnow, /* wage rate*/
	   double z,/*value of tfp for manager of current type*/  
	   double nstar, /*upper bound on labor for manager of type z*/
	   double eta);  /*enforcement parameter*/

/* function predict: returns intermediation ratio[1] and size distribution for given values of eta and lambdas*/  

void predict(double eta,  /*enforcement parameter*/
			   double lambda1, /*location of management type distribution*/
			   double lambda2, /*disperions of management type distribution*/
			   double ratio[1], /*intermediation ratio[1]*/
			   double cat1[1], /*fraction of small plants*/
	           double cat2[1], /* fraction of plants in second category*/
	           double cat3[1]); /*fraction of plants in thirs category*/ 

/*function loss: evaluates square distance between model and data dist,at eta such that yearly intermedation ratio is 2*/
double loss(double lambda[2]); 


/* main: minimize loss usind a downhill simplex method. (Amoeba function pp408-412 in numerical recipes for C). 
Downhill simplex method is robust, and requires no computation of derivatives. */
  
void main(void)
{
double lambda[2];
int ndim=2; /*dimension of the problem*/
double p[3][2]; /*function evaluated at vertices*/
double y[3];
double ftol=0.01;
int nfunk;
int i;

/*low alpha*/
/*lambda[0]=-2.194015;
lambda[1]=0.469479;*/
/*eta=0.5875;*/

/*high alpha*/
/*lambda[0]=-4.443012;
lambda[1]=2.312406;*/
/*eta=0.390625*/

/*low managerial share*/
/*lambda[0]=-3.916881;
lambda[1]=0.316743;
eta=0.531250;*/

/*high managerial share*/
/*lambda[0]=-3.817824;
lambda[1]=1.555104;*/
/*eta=0.418750;*/

/*parameters for benchmark*/

lambda[0]=-2.099260;
lambda[1]=0.657430;
/*eta=0.5;
*/

/* parameters for low beta */

 /* lambda[0]=-2.148205;
lambda[1]=0.611606;
eta=0.475*/




p[0][0]=-4;
p[0][1]=1;
p[1][0]=-2;
p[1][1]=1;
p[2][0]=-4;
p[2][1]=0.5;


/*printf("n= %f \n ",find(2,3,0.5,10));*/


/*i=0;
while(i<3)
{
lambda[0]=p[i][0];
lambda[1]=p[i][1];
y[i]=loss(lambda);

printf("i=%d y[i]=%f \n",i,y[i]);
i=i+1;
}


amoeba(p, y,ndim,ftol,loss,&nfunk);*/



loss(lambda);

}





double loss(double lambda[2])
{
double eta,etalow,etahigh,out,inceta;
double ca1[1],ca2[1],ca3[1];
double fin[1];
double lambda1,lambda2;
int it=0;


lambda1=lambda[0];
lambda2=lambda[1];
/*printf("lambda1=%f lambda2=%f \n",lambda1,lambda2);*/
fin[0]=0;
ca1[0]=0;
ca2[0]=0;
ca3[0]=0;

etalow=0;
etahigh=0.57;
inceta=etahigh/100;

/*eta=0.57;*/

printf("eta w cat1 avs totd totk toty tots ed\n");


it=1;
while(it<100)
{
eta=(double)it*inceta;
/*eta=0.57;*/
predict(eta,lambda1,lambda2,fin,ca1,ca2,ca3);



it=it+1;
}


/* old loss function below 
printf("out=%f \n",pow(.50-ca1[0],2)+pow(.33-ca2[0],2)+pow(.15-ca3[0],2));
return pow(.487-ca1[0],2)+pow(.327-ca2[0],2)+pow(.149-ca3[0],2);*/

/*new loss function: matches cat1 and avs and makes sure that ed is near 0*/

/*predict(eta,lambda1,lambda2,fin,ca1,ca2,ca3);*/


/*printf("out=%f \n",pow(50-100*ca1[0],2)+pow(50-ca3[0],2)+100000*pow(ca3[0],2));*/
if(pow(ca2[0],2)>0.01) out=100000000000;
else out=10*pow(50-100*ca1[0],2)+pow(50-ca3[0],2);

return out;
}

void predict(double eta, 
			 double lambda1, 
			 double lambda2, 
			 double ratio[1], 
			 double cat1[1], 
			 double cat2[1], 
			 double cat3[1]) 
{
/*global variables*/
double r;
double beta,betap;



double z[110];
double mu[110];
double kmax[110],nmax[110];
double cons,p,nk;
int grida=100;
double inca;
double V[110][100];
double d[110][100];
double save[110];
double save2[110];
double savew,savew2;
double n[110];
double y[110];
double k[110];
double deq[110];

double Ww,Wm[110];
double vnow,vnowp;

double wlow,whigh,wnow;
int i,j,it,itd,count;
double rat1;
int rat2,rat2p;
double weight;
double dlow,dhigh,dnow;
double slow,shigh,snow;
double anow;
double c1,c2,c3;
double know,nnow;
double ed;
double s1=9.5; /*top of first size category*/ 
double s2=49.5; /*top of second size category*/
double s3=249.5; /*top of third size category*/
double cat4,tot;
/*double *p1,*p2,*p3,*p4;*/
double up,down,res1,res2;

double totk,totd,toty,tots,totm;
 

/*printf("%f \n",alphap);*/




/*p1=&up;
p2=&down;
p3=&res1;
p4=&res2;*/

i=0;
while(i<110)
{
	
z[i]=pow(((double)i+1)/110,2);
mu[i]=1/(double)gridz;
i=i+1;
}

 totm=0;
i=0;
while(i<110)
{
	j=i-1;
	it=i+1;



	
if(i>0)
{
down=(z[i]+z[j])/2;
down=(log((down))-lambda1)/lambda2;
}



if(i<109)
{
up=(z[it]+z[i])/2;
up=(log(up)-lambda1)/lambda2;
}




if(i==0)
{
cumnor(&up,&res1,&res2);
mu[i]=res1;
}

if(i==109)
{
up=z[i]+(z[i]-z[j])/2;
up=(log(up)-lambda1)/lambda2;
/*printf("up=%f down=%f \n",up,down);*/

cumnor(&up,&res1,&res2);
mu[i]=res1;
cumnor(&down,&res1,&res2);
mu[i]=mu[i]-res1;
}

if(i>0 && i<109)
{
cumnor(&up,&res1,&res2);
mu[i]=res1;
cumnor(&down,&res1,&res2);
mu[i]=mu[i]-res1;
}



totm=totm+mu[i];
/*printf("i=%d z=%f mu=%f totm=%f\n",i,z[i],mu[i],totm);*/
i=i+1;
}


i=0;
while(i<110)
{
	mu[i]=mu[i]/totm;
	/*printf("i=%d z=%f mu=%f totm=%f\n",i,z[i],mu[i],totm);*/
	i=i+1;
}

totm=0;
wlow=0;
whigh=10;

r=pow(1.04,20)-1;
beta=1/(1+r);
beta=0.99/(1+r);
betap=beta*(1+beta);
p=r+delta;

/*whigh=.134;
wlow=.134;*/


/*begin iteratio[1]n on wages*/

it=1;
while (it<20)
  {

wnow=(whigh+wlow)/2; 

/*wnow=1;*/

inca=wnow/((double)grida-1);
nk=(1-alpha)*p/(alpha*wnow);
/*nk=pow(nk,1/(1-rho));*/

/*value function for workers*/

/*c1=1/(1+betap)*(wnow+wnow/(1+r));
c2=betap/(1+betap)*(1+r)*(wnow+wnow/(1+r));
c2=betap/(1+betap)*(1+r)*(wnow+wnow/(1+r));
Ww=log(c1)+betap*log(c2);
savew=wnow-c1;
savew2=beta/(1+beta)*c2;*/

/*printf("wnow=%f c1=%f c2=%f Ww=%f \n",wnow,c1,c2,Ww);*/

c1=1/(1+betap)*(wnow+wnow/(1+r));
c2=beta/(1+betap)*(wnow+wnow*(1+r));
c3=pow(beta,2.0)/(1+betap)*(wnow*(1+r)+wnow*pow(1+r,2.0));

if(c1>wnow) 
{
	c1=wnow;
	c2=1/(1+beta)*wnow;
	c3=beta/(1+beta)*wnow*(1+r);
}

Ww=log(c1)+beta*log(c2)+pow(beta,2.0)*log(c3);
savew=wnow-c1;
savew2=wnow+(wnow-c1)*(1+r)-c2;

/*printf("wnow=%f c1=%f c2=%f c3=%f Ww=%f \n",wnow,c1,c2,c3,log(c1));*/
/*printf("w=%f inca=%f c1=%f c2=%f Ww=%f \n",wnow,inca,c1,c2,Ww);*/
 
i=0;
while(i<110)
{

cons=alpha*A*z[i]/p*pow(theta*p/(alpha*wnow),theta);
kmax[i]=pow(cons,1/(1-alpha-theta));
nmax[i]=nk*kmax[i];

/*printf("z %f  n %f  k %f out %f  \n",z[i],nmax[i],kmax[i],outp2(z[i],kmax[i],nmax[i]));*/

/*algebra above was triple-checked*/

/*printf("w=%f z=%f cons=%f kmax=%f nmax=%f nk=%f\n",wnow,z[i],cons,kmax[i],nmax[i],nk);*/

j=0;
while(j<100)
{
	anow=j*inca;


dlow=0;
dhigh=kmax[i]-anow;

if(dhigh<0)
{
		/*printf("%d **** \n",j);*/
	V[i][j]=outp2(A*z[i],kmax[i],nmax[i])-kmax[i]*p-nmax[i]*wnow+anow*(1+r);
	d[i][j]=0;
	
		

}
else
{
count=0;
	itd=0;
while(itd<50)
{
	
dnow=(dlow+dhigh)/2;

/*if(j==109) printf("dnow=%f ir=%f \n",dnow,ir(anow,dnow,ap[i],eta));*/

if(ir(anow,dnow,wnow,z[i],nmax[i],eta)>0)
dlow=dnow;
else
{dhigh=dnow;
count=count+1;}

itd=itd+1;
}


dnow=(dhigh+dlow)/2;
if(count==itd) dnow=0; /*if ir never met, make dnow 0 instead of roughly 0*/
nnow=find(A*z[i],anow+dnow,wnow,nmax[i]);
d[i][j]=dnow;
V[i][j]=outp2(A*z[i],anow+dnow,nnow)-(anow+dnow)*p-nnow*wnow+anow*(1+r);
}

j=j+1;
}



/*lifetime value function for managers*/

slow=0;
shigh=wnow;
itd=0;
while(itd<20)
{
snow=(slow+shigh)/2;

rat1=snow/wnow*((double)grida-1);
rat2=(int)rat1;
weight=rat1-(double)rat2;
rat2p=rat2+1;

/*printf("%f %f %f %d %f\n",snow,wnow,rat1,rat2,weight);*/

vnow=weight*V[i][rat2p]+(1-weight)*V[i][rat2];
vnowp=(V[i][rat2p]-V[i][rat2])/inca;

 

if(1/(wnow-snow)>betap*vnowp/vnow)
shigh=snow;
else
slow=snow;

itd=itd+1;
}

save[i]=snow;
c1=wnow-snow;
c2=vnow*1/(1+beta);
c3=beta/(1+beta)*vnow*(1+r);
Wm[i]=log(c1)+beta*log(c2)+pow(beta,2.0)*log(c3);
save2[i]=vnow-c2;

/*printf("z=%f s=%f vnow=%f vnowp=%f c1=%f Wm[i]=%f Ww=%f cond=%f \n",z[i],save[i],vnow,vnowp,c1,Wm[i],Ww,1/c1-betap*vnowp/vnow);*/ 

if(Ww>Wm[i])
{
	n[i]=-1;
	y[i]=0;
	k[i]=0;
    deq[i]=0;
	save[i]=savew;
	save2[i]=savew2;}
else
{
dnow=weight*d[i][rat2p]+(1-weight)*d[i][rat2];
know=save[i]+dnow;
nnow=find(A*z[i],know,wnow,nmax[i]);;
n[i]=nnow;
y[i]=outp2(A*z[i],know,nnow);
k[i]=know;
deq[i]=dnow;
/*if(i==109) printf("vnow=%f wm=%f dnow=%f know=%f nnow=%f n[i]=%f \n", vnow,Wm[i],dnow,know,nnow,n[i]);*/
}

/*if(eta>0.4)
printf("w=%f z=%f n=%f k=%f kmax=%f mu=%f \n",wnow,z[i],n[i],know,kmax[i],mu[i]);*/
/*if (it==19) printf("i=%d w[i]=%f save2[i]=%f \n",i,Wm[i],save2[i]);*/
i=i+1;
}

/*calculate aggregate excess labor demand*/

ed=0;
cat1[0]=0;
cat2[0]=0;
cat3[0]=0;
cat4=0;
totk=0;
toty=0;
totd=0;
tots=0;
totm=0;
i=0;
while(i<110)
{
	ed=ed+mu[i]*(n[i]-1);
	totk=totk+mu[i]*k[i];
	toty=toty+mu[i]*y[i];
	totd=totd+mu[i]*deq[i];
	tots=tots+mu[i]*(save[i]+save2[i]);
	
	if(n[i]>0)
	{
			totm=totm+mu[i];
		if (n[i]<s1) cat1[0]=cat1[0]+mu[i];
		else if(n[i]<s2) cat2[0]=cat2[0]+mu[i];
			else if (n[i]<s3) cat3[0]=cat3[0]+mu[i];
			else cat4=cat4+mu[i];}

	i=i+1;

}

	tot=cat1[0]+cat2[0]+cat3[0]+cat4;
	cat1[0]=cat1[0]/tot;
	cat2[0]=cat2[0]/tot;
	cat2[0]=ed;
	cat3[0]=cat3[0]/tot;
	cat3[0]=(2-totm)/totm;
	cat4=cat4/tot;

if(ed>0)
wlow=wnow;
else
whigh=wnow;

/*printf("wnow=%f ed=%f eta=%f nmax=%f\n",wnow,ed,eta,nmax[49]);*/

/*i=0;
while(i<110)
{
printf("wnow=%f z=%f save=%f Wm[i]=%f Ww=%f n[i]=%f ed=%f \n",wnow,z[i],save[i],Wm[i],Ww,n[i],ed);
i=i+1;}
*/ 

it=it+1;
}


/*printf("wnow=%f z[49]=%f \n",wnow,z[49]);*/
/*j=0;
while(j<grida)
{
	i=j-1;
printf("%f %f %f %f %f %f %f %f %f %f %f\n",(double)j*inca,d[49][j],(double)j*inca+d[49][j],kmax[49],cat1[1],cat2[1],cat3[1],cat4,totd/toty,totk/toty,alpha/(1+r));
j=j+1;
}
*/
ratio[0]=totd/(toty+(tots-totk)*r);

/*printf("*****eta=%f lambda1=%f lambda2=%f ratio[1]=%f cat1[1]=%f cat2[1]=%f cat3[1]=%f %f \n",eta,lambda1,lambda2,ratio[0],cat1[0],cat2[0],cat3[0],cat4);*/

printf("%.5f %.6f %.5f %.5f %f %f %f %f %f\n",eta,wnow,cat1[0],cat3[0],totd,totk,toty,tots,ed);



/*printf("i z[i] mu[i] n[i] \n");
i=0;
while(i<110)
{
printf("%d %.5f %.6f %.5f \n",i,n[i],z[i],save[i]);
i=i+1;
}
*/

/*printf("eta=%.3f cs=%f ratio=%f cat1=%.3f  cat3=%.3f housing=%.4f \n",eta,(tots*(delta+pow(1.04,20)-1))/(toty+(tots-totk)*(delta+pow(1.04,20)-1)),ratio[0],cat1[0],cat3[0],(tots-totk)/tots);*/
}





double ir(double anow, 
	   double dnow, 
	   double wnow,
	   double z,
	   double nstar,
	   double eta)   
{
double r,p;
double know,nnow;
double out;


r=pow(1.04,20)-1;
p=r+delta;
know=anow+dnow;
nnow=find(A*z,know,wnow,nstar);

out=eta*(outp2(A*z,know,nnow)-know*p-nnow*wnow+know*(1+r))-dnow*(1+r);

return out;
}




void amoeba(double p[3][2], double y[3], int ndim, double ftol,
	double (*funk)(double []), int *nfunk)
{
	double amotry(double p[3][2], double y[3], double psum[2], int ndim,
		double (*funk)(double lamda[2]), int ihi, double fac);
	int i,ihi,ilo,inhi,j,mpts=ndim+1;
	double rtol,sum,swap,ysave,ytry,*psum;

		printf("nfunk=%d \n",nfunk);

	psum=vector(0,ndim-1);
	*nfunk=0;
	GET_PSUM
	for (;;) {
		ilo=0;
		ihi = y[0]>y[1] ? (inhi=1,0) : (inhi=0,1);
		for (i=0;i<mpts;i++) {
			if (y[i] <= y[ilo]) ilo=i;
			if (y[i] > y[ihi]) {
				inhi=ihi;
				ihi=i;
			} else if (y[i] > y[inhi] && i != ihi) inhi=i;
		}
		rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
		if (rtol < ftol) {
			SWAP(y[1],y[ilo])
			for (i=0;i<ndim;i++) SWAP(p[1][i],p[ilo][i])
			break;
		}
		if (*nfunk >= NMAX) nrerror("NMAX exceeded");
		*nfunk += 2;
		ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0);
		if (ytry <= y[ilo])
			ytry=amotry(p,y,psum,ndim,funk,ihi,2.0);
		else if (ytry >= y[inhi]) {
			ysave=y[ihi];
			ytry=amotry(p,y,psum,ndim,funk,ihi,0.5);
			if (ytry >= ysave) {
				for (i=0;i<mpts;i++) {
					if (i != ilo) {
						for (j=0;j<ndim;j++)
							p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
						y[i]=(*funk)(psum);
					}
				}
				*nfunk += ndim;
				GET_PSUM
			}
		} else --(*nfunk);
	}
	/*free_vector(psum,1,ndim);*/
}


double amotry(double p[3][2], double y[3], double psum[2], int ndim,
	double (*funk)(double lambda[2]), int ihi, double fac)
{
	int j;
	double fac1,fac2,ytry,*ptry;

	ptry=vector(0,ndim-1);
	fac1=(1.0-fac)/ndim;
	fac2=fac1-fac;
	for (j=0;j<ndim;j++) 
	{ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
	printf("j=%d ihi=%d psum=%f fac1=%f p=%f fac2=%f \n",j,ihi,psum[j],fac1,p[ihi][j],fac2);}

	ytry=(*funk)(ptry);
	if (ytry < y[ihi]) {
		y[ihi]=ytry;
		for (j=0;j<ndim;j++) {
			psum[j] += ptry[j]-p[ihi][j];
			p[ihi][j]=ptry[j];
		}
	}
	/*free_vector(ptry,1,ndim);*/
	return ytry;
}


double find(double z, double k, double w, double nstar)
{
	int it=0;
    double nlow,nhigh,nnow,cons;
	
nnow=theta*z*pow(k,alpha)/w;
nnow=pow(nnow,1/(1-theta));

return nnow;

/*this function's algebra was triple checked*/
}





double outp2(double z,double k, double n)
{
	double res;
	/*res=alpha/pow(k,-rho)+(1-alpha)/pow(n,-rho);
	return z*pow(1/res,-nu/rho);*/
	res=z*pow(k,alpha)*pow(n,theta);
	return res;
	/*this function's algebra was triple checked*/
}