/* LUCAS   ------- Simplified Lucas model
                   gama = 0     no external effect of human capital*/

library projec2;

/* Parameters of model */
alpha   = .6;               /* capital share */
rho     = .1;               /* discount rate */
theta   = .9;               /* risk aversion */
phi     = .15;              /* effectiveness of investment in human capital */
delta   = .1;               /* depreciation rate */

/* calibration on steady state */
uss     = 1-(phi-rho)/(theta*phi);
ua      = uss^(1-alpha);
aa      = (phi+delta)/ua/alpha;
css     = (delta+phi)/alpha-(phi-rho)/theta-delta;

nst     = 1;    /* number of state variables */
ncon    = 2;    /* number of control variables */
nfc     = 10;    /* degree of approximation */
mcons   = nst~nfc~ncon;
case    = 1;    /* case = 1: normal solution;
                          2: theta=alpha */
x       = { .5, 1.5 }; /* interval of state variables */
/*a0      = zeros(nst+1,ncon); /* initial linear guess */*/
a0      = zeros(4,1);

/* procedure for residual function */
proc _FRES(x,a);
local par,R1,R2,h0,h1,h2,h3,h4,h5,h6,h7,h8,h9;
par = reshape(a,2,nfc)';
R1  = 0; R2 = 0;
h1  = __APROX(x,par[.,1]);                      /* c(k) */
h2  = __APROX(x,par[.,2]);                      /* u(k) */
h3  = __DAPROX(x,par[.,1],1);                  /* c'(k) */
h4  = __DAPROX(x,par[.,2],1);                  /* u'(k) */
h5  = aa*__GPOW(h2/x,1-alpha);                  /* AA*(h2/x)^(1-alpha) */
if case == 1;
    h6 = (h5-h1-delta-phi*(1-h2))*x*h3;              /* 1.term of 1.equation */
    h7 = ((alpha*h5-rho-delta)/theta-h5+h1+delta)*h1;/* 2.term of 1.equation */
    h8 = (h5-h1-delta-phi*(1-h2))*x*h4;              /* 1.term of 2.equation */
    h9 = ((delta+phi)/alpha-h1-delta-phi*(1-h2))*h2; /* 2.term of 2.equation */
else;
    h6 = (h5-h1-delta-phi*(1-h2))*x*h3;              /* 1.term of 1.equation */
    h7 = ((-rho-delta)/theta+h1+delta)*h1;           /* 2.term of 1.equation */
    h8 = (h5-h1-delta-phi*(1-h2))*x*h4;              /* 1.term of 2.equation */
    h9 = ((delta+phi)/alpha-h1-delta-phi*(1-h2))*h2; /* 2.term of 2.equation */
endif;
R1  = h6-h7;
R2  = h8-h9;
retp(R1|R2);
endp;

/* initial value of parameters */
a0    = a0 + .5;

PROJSET;
_prmeth = 2;  /* it is necessary for convergence to start with STEEP */
_prsave = 1;
_prldfn = "a_2";
_prsvfn = "a_10";
_prinit = 1;
{a,ret} = PROJEC(&_FRES,mcons,x,a0);