@-GAUSS CODE for A NON-STOCHASTIC GROWTH MODEL-@
@-DEFINE PARAMETERS, DEFINE SOME USEFUL FUNCTIONS, COMPUTE STEADY STATE-@

new;

library nlsysmt, pgraph;

#include nlsysmt.sdf

struct nlControl nlc;
struct nlOut nlo;
nlc = nlControlCreate;
nlc.fvtol = 0.0000000001;
nlc.output = 0;

graphset;
_pdate=0; _pnum=2; _plwidth = 6; 
_plclr = 6; _pltype={6,3,2}; _pgrid=1|0;

@ Define parameter values @

beta  = 0.99;  
delta = 0.02; 
theta = 0.36; 
gama  = 1.5;  

nvar  = 2;

@--  Order of variables
x[1] = capital   
x[2] = cons     --@

@ Some useful functions @

fn rrate(x) = theta*x[1]^(theta-1);
fn prodn(x) = x[1]^theta;
fn muc(x) = x[2]^(-gama);

@--Define the left and right hand sides of the First Order Conditions--@

fn lhs1(x)=x[1];     
fn rhs1(x)=(1-delta)*x[1]+prodn(x)-x[2];

fn lhs2(x)=beta*muc(x)*(1-delta+rrate(x));  
fn rhs2(x)=muc(x); 

@-Define a procedure to check whether FOCs are satisfied-@

proc sse(x);

local e;
e = zeros(nvar,1);

e[1] = lhs1(x)-rhs1(x);
e[2] = lhs2(x)-rhs2(x);

retp(e); endp;


@-Make an initial guess for the steady state-@

x0=zeros(nvar,1);

x0[1] = 58;           
x0[2] = 4;        

@-Compute Steady State-@

{nlc, nlo}=nlsys(nlc, &sse,x0);

xbar=nlo.xp;

print "steady state capital, consumption";
xbar';


@-COMPUTE DECISION RULES-@

psi=(gradp(&lhs1,xbar)|gradp(&lhs2,xbar));
phi=(gradp(&rhs1,xbar)|gradp(&rhs2,xbar));


A = inv(psi)*phi;
{ L,Q } = eigv(A);

IQ = inv(Q);

if (sumc(abs(L).>1)/=1);
 print "either too few or too many unstable roots";
 end;
endif;

if abs(L[1])>1;
 ck = -IQ[1,1]/IQ[1,2];
elseif abs(L[2])>1;
 ck = -IQ[2,1]/IQ[2,2];
endif;
 
con_rule = ck;
print "coefficient on capital in consumption rule " con_rule;

cap_lom = A[1,1] + A[1,2]*con_rule;
print "coefficient on capital in capital lom " cap_lom;


@-Compute nonstochastic dynamics to steady state-@

horizon = 200;

cap_seq = zeros(horizon+1,1);
con_seq = zeros(horizon,1);

cap_seq[1] = 0.1*xbar[1];

i=1; do while i<=horizon;

	cap_seq[i+1] = cap_lom* (cap_seq[i]-xbar[1]) + xbar[1];
	con_seq[i] = con_rule* (cap_seq[i]-xbar[1]) + xbar[2];

i=i+1; endo;

@
title("path for capital");
xy(0,cap_seq);

title("path for consumption");
xy(0,con_seq);
@

@ Define eigenvectors in a graphable form @

t1=-Q[2,1]/Q[1,1]|0|Q[2,1]/Q[1,1];
t2=-Q[2,2]/Q[1,2]|0|Q[2,2]/Q[1,2];

@ Compute the dynamics for capital and consumption for an arbitrary initial starting point @

kseq=zeros(100,1);
cseq=zeros(100,1);

kseq[1]=0.1;
cseq[1]=0.01;

i=2; do while i<=100;

t=A*(kseq[i-1]|cseq[i-1]);

kseq[i]=t[1];
cseq[i]=t[2];

i=i+1;
endo;

@ Produce the picture I showed in class @

window(1,1,0);
setwind(1);

title("eigenvalues and eigenvectors");
xlabel("khat");
ylabel("chat");
xtics(-1,1,0.5,1);
ytics(-0.04,0.04,0.02,1);
xy((-1|0|1),t1~t2);

xy(kseq,cseq);

endwind;