%Chris House, Macro 607, Pset5
%Question 3
%April 9, 2008

clc;
clear;
close all;
%addpath('U:\Public\html\files\matlab\functions');
addpath('\\afs\umich.edu\user\e\c\econandy\Public\html\files\matlab\functions');

alpha = 0.32;
beta = 0.96;
delta = 0.1;
sigma = 2;
rho = 0.95;
%rho = 0.42; %In Schmitt-Grohe
omega = 1.455;
p = 0.04;
dbar = 0.7442;
psi = 0.000742;

%Calculate the nonstochastic steady-state values of the variables
%From the EE equation
rbar= 1/beta - 1 + delta;
%From the AR(1) process for tech growth
abar = 1;
%From MPL = wage, and rbar
%Original kbar and nbar
%nbar = ((rbar/alpha)*(abar^(alpha/(2*alpha - 1))))^(1/(1-alpha-(omega/alpha)));
%kbar = (abar*nbar^(1-alpha-omega))^(1/alpha);

%New kbar
gamma = (1-alpha)/(omega-1+alpha);
epsilon = 1/(alpha-1+alpha*gamma);
kbar = (rbar/(alpha*abar*((1-alpha)*abar)^gamma))^(epsilon);
nbar = ((1-alpha)*abar*kbar^alpha)^(1/(omega-1+alpha));

ybar = abar*kbar^(alpha)*nbar^(1-alpha);
cbar = ybar - delta*kbar;
xbar = cbar - (nbar^omega)/omega;

%Input the B1 Matrix of coefficients
%Our variable ordering C, K, b, A, w, r, p, N, Y, I;
B1_11 = sigma*beta*rbar*xbar^(-sigma)*cbar;
B1_12 = (1-alpha)*beta*rbar*xbar^(-sigma);
B1_13 = -beta*rbar*xbar^(-sigma);
B1_16 = sigma*xbar^(-sigma-1)*nbar^(-omega)-beta*(1-alpha)*rbar*xbar^(-sigma);

