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

clc;
clear;
close all;
text('Interpreter','latex');

addpath('\\afs\umich.edu\user\e\c\econandy\Public\html\files\matlab\functions');
%Set parameter values:
alpha = 0.32;
beta = 0.96;
delta = 0.1;
sigma = 2;
rho = 0.95;
omega = 1.455;

%Calculate the nonstochastic steady-state values of the variables
rbar = 1/beta-1+delta;
abar = 1;
gamma = (1-alpha)/(omega+alpha-1);
kbar=(((1/beta)-(1-delta))*(1/(alpha*((1-alpha)*abar)^gamma)))...
    ^(1/(alpha*(1+gamma)-1)); %double checked
nbar = ((1-alpha)*abar*(kbar^alpha))^(1/(omega+alpha-1)); %double checked
ybar = abar*(kbar^(alpha))*(nbar^(1-alpha));
cbar = ybar - delta*kbar;
xbar = cbar - (nbar^omega)/omega;

%Input the A Matrix of coefficients
A_11 = -sigma*(xbar^(-sigma-1))*cbar;
A_12 = beta*(xbar^(-sigma))*rbar*(alpha-1);
A_13 = beta*(xbar^(-sigma))*rbar;
A_14 = beta*(xbar^(-sigma))*rbar*(1-alpha)+sigma*(xbar^(-sigma-1))*(nbar^(omega));
%My variable ordering C, K, A, N, Y, R, W, I;
A(1,1:8) = [A_11 A_12 A_13 A_14 0 0 0 0];
A(2,1:8) = [0 kbar 0 0 0 0 0 0];
A(3,1:8) = [0 0 1 0 0 0 0 0];
A(4,1:8) = [0 0 0 0 0 0 0 0];
A(5,1:8) = [0 0 0 0 0 0 0 0];
A(6,1:8) = [0 0 0 0 0 0 0 0];
A(7,1:8) = [0 0 0 0 0 0 0 0];
A(8,1:8) = [0 0 0 0 0 0 0 0];

%Input the B Matrix of coefficients
B_11 = -sigma*(xbar^(-sigma-1))*cbar;
B_14 = sigma*(xbar^(-sigma-1))*(nbar^(omega));

B_22 = (1-delta)*kbar+alpha*ybar;
B_24 = ybar*(1-alpha);

%My variable ordering C, K, A, N, Y, R, W, I;
B(1,1:8) = [B_11 0 0 B_14 0 0 0 0];  %Euler
B(2,1:8) = [-cbar B_22 0 B_24 0 0 0 0]; %BC
B(3,1:8) = [0 0 rho 0 0 0 0 0];  %A process
B(4,1:8) = [0 alpha 1 (1-alpha) -1 0 0 0]; %Y = F(A,K,N)
B(5,1:8) = [0 -alpha -1 (omega-1+alpha) 0 0 0 0]; %Optimal Labor
B(6,1:8) = [0 (alpha-1) 1 (1-alpha) 0 -1 0 0]; %R = MPK
B(7,1:8) = [0 alpha 1 -alpha 0 0 -1 0]; %W = MPN
B(8,1:8) = [cbar 0 0 0 -ybar 0 0 delta*kbar]; %Y = C+I

%Splitting the A matrix to account for the redundant variables
C_11 = A(1:3,1:3);
C_12 = A(1:3,4:8);
%Splitting the B matrix to account for the redundant variables
D_11=B(1:3,1:3);
D_12=B(1:3,4:8);
D_21=B(4:8,1:3);
D_22=B(4:8,4:8);
%Calculate F which "solves" for the redundant variables
F=-inv(D_22)*D_21;
eig(F*F');
%Calculate the M matrix
M=inv(C_11+C_12*F)*(D_11 +D_12*F);
eig(C_11+C_12*F)
%Put the M matrix into our eigorder function
[L, G, i] = myeigorder(M);

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

%Number of state variables = number of stable eigenvalues, K(t), A(t)
m = i;

%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:8)=zeros(3,5);
Afinal(4:8,4:8)=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:8)=zeros(3,5);
Bfinal(4:8,4:8)=zeros(5,5);

%Subject the economy to a +1% technology shock
Y = 99*ones(100,8);
for i = 1:1:size(Y,1);
    if i == 1
        epsilon = [zeros(2,1); 1; zeros(5,1)];
        Y(i,1:8) = Afinal*zeros(8,1) + Bfinal*epsilon;
    else
        Y(i,1:8) = Afinal*Y((i-1),1:8)';
    end
end

%Plot the impulse response functions for 100 periods
YIRF = Y;
hold on;
plot(YIRF(:,1),'k-','LineWidth',2.5);
plot(YIRF(:,2),'k:','LineWidth',2.5);
plot(YIRF(:,3),'k-.','LineWidth',2.5);
plot(YIRF(:,4),'k-+','LineWidth',2.5);
plot(YIRF(:,5),'k-o','LineWidth',2.5);
title('Log deviation response to +1% shock in A','FontName','Times New Roman','FontSize',20);
legend('C','K','A, \rho = 0.95','N','Y','Location','NorthEast');
set(gca,'FontName','Times New Roman','FontSize',20);
hold off;