/* determgrowth.prg: uses a projection algorithm to solve the deterministic growth model 
				approximation based on a weighted least squares objective function */

new;
library optmum,pgraph;
graphset;
_pdate = "";
_psymsiz =  3;
_paxht   = .2;
_pnumht  = .2;
_ptitlht = .3; 
_paxes   = { 1, 1 }; 

@ establish parameters @

alp = 0.33;
bet = 0.96;
del = 0.1;
phi = 2;
p = alp|bet|del|phi;

ss = steady(p);
kss = ss[4];

/* establish range for state */

wdth = 0.2;
kl = kss-wdth*kss;
kh = kss+wdth*kss;
ksteps = 50;
step = (kh-kl)/ksteps;
krange = seqa(kl,step,ksteps+1);
knorm = (krange-kss)/(wdth*kss);
krange~knorm;

/* starting values for gam */

gam0 = 1.2|0.1|0|0|0|0|0|0|0|0|0;

/* establish weighting function g(k) */

gee = ones(rows(krange),1);
gee = 1-(abs(krange-kss)/kss)^.1;

{ gamopt,fopt,gopt,retcode } = optmum(&feval,gam0);

ceeopt = ceeofk(knorm,gamopt);

policy = ceeopt~krange;

/* calculate slope of policy function */

dcdk = gradp(&ck,knorm);
dcdk = diag(dcdk)/(wdth*kss);

/* graph results */

begwind;
window(2,1,0);

hor = krange[26];
ver = ceeopt[26];
_pmsgctl = { 3.53 1.14 .2 0 1 1 0 };
_pmsgstr = "(C[*], K[*])";
_psym = { 3.53 1.163352 7 7 1 1 0};
xlabel("K");
ylabel("C");
xtics(2.8,4.2,.16,0);


xy(krange,ceeopt);
nextwind;

_pmsgctl = { 3.54 .14 .2 0 1 1 0 };
_pmsgstr = "(C[*], K[*])";
_psym = { 3.53 .139326 7 7 1 1 0};
ylabel("dC/dK");

xy(krange,dcdk);

endwind;

/* collect procedures */

proc steady(p);
    /* calculates the steady state of the model */
    local alp,bet,del,kss,yss,iss,css;
    alp = p[1];
    bet = p[2];
    del = p[3];

    kss = (alp/(1/bet-1+del))^(1/(1-alp));
    yss = kss^alp;
    iss = del*kss;
    css = yss - iss;

    retp(yss|css|iss|kss);
endp;

proc ceeofk(k,gam);
    /* evaluates c(k) */
    local  ntees,nks,tees,j,cee;

    ntees = rows(gam);
    nks = rows(k);   
    tees = zeros(nks,ntees);
    tees[.,1] = ones(nks,1);
    tees[.,2] = k;
    j=3; do while j<=ntees;
        tees[.,j] = 2*k.*tees[.,j-1]-tees[.,j-2];
    j=j+1; endo;    
    cee = tees*gam;

    retp(cee);
endp;

/* function evaluating proc */

proc feval(gam);

    local nrows,cee,cp,f,kp,kpnorm,j,mtrk;

nrows = rows(knorm);
cee = zeros(nrows,1);
f = zeros(nrows,1);

j = 1; do while j<=rows(knorm);
    cee[j] = ceeofk(knorm[j],gam);
    kp = krange[j]^alp-cee[j]+(1-del)*krange[j];
    kpnorm = (kp-kss)/(wdth*kss);
    cp = ceeofk(kpnorm,gam);

    f[j] = cee[j]^(-phi) - maxc(( bet*(alp*(kp^(alp-1))+1-del)*(cp^(-phi)) )|(kp^(-alp*phi)));
j = j+1; endo;
    
    mtrk = f.*f.*gee;

    retp(sumc(mtrk));
endp;

proc ck(k);
    /* used to calculate the slope of c(k) */
    local  ntees,nks,tees,j,cee;

    ntees = rows(gamopt);
    nks = rows(k);   
    tees = zeros(nks,ntees);
    tees[.,1] = ones(nks,1);
    tees[.,2] = k;
    j=3; do while j<=ntees;
        tees[.,j] = 2*k.*tees[.,j-1]-tees[.,j-2];
    j=j+1; endo;    
    cee = tees*gamopt;

    retp(cee);
endp;