library pgraph;
pqgwin many;

beta = 0.96;
theta = 0.36;
delta = 0.08;

fn u(c)=ln(c);
fn r(k)=theta*k^(theta-1)-delta;
fn y(k)=k^theta;

maxgap=0.00001;
ngpk=100;           @ No. OF GRID POINTS for CAPITAL @
capt=100;           @ LENGTH OF SIMULATION @

kmin=0.01;          @ BOUNDS for GRID ON CAPITAL @
kmax=10;

@ INITIALIZE VALUE FUNCTION @

v=zeros(ngpk,1);    @ VALUE FUNCTION @
vn=zeros(ngpk,1);
k=zeros(ngpk,1);    @ CAPITAL GRID @
ikn=zeros(ngpk,1);  @ OPTIMAL DECISION RULE @

iseq=zeros(capt,1);

@ CONSTRUCT GRID ON CAPITAL @

k[1]=kmin;
i=2; do while i<=ngpk;
    k[i]=k[i-1]+(kmax-kmin)/(ngpk-1);
i=i+1; endo;

gapv=1;

it=1; do while gapv>=maxgap;

@ TRY EACH POSSIBLE GRID POINT for NEXT PERIOD CAPITAL @
@ BEGIN WITH LOWEST POINT IN THE GRID @

i=1; do while i<=ngpk;

consu=y(k[i])+(1-delta)*k[i]-k[1];
vstar=u(consu)+beta*v[1];
ikstar=1;

j=2; do while j<=ngpk;

consu=y(k[i])+(1-delta)*k[i]-k[j];

if (consu>0.0);
vstarn=u(consu)+beta*v[j];

@ if vstarn BEATS THE OLD vstar, WE HAVE A new WINNER @

    if (vstarn>vstar);
        vstar=vstarn;
        ikstar=j;
    endif;
endif;

j=j+1; endo;

@ save new VALUE for v[i] and ASSOCIATED DECISION RULE @

vn[i]=vstar;
ikn[i]=ikstar;

i=i+1; endo;

@ CHECK WHETHER new VALUE FUNCTION MUCH DIFFERENT TO THE OLD @

gapv=sumc(abs(vn-v));

@ UPDATE VALUE FUNCTION @

v=vn;

it=it+1; endo;

@ run A SHORT SIMULATION @

iseq[1]=10;

i=2; do while i<=capt;
iseq[i]=ikn[iseq[i-1]];
i=i+1; endo;

kseq=k[iseq];

title("capital sequence");
xy(0,kseq);

title("net savings");
xlabel("current capital");
ylabel("k' - k");
xy(k,k[ikn]-k);