B1(1,1:10) = [B1_11 B1_12 B1_13 0 0 B1_16 0 0 0 0];
B1(2,1:10) = [0 kbar dbar 0 0 0 0 0 0 0];
B1(3,1:10) = [0 0 0 1 0 0 0 0 0 0];
B1(4,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(5,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(6,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(7,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(8,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(9,1:10) = [0 0 0 0 0 0 0 0 0 0];
B1(10,1:10) = [0 0 0 0 0 0 0 0 0 0];

%Input the B2 Matrix of coefficients
B2_11 = sigma*xbar^(-sigma-1)*cbar;
B2_12 = sigma*xbar^(-sigma-1)*nbar^(omega);

B2(1,1:10) = [B2_11 B2_12 0 0 0 0 0 0 0 0];  %euler
B2(2,1:10) = [-cbar kbar/beta (1+p+psi*dbar)*dbar ybar/kbar 0 0 0 0 0 0]; %bc
B2(3,1:10) = [0 0 0 rho 0 0 0 0 0 0];  %A stochastic
B2(4,1:10) = [0 alpha 1 -1 0 -alpha 0 0 0 0]; %pm wage eqn
B2(5,1:10) = [0 -(1-alpha) 1 0 -1 (1-alpha) 0 0 0 0]; %pm rent eqn
B2(6,1:10) = [0 alpha 1 0 0 1-alpha-omega 0 0 0 0];  %labor supply eqn
B2(7,1:10) = [0 alpha 1 0 0 (1-alpha) -1 0 0 0]; % Y=F(k,n,a)
B2(8,1:10) = [cbar 0 0 0 0 0 -ybar delta*kbar 0 0]; %accounting identity
B2(9,1:10) = [0 rbar*(1-alpha) psi*dbar -rbar 0 0 0 -rbar*(1-alpha) 0 0]; %jorgenson
B2(10,1:10) = [0 0 0 0 0 0 -1 1 0 0]; %don't know what to do here

%Splitting the B1 matrix to account for the redundant variables
b1_11 = B1(1:3,1:4);
b1_12 = B1(1:3,5:10);
%Splitting the B2 matrix to account for the redundant variables
b2_11=B2(1:3,1:4);
b2_12=B2(1:3,5:10);
b2_21=B2(4:10,1:4);
b2_22=B2(4:10,5:10);
%Calculate F which "solves" for the redundant variables
F=-inv(b2_22)*b2_21;
%Calculate the M matrix
M=inv(b1_11+b1_12*F)*(b2_11 +b2_12*F);
ttt = eig(M);
%Put the M matrix into our eigorder function
[L, G, i] = myeigorder(M);
%[L, G, i] = eig_order(M);

% Number of control variables = number of unstable eigenvalues , C(t), N(t)
j = size(L,1)-i;
n = j;

%Number of state variables = number of stable eigenvalues, k(t), delta(t)
m = i;
%This mimics our work in 2 which solves to take care of unstable
%eigenvalues
invG = inv(G);
G11 = invG(1:i,1:n);
G12 = invG(1:i,(n+1):(n+m));
G21 = invG((i+1):(i+j),1:n);
G22 = invG((i+1):(i+j),(n+1):(n+m));

H11 = M(1:n,1:n);
H12 = M(1:n,(n+1):(n+m));
H21 = M((n+1):(n+m),1:n);
H22 = M((n+1):(n+m),(n+1):(n+m));

[a1, a2] = size(-inv(G21)*G22*(-H21*inv(G21)*G22+H22));
[a3, a4] = size(-H21*inv(G21)*G22+H22);

s = max(size(9*ones(a1+a3,a2+abs(a4-a2))));
A = 9*ones(s,s);
A(1:a1,1:size(A,2)) = [zeros(a1,size(A,2)-a2) -inv(G21)*G22*(-H21*inv(G21)*G22+H22)];
A((a1+1):size(A,1),1:size(A,2)) = [zeros(a3,size(A,2)-a4) -H21*inv(G21)*G22+H22];
Z=F*A;
%Afinal is the final Matrix on the LHS
Afinal(1:size(M,1),1:size(M,2))=A;
Afinal(size(M,1)+1:size(M,1)+size(F,1),1:size(Z,2))=F*A;
Afinal(1:3,4:10)=zeros(3,5);
Afinal(4:10,4:10)=zeros(5,5);

B = 9*ones(s,s);
B(1:a1,1:size(A,2)) = [zeros(a1,size(A,2)-a2) -inv(G21)*G22];
B((a1+1):size(A,1),1:size(A,2)) = [zeros(a3,size(A,2)-a4) eye(a4,a4)];
%Bfinal is the final Matrix on the RHS 
Bfinal(1:size(B,1),1:size(B,2))=B;
Bfinal(size(B,1)+1:size(F,1)+size(B,1),1:size(B,2))=F*B;
Bfinal(1:3,4:10)=zeros(3,5);
Bfinal(4:10,4:10)=zeros(5,5);

%2.e Subject the economy to a 1% tech shock
Y = 99*ones(100,10);
for i = 1:1:size(Y,1);
    if i == 1
        epsilon = [zeros(3,1); 1; zeros(6,1)];
        Y(i,1:10) = Afinal*zeros(10,1) + Bfinal*epsilon;
    else
        Y(i,1:10) = Afinal*Y((i-1),1:10)';
    end
end

%Plot the impulse response functions (for 100 periods instead of 30)
YIRF = Y;
figure;
%c
plot(YIRF(:,1),'k-','LineWidth',1.5);
hold on;
plot(YIRF(:,2),'k:','LineWidth',1.5);
plot(YIRF(:,6),'k-+','LineWidth',1.5);
plot(YIRF(:,7),'k-o','LineWidth',1.5);
plot(zeros(size(YIRF(:,2),1)),'k-','LineWidth',1);
title('Impulse response to 1% shock in A')
legend('C tilda','K tilda', 'N tilda', 'Y tilda', 'Location','NorthEast');