/* ================================================================= 
    Here is some code that answers question 1b from Assignment 3

     David Andolfatto
     Simon Fraser University
     March 2001 
  =================================================================*/

library nlsys; 

/* 	PARAMETER VALUES */

alpha = 0.35;
beta   = 0.96;

kstar = (alpha*beta)^(1/(1-alpha));  /* solve for the steady state capital stock */

/* Now let's try to figure out in which period the steady state is almost reached */

ks = zeros(100,1);

	ks[1] = 0.01*kstar;

	j = 1;
	do while j le 100;

		if abs( ln(ks[j]) - ln(kstar) ) le 0.01;
		goto SSDATE;
		endif;

		ks[j+1] = alpha*beta*ks[j]^alpha;

	j = j+1;
	endo;

SSDATE:
length = j;   
truesoln = ks[1:length];

" "; 
"We are 1% away from the steady state in: " length " periods."; 
" ";
"The true solution is given by: ";
truesoln; " "; 
"Press ENTER to continue."; cons; /* This is an easy way to pause the output
							    to your screen. */ 

/* Now, let's solve for the optimal path. */

guess = ones(length,1)*(1/2)*kstar;
soln    = guess;
soln[1] = 0.01*kstar;
soln[length] = kstar;
guess[1] = soln[1]; 
converge = 0;

	iteration = 1;
	do while converge == 0;

		"Iteration: " iteration; print;

		__output=0;
		i = 1;
		do while i le length-2; 

			x0 = soln[i+1];
			{x,f,g,h} = NLSYS(&FOC,x0);

			soln[i+1] = x; 

	        i = i+1;
		endo;	

		if abs( ln(soln) - ln(guess) ) le 0.0001;
		converge = 1;
		endif; 

		guess = soln;  /* Update your guess */

	iteration = iteration + 1;
	endo;

print; 
"Convergence achieved in " iteration " iterations."; 
print;
print;

format /rd 20,4;
"        True Solution    Approximate Solution "; 
truesoln~soln;
print;


end;



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

	fn1 = x[1]^alpha-guess[i+2] - 
		alpha*beta*x[1]^(alpha-1)*(guess[i]^alpha - x[1]);

	retp( fn1 );
endp;

/*================================================================================*/
/* 	Interpolation Procedure: Here, I want a function that takes as input some arbitrary shock 
	value (not necessarily on the grid) and return the value implied by the decision rule (an 
	NZx1 vector of points).  */

proc INTERP(rule,zz);
	local i, w1;

	i = maxindc( zshock .ge zz );

	w1 = (1/step)*( zz - zshock[i-1] );

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

endp;

/*================================================================================*/