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

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

alpha = 0.35;
beta = 0.99;
delta = 0.025;
rho = 0.95;
%Initial state variable values

b=1;
rbar= 1/beta - 1 +delta;

%Input the B1 Matrix of coefficients
B1 = 9*ones(8,8);
B1(1,1:8) = [1 -beta*rbar*(alpha-1) -beta*rbar 0 0 -beta*rbar*(1-alpha) 0 0];
B1(2,1:8) = [0 1 0 0 0 0 0 0];
B1(3,1:8) = [0 0 1 0 0 0 0 0];
B1(4,1:8) = [0 0 0 0 0 0 0 0];
B1(5,1:8) = [0 0 0 0 0 0 0 0];
B1(6,1:8) = [0 0 0 0 0 0 0 0];
B1(7,1:8) = [0 0 0 0 0 0 0 0];
B1(8,1:8) = [0 0 0 0 0 0 0 0];

%Input the B2 Matrix of coefficients
B2 = 9*ones(8,8);
B2(1,1:8) = [1 0 0 0 0 0 0 0];
B2(2,1:8) = [(-rbar/alpha + delta) 1/beta rbar/alpha 0 0 (rbar/alpha)*(1-alpha) 0 0];
B2(3,1:8) = [0 0 rho 0 0 0 0 0];
B2(4,1:8) = -1*[0 (1-alpha) -1 1 0 (alpha-1) 0 0];
B2(5,1:8) = -1*[0 -alpha -1 0 1 alpha 0 0];
B2(6,1:8) = -alpha*[b/alpha -1 -1/alpha 0 0 1 0 0];
B2(7,1:8) = -1*[0 -alpha -1 0 0 (alpha-1) 1 0];
B2(8,1:8) = -1*[((rbar/(alpha*delta))-1) -rbar/delta -rbar/(alpha*delta) 0 0 -rbar*(1-alpha)/(delta*alpha) 0  1];

%Splitting the B1 matrix to account for the redundant variables
b1_11 = B1(1:3,1:3);
b1_12 = B1(1:3,4:8);
%Splitting the B2 matrix to account for the redundant variables
b2_11=B2(1:3,1:3);
b2_12=B2(1:3,4:8);
b2_21=B2(4:8,1:3);
b2_22=B2(4:8,4:8);
%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);
%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)
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: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);

randn('seed',1234);
Y = 99*ones(1000,8);
epsilonvec = 99*ones(1000,8);
for i = 1:1:size(Y,1);
    epsilon = [zeros(2,1); normrnd(0,1);zeros(5,1)];
    epsilonvec(i,1:8) = epsilon';
    if i == 1
        Y(i,1:8) = Afinal*[0 0 0 0 0 0 0 0]' + Bfinal*epsilon;
    else
        Y(i,1:8) = Afinal*Y((i-1),1:8)'+Bfinal*epsilon;
    end
end

Yfinal = Y((size(Y,1)-999):size(Y,1),(1:8));
efinal = epsilonvec((size(Y,1)-999):size(Y,1),(1:8));

%Similulated based on fake Y
Vus = cov(Yfinal);
VCV = cov(efinal);
Vhouse = cvar(Afinal,Bfinal,VCV);
test = Vus - Vhouse;

%Analytical Moments
AVCV=zeros(8,8);
AVCV(3,3) = 1;
AVhouse = cvar(Afinal,Bfinal,AVCV);

% e)
%Now we get the simulated data for Y and not just growth rate.  
Cbar=(1-alpha)*(rbar/alpha)^(alpha/(alpha-1));
Kbar= Cbar/(rbar/alpha - delta);
Wbar = Cbar;
Nbar = Kbar*(rbar/alpha)^(1/(1-alpha));
Ybar = Kbar^(alpha)*Nbar^(1-alpha);
Ibar = delta*Kbar;
Ytotalbar=[Cbar Kbar 1 rbar Wbar Nbar Ybar Ibar];
logYlevels=repmat(Ytotalbar,1000,1).*(1+Yfinal);

%Now compare the hpfiltered data with the TRUE log deviations. 
mu = 1600;
stalogy = 9*ones(1000,8);
for i = 1:1:8;
    [tau, F2] = myHPfilter(logYlevels(:,i),mu);
    stalogy(:,i) = logYlevels(:,i) - tau;
end
date = 1:1000;
z = zeros(1000,1);

Difference=Yfinal-stalogy;
titlevec = {'C(t)' 'K(t)' 'A(t)' 'w(t)' 'r(t)' 'N(t)' 'Y(t)' 'I(t)'};
figure;
hold on;
for i= 1:8
    subplot(3,3,i)
    plot(date,Difference(:,i),'k-');
    hold on;
    plot(date,zeros(1000,1),'k-');
    title(titlevec(1,i));
   axis([0 1000 -40 40]);
end
% annotation(gcf,'textbox','FitBoxToText','off','Position',...
%     [3 3 1 0.1],'LineStyle','none','String',...
%     {'Log deviations from SS minus filtered simulated variables'});

%Now calculate the moments of the detrended log variables and compare them
%to the analytical moments from cvar.m
VHP = cov(stalogy);
covVHP = eig(VHP)
covVUS = eig(Vus)

covAVhouse = eig(AVhouse)
covVhouse = eig(Vhouse)