%--------------------NON-STOCHASTIC GROWTH MODEL---------------------------
%%% THIS FILE SOLVES THE NON-STOCHASTIC GROWTH MODEL PRESENTED
%%%  IN LECTURE TWO OF RECURSIVE MACROECONOMICS OF MARTIN ELLISON 
%%% A SIMILAR VERSION IS AVAILABLE AT http://www.eco.utexas.edu/~cooper/dynprog/dynprog1.html
%--------------------------------------------------------------------------
%The production function is kapa to the power of alfa
%The utility function is CRRA with signa the coefficient of relative risk aversion  
%It includes transition dynamics only for K

%--------------------------------------------------------------------------
%                           DEFINE PARAMETERS
%--------------------------------------------------------------------------

clear all

     sigma=1;
     beta=.9;
     alpha=.75;
     delta=0.3;                %depreciation factor

N=100;
Kstar=(((1/(alpha*beta))-((1-delta)/alpha)))^(1/(alpha-1));
Klo=Kstar*0.9;
Khi=Kstar*1.1;

step=(Khi-Klo)/N;
K=Klo:step:Khi;

n=length(K);                   % n is the true length of the state space

%--------------------------------------------------------------------------
%               CONSTRUCT FIRST GUESS OF THE VALUE FUNCTION
%--------------------------------------------------------------------------

ytot = K.^alpha+(1-delta)*K;   % first guess for consumption= product of capital plus its depreciated value 
ytot=ytot'*ones(1,n);          % put the vector above as columns in a matrix 

if sigma==1
    v=log(ytot);               % first guess of the value function- the level of utility associated with consuming current capital immediately
else
    v=ytot.^(1-sigma)/(1-sigma);
end

%--------------------------------------------------------------------------
%               CALCULATE CONSUMPTION AND GET ASSOCIATED UTILITY
%--------------------------------------------------------------------------
for i=1:n
    for j=1:n
        C(i,j) = K(j)^alpha + (1-delta)*K(j) - K(i);
    end
end

if sigma==1
    U=log(C);
else
    U=(C.^(1-sigma))/(1-sigma);
end

%--------------------------------------------------------------------------
%                          VALUE FUNCTION ITERATION
%--------------------------------------------------------------------------

T=100
for j=1:T
        w=U+beta*v;
        v1=max(w);
        v=v1'*ones(1,n);
end
[val,ind]=max(w);
optk = K(ind);

%--------------------------------------------------------------------------
%                              GRAPH RESULTS
%--------------------------------------------------------------------------

figure(1)
plot(K,optk','LineWidth',2)
hold on
plot(K,K','--r',...
    'LineWidth',2)
xlabel('K');
ylabel('K`');
legend('Policy function','45 degree line',4);
text(0.4,0.65,'45 degree line','FontSize',18)
text(0.4,0.13,'K`','FontSize',18)

figure(2)
plot(K,(optk-K)','LineWidth',2)
xlabel('K');
ylabel('K` - K');

%--------------------------------------------------------------------------
%                           TRANSITION DYNAMICS
%--------------------------------------------------------------------------

p=50;
mi=zeros(p,1);
m=zeros(p,1);
mi(1)=N;
m(1)=K(mi(1));
for i=2:p
    mi(i)=ind(mi(i-1));
    m(i) = K(mi(i));
end
t=1:1:50

figure(3)
plot(t,m)
xlabel('t');
ylabel('k(t)');