%  FARMER.M solves problem 4 in Farmer, p.61

%  R. Chang, 9/26, 2005


clear;

% Parameters

r = 0.0163;
a = 0.4;
lam = 0.95;
t = 0.9;
B = 0.9957;
d = 0.0272;

gss = B^(1/(1-t));

x = (a / ( r + d ) )^(1/ (1-a));

lss = ( (1-a)*(x^a) + gss)/ (x^a - d*x );

cy = (1-a)/lss;

gy = 1 - cy - d*(x^(1-a)) ; 

% Now write the model as 

% aa * E_t X(t+1) = bb * X(t) 

%   where X(t) = (k(t), a(t), g(t), c(t), i(t), l(t))

% Recall that state (endogenous and exogenous) variables go first in X(t),
% ie X(t) = (S(t)', J(t)')' in the notation below

aa = zeros(6,6); bb = zeros(6,6);

aa(1,4) = -1;
aa(1,2) = 1 - (1+r)*(1-d);
aa(1,6) = (1-a)*aa(1,2);
aa(1,1) = -aa(1,5);
bb(1,4) = -1;

bb(2,2) = 1;
bb(2,1) = a;
bb(2,6) = - a ;
bb(2,4) = -1;

bb(3,4) = cy;
bb(3,5) = d*x^(1-a);
bb(3,3) = gy;
bb(3,2) = -1;
bb(3,1) = -1;
bb(3,6) = - 1 + a;

aa(4,1) = 1;
bb(4,1) = 1-d;
bb(4,5) = d;

aa(5,2) = 1;
bb(5,2) = lam;

aa(6,3) = 1;
bb(6,3) = t;


% Now we simply call the solab procedure, which calculates the solution as
%
%       J(t) = ff*S(t)
%       E_t S(t+1) = gg*S(t)

% where S(t) = (k(t),a(t),g(t))' is the vector of state variables
%   and J(t)= (c(t), i(t), l(t))' the vector of jumping variables

[ff gg] = solab(aa,bb,3);

% Now we compute impulse responses to a (one st. dev.)  tech shock 

sresponses = zeros(3, 40) ;  % This will hold the responses of the state variables 

sresponses(2,1) = sqrt(0.07);

for tcounter = 2:40 ,
    sresponses(:,tcounter) = gg*sresponses(:, tcounter-1);
end;

t_axis = 1:1:40;

plot(t_axis,sresponses(1,:),t_axis,sresponses(2,:));

title( 'Responses of state variables')
legend('capital', 'technology')
pause

close all


Jresponses = zeros(3,40); % This matrix will hold the responses of jumping variables

for t_counter = 1:40,
    Jresponses(:,t_counter) = ff*sresponses(:,t_counter);
end;

plot(t_axis,Jresponses(1,:),t_axis,Jresponses(2,:),t_axis,Jresponses(3,:));
title( 'Responses of jumping variables')
legend('consumption', 'investment','hours')

