%
%  divislabor.m  -- solves Hansen's(1985) Divisible Labor Model
%
%  George Hall
%  Brandeis University
%  October 2006

close all
clear
global TE nx

nx = 3;     % number of state varibles
nu = 2;     % number of control variables 
nz = nx+nu;
%
%  set initial parameter values 
%                       
TE     = 1;
theta  = 0.36;         
delta  = 0.025;         
beta   = 0.99;    
A      = 2;         
rho   = 0.95;         
sigma = 0.007;         
%
%   compute steady state for detrended (but still in levels) economy
%
pstar = [theta; delta; beta; A; rho ];
[zstar] = sstate(pstar);
%
%     Create State Space Representation and map
%     problem into the Optimal Linear Regulator
%     Framework.
%
%        max  X'qX + U'rU + 2U'wX
%
%        s.t  X(t+1) = A1 X(t) + B U(t) + C w(t+1)
%
%     where X' = [ 1  Z(t) K(t) ]
%
%     and  U' = [ I(t+1) H(t) ]
%
[r,dr_dz,d2r_dz2 ] = derret(zstar,pstar);
%
%   Use formula from Hansen, McGrattan & Sargent to get matricies 
%
zSS = zstar;
e = [ 1 0 0 0 0 ]';
M = e*( r - dr_dz'*zSS + .5*zSS'*d2r_dz2'*zSS)*e' + ...
    .5*( e*dr_dz' + dr_dz*e' - e*zSS'*d2r_dz2 - d2r_dz2*zSS*e' + d2r_dz2 );
Qmat = M(1:nx,1:nx);
Rmat = M(nx+1:nz,nx+1:nz);
Wmat = M(1:nx,nx+1:nz);
%
%   form the A,B,C,G,D matricies
%
A1=  [  1         0          0  ;
      (1-rho)    rho         0  ;
        0         0      (1-delta) ];
        
B = [       0         0;
            0         0;
            1         0 ];   

C = [      0    ; 
         sigma ;
           0   ];

%
%    solve for decison rules:  U(t) = -F X(t)
%
[F,P] = olrp(beta,A1,B,Qmat,Rmat,Wmat);
%
%
%
A = A1 - B*F;
%
%    Let's print out some impulse response functions    
%
G = [ A; -F ];

disp('Hansen: Technology Shock' );

C = [      0;
          sigma;
           0 ];

ny = nz;

[Y2sim,X] = bimpulse(A,C,G,zeros(ny,1),1,50);
kss = zSS(3);
iss = zSS(4);
hss = zSS(5);
qss  = 1*kss^(theta)*hss^(1-theta);
css  = qss - iss;
Y2sim(:,6) = (1+Y2sim(:,2)).*((kss+Y2sim(:,3)).^theta).*((hss+Y2sim(:,5)).^(1-theta)) - qss;
Y2sim(:,7) = Y2sim(:,6) - Y2sim(:,4);
Y2sim = 100*Y2sim./repmat([zSS' qss css ],50,1);


figure(1)
plot([2:50]',Y2sim(2:50,[2 5 3]));
title('Hansen Model: Tech. Shock on labor and capital');
xlabel('periods');
ylabel('percentage'); 
legend('lambda','labor','capital');
print -dpsc2 hansen_imp_resp.ps

figure(2)
plot([2:50]',Y2sim(2:50,[2 6 7 4]));
title('Hansen Model: Tech. Shock on q,c,i');
xlabel('periods');
ylabel('percentage'); 
legend('lambda','output','consumption','investment');
print -dpsc2 -append hansen_imp_resp.ps