/* ================================================================= 
    This codes solves for the optimal trajectory for the capital stock in a 
     deterministic optimal growth model. 
 
     Preferences: additive log separable (with constant discount factor).
     Technology: Cobb-Douglas
     Constant depreciation rate.

     David Andolfatto
     Simon Fraser University
     February 2004
  =================================================================*/

library nlsys; 

/* 	PARAMETERS */

alpha = 0.35; /* capital's share of income */
beta   = 0.96; /* discount factor */
delta  = 0.10; /* depreciation rate */

/* Solve for the steady-state capital stock and choose the tfp parameter to normalize
    steady-state ouput to unity. */

kstar = alpha/( 1/beta -1 + delta); /* Here, I have used the fact that f' = alpha*y/k, where y = 1 */

tfp = 1/kstar^alpha;

length = 100;

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

guess = ones(length,1)*(1/2)*kstar; /* Initial Guess */
soln    = guess;
soln[1] = 0.25*kstar; /* Initial Condition */
guess[1] = soln[1];
soln[length] = kstar;
guess[length] = kstar;

converge = 0;

	iteration = 1;
	do while converge == 0;

		cls;
		"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;	

                diff = abs( ln(soln) - ln(guess) );

		print;
		"Sum of Differences: " sumc(diff);
		print;

		if diff le 0.0001;
		converge = 1;
		endif; 

		guess = soln;  /* Update your guess */

	iteration = iteration + 1;
	endo;

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

/* Generate optimal GDP path */

gdp = zeros(length,1);

t = 1;
do while t le length;

gdp[t] = tfp*soln[t]^alpha;

t = t+1;
endo;

end;



/*===============================================================================*/
proc FOC(x);
	local fn1, c1, c2;

	c1 = tfp*guess[i]^alpha + (1-delta)*guess[i] - x[1];
        c2 = tfp*x[1]^alpha + (1-delta)*x[1] - guess[i+2];

	fn1 = c2 - ( alpha*tfp*x[1]^(alpha-1) + 1 - delta )*beta*c1;

	retp( fn1 );
endp;