@Sample GAUSS Program Illustrating how to find optimal decision rules for
     the Basic Model.  The following functional forms are used:
     U(c,l) = log c + a log l
     F(k,h) = k^theta * h^(1-theta)
     A(z) = 1 - rho + rho z         @
new;
format 9,5;

@ Step 1.  Setting up the model @

@ Parameter Values :@
     theta=.4;
     delta=.025;
     rho=.95;
     beta=.99;
     a=2;

ns=3;  @Number of state variables (1,z,k)@
nd=2;  @Number of decision variables (h,i)@

@ Procedure that defines the return function, r(z,k,h,i), where
     y[1,1] = z
     y[2,1] = k
     y[3,1] = h
     y[4,1] = i  @
proc retf(y);
     local r;
     r=ln(y[1,1]*y[2,1]^theta*y[3,1]^(1-theta)-y[4,1])
                  +a*ln(1-y[3,1]);
     retp(r);
endp;

@ Compute steady state values of x which will be used to form
     a quadratic approximation of the return function. @

ys=zeros(ns+nd-1,1);
disrat=1/beta-1;
ys[1,1]= 1;
ys[3,1]=(1-theta)*(disrat+delta)/((a+1)*(disrat+delta)-theta
                  *(disrat+(a+1)*delta));
ys[2,1]=(((disrat+delta)/theta)^(1/(theta-1)))*ys[3,1];
ys[4,1]=ys[2,1]*delta;


@ Define a matrix c containing the laws of motion for all state variables:
   1'= 1
   z'= 1 - rho + rho * z
   k'= (1-Delta)*k + i     @

c=zeros(ns,2*ns+nd);
c[1,1]=1;
c[2,1]=1-rho;
c[2,2]=rho;
c[3,3]=1-delta;
c[3,5]=1;

@ Step 2.  Compute quadratic approximation of return function. @

q=quad(ys,0.000001,&retf);

@ Step 3.  Solve dynamic program.  @

test=10;
n=ns+nd;
v=ones(ns,ns)*(-1000); format /rd 8,5;
iter=0;
do until test lt 1E-8;
     tv=q~zeros(n,ns)|zeros(ns,n)~(v*beta);
     nv=rows(tv);

     @Reduce out laws of motion for state variables.@

     i=1;
     do until i>ns;
          tv=reduce(tv,c[ns-i+1,1:nv-i]);
          i=i+1;
     endo;

     @Reduce out first order conditions.@
     dsave=tv[ns+1:n,1:n];
     i=1;
     do until i>nd;
          d=-tv[n-i+1,1:n-i]./tv[n-i+1,n-i+1];
          if tv[n-i+1,n-i+1]>0; format 3,0; "Second order condition fails!";;
             "  Decision variable number ";; i;;"at iteration ";; iter+1;
             format 9,5; "Number checked is equal to ";; tv[n-i+1,n-i+1];
          endif;
          tv=reduce(tv,d);
          i=i+1;
     endo;
     test=abs(tv-v); test[1,1]=0;
     v=tv;
     iter=iter+1;
endo; format 4,0;
"Dynamic program required ";; iter;; "iterations."; ?;

@ Step 4.  Compute decision rules. @

decis=-dsave[.,1:ns]/dsave[.,ns+1:n];

@ Step 5.  Compute steady states from decision rules and compare
          with original steady states. @

ssst=1|1|(decis[2,1]+decis[2,2])/(delta-decis[2,3]);
ss=decis*ssst;
ss=ssst|ss;
format 10,5;
"Original Steady State: ";;ys;
"Steady State from Approximation: ";;ss;
end;