% dpexample.m

clear
close all

% Stochastic dynamic programming
% Basic example of value function iteration

% 1. Setting the grid for capital and related information

% Determine stationary capital stock
% from the condition f'(k)=rr+dep;

rr=.06;
dep=.14;
a=1;

% production function is y=ak^gam

gam=.333

kbar=((a*gam)/(rr+dep))^(1/(1-gam));

b=1/(1+rr);

% capital grid

N=500;

index=1:N;

upper=1.1;
lower=.1;

kgrid=kbar*lower+kbar*(upper-lower)*(index/N);

% Gross output at each capital stock level.

ygrid=a*kgrid.^gam + (1-dep)*kgrid;

% 2. utility function curvature
% (and associated plotting parameters)
uopt=1;

if uopt==1;
    
theta=3;
ymin=-20;
ymax=0;

else
    
theta=.9;
ymin=0;
ymax=200;

end

% Initial value function at each capital stock level.

v=(ygrid.^(1-theta))/(1-theta);;

% change in value function

delv=v-0;

% Initial policy (no capital accumulation)

kpol=zeros(N,1);

% Matrix for holding sequence of value functions

it=0;
itmax=20;

% 3. Iterative scheme

% continue iteration so long as change in v is small and number o
% iterations is less than maximum

while abs(max(delv))>.0001 & it<itmax;    
    
% At each point in the k grid

for i=1:N;
    
    cons=ygrid(i)-kgrid;
    % require that consumption is positive 
    testp=cons>0;
    cons=testp.*cons+(1-testp)*(.0000001);
    utility=(cons.^(1-theta))/(1-theta);
    % Calculate maximum value and associated action
    [vnewi poli]=max(utility+b*v);
    vnew(i)=vnewi;
    kpolnew(i)=kgrid(poli);
end

subplot(211),plot(kgrid/kbar,v,'r-',kgrid/kbar,vnew,'k-')
title(['Value Function Updating at Iteration ',num2str(it)])
legend('v','vnew','Location','Southeast')
xlabel('Capital (fraction of ss)')
ylabel('Value')
axis([lower upper ymin ymax])

subplot(212),plot(kgrid/kbar,kpol/kbar,'r-',kgrid/kbar,kpolnew/kbar,'k-')
title(['Policy Function Updating at Iteration ',num2str(it)])
legend('pol','polnew','Location','Southeast')
xlabel('Capital (fraction of ss)')
ylabel('Kprime')

axis([lower upper lower upper])

if (it==1)
    pause
else
    pause(.1)
end

kpol=kpolnew;
v=vnew;
it=it+1;

end


