%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   stochasticGrowth.m:  A Matlab program to solve a simple stochastic growth 
%               model 			
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
! rm sgrowth.out
diary sgrowth.out; 
time_start=cputime;
disp('A SIMPLE STOCHASTIC GROWTH MODEL');
disp('');
%
%  set parameter values
%
alpha  = 0.33;              % production parameter
beta   = 0.95;              % subjective discount factor 
prob   = [ .3 .7; .3 .7];   % prob(i,j) = probability (A(t+1)=Aj | A(t) = Ai)
delta  = .92;               % 1 -depreciation rate
A_high = 1.1;               % high value for technology
A_low  = 0.9;               % low value for technology
%
%   form capital grid
%   
mink =   0;                     % minimum value of the capital grid
maxk =  6;                     % maximum value of the capital grid   
ink  =  .01;                     % size of capital grid increments
nk   = round((maxk-mink)/ink+1);   % number of grid points
kgrid = [ mink:ink:maxk ]';


%
%  fill in consumption matrices and single period returns
%

kp = repmat(kgrid,1,nk);
k  = repmat(kgrid,1,nk)';

invest = kp-delta*k;

cons_high = A_high*(k.^(alpha)) - invest;
cons_low  = A_low*(k.^(alpha)) - invest;

util_low    = log(cons_low);
util_high   = log(cons_high);

util_low(find(cons_low<=0)) = -inf;
util_high(find(cons_high<=0)) = -inf;

%
%  initialize some variables
%

v = zeros(nk,2);                  % value functions
v(find(kgrid< 0),:) = -inf;
decis   = zeros(nk,2);            % policy functions 
clear k kp
iter    = 0;
metric1 = 10;
metric2 = 10;
[rs,cs] = size(util_low);
%
%  iterate on Bellman's equation and get the decision 
%  rules and the value function at the optimum         
%
while ((metric1 ~=0) | (metric2 > 0.01));
    
   [tv1,tdecis1]=max(util_low + beta*repmat(v(:,:)*prob(1,:)',1,nk));
   [tv2,tdecis2]=max(util_high + beta*repmat(v(:,:)*prob(2,:)',1,nk));

    tdecis=[tdecis1' tdecis2'];
    tv=[tv1' tv2'];

    metric1 = max(any(tdecis-decis));
    metric2 = max(abs((v(:)-tv(:))./v(:)));
    v=tv;
    decis=tdecis;
    iter=iter+1;
    disp([ iter metric1 metric2 ]);
end;

decis=(decis-1)*ink + mink;

%
%   form transition matrix
%   trans is the transition matrix from state at t (row)
%   to the state at t+1 (column) 
%   The eigenvector associated with the unit eigenvalue
%   of trans' is  the stationary distribution. 
% 

g2=zeros(cs,cs);
g1=zeros(cs,cs);
for i=1:cs
    g1(i,tdecis1(i))=1;
    g2(i,tdecis2(i))=1;
end
trans=[ prob(1,1)*g1 prob(1,2)*g1; prob(2,1)*g2 prob(2,2)*g2];
trans= trans';
probst = (1/(2*nk))*ones(2*nk,1);
test = 1;
while test > 10^(-12);
   probst1 = trans*probst;
   test=max(abs(probst1-probst));
   probst = probst1;
end;
%
%   vectorize the decision rule to be conformable with probst
%   calculate mean level of capital
%
kk=decis(:);
meanK=probst'*kk;
%
%  calculate measure over (k,A) pairs
%  lambda has same dimensions as decis
%
lambda=zeros(cs,2);
lambda(:)=probst;
%
%   calculate stationary distribution of capital 
%
probk=sum(lambda');     
probk=probk';
%
%   print out results
%
disp('PARAMETER VALUES');
disp('');
disp('    alpha      beta       '); 
disp([ alpha beta ]);
disp(''); 
disp('RESULTS ');
disp('');
disp('      mean of K ');
disp([ meanK ]);

time_spent=cputime-time_start;
time_spent

figure(1)
plot(kgrid',v(:,1),'-',kgrid',v(:,2),':');
title('STOCH GROWTH MODEL: VALUE FUNCTION');
print value.ps

figure(2)
plot(kgrid',decis(:,1),'.',kgrid',decis(:,2),':',kgrid',kgrid','-');
title('STOCH GROWTH MODEL: POLICY FUNCTION');
axis([ 0 maxk 0 maxk ]);
print policy.ps

figure(3)
plot(kgrid,probk);
title('DISTRIBUTION OF CAPITAL');
xlabel('CAPITAL');
ylabel('FRACTION OF AGENTS');
print capdist.ps



consopt_high=zeros(1,nk);
consopt_low=zeros(1,nk);

for i=1:nk;
consopt_high(i) = cons_low(tdecis(i,1),i);
consopt_low(i)  = cons_high(tdecis(i,1),i);
end;

figure(4)
plot(kgrid',consopt_high,'-',kgrid',consopt_low,':');
title('STOCH GROWTH MODEL: POLICY FUNCTION FOR CONSUMPTION');

