/*===============================================================
   This code solves and simulates a basic Pissarides-style labor market search
   model.
   David Andolfatto (January 2006).
================================================================*/

library nlsys; /* Call the nonlinear equations solver */

/*=============
   PARAMETERS 
==============*/

/* Discount Factor */
beta = 0.9959; 
/* Separation Rate */
sep = 0.034;
/* Unemployment Rate */
ustar = 0.07;
/* Job-Finding Rate */
pstar = sep*(1-ustar)/ustar;
/* Matching Elasticity */
alpha = 0.50;
/* Scale Parameter on Matching Technology */
mtech = pstar;
/* Bargaining Parameter */
zeta = alpha;
/* Leisure Parameter */
zed = 0.40;
/* Vacancy Cost */
kappa = ( beta*pstar*(1-zeta) )*(1 - beta*(1-sep-zeta*pstar) )^(-1)*(1-zed);

/* Stochastic Process for Productivity */
rho = 0.97;
sigma = 0.007;

/*===================
   GRID SPECIFICATION 
====================*/

lbound = -2.5*sigma*(1-rho^2)^(-1/2);
ubound = -lbound;
ngrid = 11;
step = (ubound-lbound)/(ngrid-1);
shock = seqa(lbound,step,ngrid);

/*==================================================
   CONSISTENCY CHECK: 
   Let's make sure that the match probabilities are well-defined
   at each point in the grid.
==================================================*/

	i = 1;
	do while i le ngrid;

		__output=0;
		x0 = 1;
		{x,f,g,h} = nlsys(&steady,x0);

		pp(x)~qq(x);

	i = i+1;
	endo;
	print; "Press ENTER to Continue."; cons; 

/*===================================================
  SOLVE FOR POLICY FUNCTIONS
===================================================*/

sss = (1 - beta*(1-sep-zeta*pstar) )^(-1)*(1 - zed); /* steady-state surplus */

surplus = ones(ngrid,1)*sss;
xsurplus = surplus;
theta = ones(ngrid,1);
xtheta = theta;

tolerance = 0.0000001;
frac = 0.1;

        /*===========================================*/
	/* Compute theta(x), conditional on value function S(x) */
        /*===========================================*/

	converge1 = 0;
	loop = 1;
	do while converge1 == 0;
	"Loop: " loop; print;

	        /* Compute value function S(x), conditional on theta(x) */

		converge2 = 0;
		sloop = 1;
		do while converge2 == 0;
		/* "sloop: " sloop;  */

		i = 1;
		do while i le ngrid;

   			ESURP = (0.5)*INTERP(surplus,rho*shock[i]+sigma) + (0.5)*INTERP(surplus,rho*shock[i]-sigma); 
			xsurplus[i] = exp(shock[i]) - zed + beta*(1 - sep - zeta*pp(theta[i]))*ESURP;
	
		i = i+1;
		endo;

		/* Check for converge of value function */
		
		/* "Value Function Difference: " abs( sumc( (ln(xsurplus)-ln(surplus))/ngrid ) ); */
		if abs( sumc( (ln(xsurplus)-ln(surplus))/ngrid ) ) le tolerance; converge2 = 1; endif;

		/* Update value function */
		surplus = xsurplus;

		sloop = sloop+1;
		endo;

	i = 1;
	do while i le ngrid;

		__output = 0;
		x0 = theta[i];
		{x,f,g,h} = nlsys(&foc,x0);

		xtheta[i] = x;
 
	i = i+1;
	endo;


	/* Check for converge of policy function */
		
	"Policy Function Difference: " abs( sumc( (ln(xtheta)-ln(theta))/ngrid ) );
	if abs( sumc( (ln(xtheta)-ln(theta))/ngrid ) ) le tolerance; converge1 = 1; endif;

	/* Update policy function */
	theta = (1-frac)*xtheta+frac*theta;

	loop = loop+1;
	endo;

/*=================
  SIMULATIONS
=================*/
 
slength = 10000;
sshock = zeros(slength+1,1);
vac = zeros(slength,1);
ump = sshock;
prob = vac;

	print;
	"Simulating the model..."; 
	print;

	/* Generate a sequence of shocks */
	j = 1;
	do while j le slength;

		if rndu(1,1) le 0.5;
		innovation = -sigma;
		else;
		innovation = sigma;
		endif;
	
		sshock[j+1] = rho*sshock[j] + innovation;

	j = j+1;
	endo;

	/* Use shock sequence to generate time-paths for variables of interest */

	ump[1] = ustar;
	j = 1;
	do while j le slength;

		stheta = INTERP(theta,sshock[j]);
		prob[j] = pp(stheta);
		ump[j+1] = ump[j] + sep*(1-ump[j]) - prob[j]*ump[j];
		vac[j] = stheta*ump[j];

	j = j+1;
	endo;

/* You may at times wish to plot your simulated data. Suppose, for example, we wish have a file
    containing vacancies and unemployment (the first 100 simulations). You can write this data
    to a file as follows:

    output on;
    output file = path\filename reset; (e.g., output file = c:\gauss\code\output.txt reset; )
    vac[1:100]~ump[1:100];
    output off;
==============================================================================*/


	/* Compute Statistics */

	stdvac = 100*stdc( ln(vac) );
	stdump = 100*stdc( ump[1:slength] );	
	corv = corrx(vac~ump[1:slength]); 

"Std Dev (%) Vacancies:         " stdvac;
"Std Dev (%) Unemployment: " stdump; 
"Correlation(u,v):                    " corv[1,2];
print;

"Program Terminated."; print;

end;

/*========================================================================*/
proc foc(x);
local fn1, ESURP;

ESURP = (0.5)*INTERP(surplus,rho*shock[i]+sigma) + (0.5)*INTERP(surplus,rho*shock[i]-sigma); 

	fn1 = beta*(1-zeta)*qq(x[1])*ESURP - kappa;

retp(fn1);
endp;

/*========================================================================*/
proc steady(x);
local fn1;

	fn1 = kappa*(1-beta*(1-sep-zeta*pp(x[1]) )) - (exp(shock[i]) - zed)*beta*qq(x[1])*(1-zeta);

retp(fn1);
endp;


/*========================================================================*/
proc pp(x);
retp( mtech*x^alpha );
endp;


/*========================================================================*/
proc qq(x);
retp( mtech*x^(alpha-1) );
endp;

/*========================================================================*/
proc INTERP(rule,xx);
local i, w1;

	i = maxindc( shock .ge xx );
	i = maxc(2|i);
	w1 = (1/step)*( xx - shock[i-1] );

	retp( (1-w1)*rule[i-1] + w1*rule[i] );
endp;