%Chris House, Macro 607, Pset5
%Question 2 f
%An RBC Model with depreciation shocks
%Andrew McCallum, Nathan Seegert, Evan Starr
%April 9, 2008

clc;
clear;
close all;
addpath('U:\Public\html\files\matlab\functions');

%2.f.1
%Enter parameters
alpha = 0.35;
beta = 0.99;
deltabar = 0.025;
rho = 0.2;

%2.f.2
%Steady state values of capital and consumption
kbar = ((1/alpha)*((1/beta)-(1-deltabar)))^(1/(alpha-1));
cbar = kbar*(1-deltabar) - kbar + kbar^(alpha);

%Define the B1 and B2 matrices
B1 = 9*ones(3,3);
B1(1,1:3) = [1 -alpha*beta*(alpha-1)*kbar^(alpha-1) beta*deltabar];
B1(2,1:3) = [0 1 0];
B1(3,1:3) = [0 0 1];

B2 = 9*ones(3,3);
B2(1,1:3) = [1 0 0]; 
B2(2,1:3) = [-cbar/kbar 1/beta -deltabar];
B2(3,1:3) = [0 0  rho];

%...and construct the M matrix
M = inv(B1)*B2;

%2.f.3
%Returns the spectral decomposition of a matrix with our desired ordering
%[L,X,j] = eig_order(M);
%For practice we wanted to code our own 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), delta(t)
m = i;

%2.f.4
%Solve for the policy function
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));

%Generalize dimensions of these matrices for constructing A and B
[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))));

%2.f.5 Construct the VAR form of the solution:
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];

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)];

%Sidetrack: Create a simulated data series using the shocks and A, B
%vectors
%Create fake Y(t)
randn('seed',1234);
Y = 99*ones(10000,3);
epsilonvec = 99*ones(10000,3);
for i = 1:1:size(Y,1);
    epsilon = [zeros(2,1); normrnd(0,0.01)];
    epsilonvec(i,1:3) = epsilon';
    if i == 1
        Y(i,1:3) = A*[0 0 0]' + B*epsilon;
    else
        Y(i,1:3) = A*Y((i-1),1:3)'+B*epsilon;
    end
end

Yfinal = Y((size(Y,1)-999):size(Y,1),(1:3));
evecfinal = epsilonvec((size(Y,1)-999):size(Y,1),(1:3));

%2.f Set the standard deviation of epsilon to 1%....
Vhouse = cvar(A,B,cov(evecfinal));
volhouse = sqrt(diag(Vhouse));
Vus = cov(Y);
volus = (1000/999)*sqrt(diag(Vus));
test=Vus - Vhouse;


%2.f Subject the economy to a 1% depreciation shock
Y = 99*ones(100,3);
for i = 1:1:size(Y,1);
    if i == 1
        epsilon = [zeros(2,1); 0.01];
        Y(i,1:3) = A*[0 0 0]' + B*epsilon;
    else
        Y(i,1:3) = A*Y((i-1),1:3)';
    end
end

%Plot the impulse response functions (for 100 periods instead of 30)
YIRF = Y;
figure;
subplot(2,1,1);
plot(YIRF(:,1),'k-','LineWidth',1.5);
hold on;
plot(YIRF(:,2),'k:','LineWidth',1.5);
plot(zeros(size(YIRF(:,2),1)),'k-','LineWidth',1);
title('Impulse Response')
legend('C tilda','K tilda','Location','SouthEast');
subplot(2,1,2);
plot(YIRF(:,3),'k-','LineWidth',1.5);
title('\delta = 0.01 shock')