clear all;

%%% Parameters %%%
beta=0.95; %discount factor%
delta=.05; %depreciation rate%
alpha=.33; %capital share%
gamma=5; %risk aversion coefficient%


%%% steady state %%%

kss=((1/beta-(1-delta))/alpha)^(1/(alpha-1));

%%% Grid %%%

incr=.01; %distance between gridpoints%

kgrid=[incr:incr:1.5*kss];
nk=length(kgrid); %number of gridpoints%


%%% Compute instantaneous utility %%%

u=zeros(nk,nk)+NaN; %set vector for instantaneous utility%

for i=1:nk; %loop for k%
    for j=1:nk; %loop for k'%
    k=kgrid(i);
    kpr=kgrid(j);
    c=k^alpha+(1-delta)*k-kpr; %consumption%
    if c>=0 u(j,i)=(c^(1-gamma)-1)/(1-gamma);
    end;
    end;
end;


%%% Value Function Iteration %%%

V=zeros(1,nk); % initial guess for V(k) % 
dist=10;
iter=1;
while dist>exp(-5) & iter<500;
    W=u+beta*repmat(V',1,nk); %total utility%
    [Vpr, ikpr]=max(W); %take max across k' for each k%
    dist=max(abs(Vpr-V)); %check convergence%
    V=Vpr; %update V(k)%
    iter=iter+1 %iteration%
end

kpr=kgrid(ikpr); %policy function%


%%% Figures %%%

figure(1); %value function%
plot(kgrid,V);
title('Value Function');
xlabel('k');
ylabel('V(k)');

figure(2); %policy function%
plot(kgrid,kpr), hold;
plot(kgrid,kgrid,':'), hold;
title('Policy Function');
xlabel('k');
ylabel('k´');


%%% Simulation %%%

k0=.4*kss; %initial state%

T=100; %number of periods%

ksim=zeros(1,T); %set vector for sequence of capital stocks%

[M,ik]=min(abs(k0-kgrid));

for t=1:T;
    ksim(t)=kgrid(ik);
    ik=ikpr(ik);
end;

csim=ksim(1:T-1).^alpha+(1-delta)*ksim(1:T-1)-ksim(2:T); %simulated consumption%

% Figures % 

figure(3);

subplot(2,1,1);
plot(ksim);
title('capital');

subplot(2,1,2);
plot(csim);
title('consumption');