/* INPUT EX4 -----  stochastic growth model from Judd,
                    two states, one approximated control */

library projec2;

/* Parameters of model */
sig     = .001;         /* dispersion */
beta    = 1/(1+.05);    /* subj. discount factor */
gam     = -.9;
alpha   = .25;          /* production f(k)= AA*k^alpha */
rho     = .8;           /* productivity shock */

/* calibration on steady state */
kss     = 1;                        /* normalization of states */
thss    = 1;
AA      = 1/(alpha*beta*thss);      /* productivity */
css     = thss*AA*kss^alpha-kss;

nst     = 2;                    /* number of state variables */
ncon    = 1;                    /* number of control variables */
nfc     = {2 2};                /* degrees of approximation */
nnn     = nst~nfc~ncon;         /* basic constants of model */
x       = {.3 .3, 1.7 1.7};     /* intervals of state variables */
par     = MULT(1,nfc);
a0      = zeros(3,ncon);        /* matrix of initial values of parameters */
                                /* linear approximation */

proc INTEG(z,x,a);
local k1,th1,h1,mpk1,y,k,th,h;
    k=x[1]; th=x[2];
    h       = __APROX(x,a);                              /* h(k,th) */
    k1      = th*aa*k^alpha-h;                          /* k(t+1) */
    th1     = exp(sig*sqrt(2)*z)*th^rho;                /* th(t+1) */
    h1      = __APROX(k1~th1,a);                         /* h(k1,th1) */
    mpk1    = alpha*AA*__GPOW(k1,alpha-1);               /* mpk(k1) */
    y       = __GPOW(h1,gam)*th1*mpk1/sqrt(pi);
    retp(y);
endp;

/* residual function for stochastic growth */
proc _FRES(x,a);
local inte,r,k,h;
    k=x[1];
    h       = __APROX(x,a);              /* h(k,th) */
    inte    = GH_QUAD(&INTEG,x,a,6);    /* approx. of expect. value */
    r       = h-(beta*inte)^(1/gam);    /* residuum */
    retp(r);
endp;

/* initial value of parameters */
a0[2] = (css-css/2)/thss;
a0[3] = css/(2*kss);

PROJSET;
_prsave = 1;
_prsvfn = "a_2-2";
{a,ret} = PROJEC(&_FRES,nnn,x,a0);