% This is a sample file for deterministic neoclassical growth theory
% There is no population growth or technology progress
% Labor supply is inelastic
% This is discrete time discrete state space model

% Utility function is util.m
% Production function is cobb.m

% ******Parameter specification******

% Parameter of the model
% risk aversion coefficient 	gamma
% capital share					alpha
% Scale factor in production	A
% discount factor				beta
% depreciation rate				delta

clear all
global gamma alpha A

gamma=1;
A=1;
alpha=0.35;
beta=0.9;
delta=0.1;

% Parameter of computation
% number of grid points				nk
% number of maximum iteration 	maxloop
% tolerence level for value function iteration	tol

nk=50;
maxloop=50;
tol=10^(-1);

% ******Low and upper bound of state variable k ******
% Steady state value 		kstar
% f'(kstar)=1/beta-1+delta

kstar=(((1/(alpha*beta*A))-((1-delta)/alpha*A)))^(1/(alpha-1));
lowk=0.9*kstar;
upk=1.1*kstar;

% ******Vectors and Matrices ******
% vector of the state
stk=(upk-lowk)/(nk-1);
Mk=lowk:stk:upk;
Mk=Mk';
% matric of feasible consumptions and one period utilities
Mu=zeros(nk,nk);
for i=1:nk
   k=Mk(i);
   for j=1:nk
      kprime=Mk(j);
      c=cobb(k)+(1-delta)*k-kprime;
      Mu(j,i)=poweru(c);
     
   end
end

% vectors of value function

Mv=zeros(nk,1);
Mvn=zeros(nk,1);

% vectors of policy function
Polc=zeros(nk,1);
Polkp=zeros(nk,1);

%****** Value Function iteration ******

% Initialize the value function
Mv=poweru(cobb(Mk-delta*Mk));

for kk=1:maxloop
% Maximization
	Mvn=max(Mu+beta*Mv*ones(1,nk))';

% Converge?
	dd=max(abs(Mvn-Mv)./(Mvn));
	if dd<tol

   	break
   end
   Mv=Mvn;
end

% ****** Report the result ******
for i=1:nk
   [Mv(i), pos]=max(Mu(:,i)+beta*Mv);
   Polc(i)=cobb(Mk(i))+(1-delta)*Mk(i)-Mk(pos);
   Polkp(i)=Mk(pos);
end

figure(1)
plot(Mk,Mv)
title('Value Function')
xlabel('capital')
ylabel('value')

print -dps value.ps

figure(2)
subplot(2,1,1)
plot(Mk,Polc)
title('Policy Function-Consumption')
xlabel('capital')
ylabel('consumption')
subplot(2,1,2)
plot(Mk,Polkp,'-',Mk,Mk,'--')
title('Policy Function-Capital')
xlabel('capital')
ylabel('capital next period')

