/* ================================================================= 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; /*================================================================================*/