%  optgrowth2.m : Solution of the stochastic optimal growth problem
%  R. Chang, 10/5/2005, 2005
% Edited, oct. 2009


clear;

% Parameters

a = 0.36;   % capital share in prod function (alpha)
b = 0.95;   % discount factor (beta)
d = 0.025;  % depreciation rate (delta)
rho = 0.95; % autocorrelation of productivity factor
s1 = 2;     % coefficient of relative risk aversion (sigma)


% Steady state values
% We assume that the production function is Cobb Douglas and 
% capital depreciation is incomplete

kss = (a / ( (1/b) - (1-d)) )^(1/(1-a)) ;   % ss capital

ky = kss^(1-a) ;                            % ss capital output ratio

cy = 1 - d*ky ;                             % ss consumption output ratio



% Now write the model as aa * E_t X(t+1) = bb * X(t) 
%   where X(t) = (k(t), s(t), c(t))
% Note that the states (both endog and exog) come first in X(t)

aa = zeros(3,3); bb = zeros(3,3);

% Eq. 1:  ky k(t+1) = s(t) + (a + (1-d)ky) k(t) - cy c(t)

aa(1,1) = ky;
bb(1,1) = a + (1-d)*ky;
bb(1,2) = 1;
bb(1,3) = -cy;

% Eq. 2: b Et {-s1 c(t+1) + (1-b*(1-d))*(s(t) +(a-1) k(t+1) )} = -s1 c(t)
aa(2,1) = (1-b*(1-d))*(a-1);
aa(2,2) = (1-b*(1-d));
aa(2,3) = -s1;
bb(2,3) = -s1;

% Eq. 3: E_t s(t+1) = rho s(t)
aa(3,2) = 1;
bb(3,2) = rho;


% Now we simply call the solab procedure, which calculates the solution as
%
%       J(t) = ff*S(t) ...(1)
%       E_t S(t+1) = gg*S(t)...(2)

% where S(t) = (k(t),s(t))' is the vector of state variables
%   and J(t)= (c(t))' the vector of jumping variables

[ff gg] = solab(aa,bb,2);

% Now we compute impulse responses 

% First, examine the adjustment process if capital is above steady state

 sresponses = zeros(3, 40) ;  % This will hold the responses of the three variables 

 sresponses(1,1) = 1; % set the initial value of k to one
 
 % Use (2) recursively to generate the response path of S(t)

 for tcounter = 2:40 ,
    sresponses(1:2,tcounter) = gg*sresponses(1:2, tcounter-1);
 end;

 % Use (1) to generate the path of J(t)
 
 for tcounter = 1:40,
  sresponses(3, tcounter) = ff*sresponses(1:2, tcounter);
end;    

% Finally. graph the responses

t_axis = 1:1:40;


plot(t_axis,sresponses(1,:),'-',t_axis,sresponses(2,:),'o', t_axis,sresponses(3,:),'* ');

title('Adjustment Path to capital initially outside steady state') 
legend('capital', 'productivity', 'consumption','Location','SouthOutside')


disp('Inspect figure. Hit key when ready...');
            pause;

% Now, examine impulse response to a productivity shock of one unit

sresponses = zeros(3, 40) ;  % This will hold the responses of the three variables 

sresponses(2,1) = 1; % set the initial value of s(t) to one

for tcounter = 2:40 ,
    sresponses(1:2,tcounter) = gg*sresponses(1:2, tcounter-1);
end;

for tcounter = 1:40,
    sresponses(3, tcounter) = ff*sresponses(1:2, tcounter);
end;    


t_axis = 1:1:40;

plot(t_axis,sresponses(1,:),'-', t_axis,sresponses(2,:),'o', t_axis,sresponses(3,:),'*');

title('Response to productivity shock')
legend('capital', 'productivity', 'consumption', 'Location','SouthOutside')