clc; format short g; clear all; close all;

% calibrate structural parameters
alpha = .36;
delta = .019;
beta = .989;
Phi = 2;
eta = 1;
a = .95;
b = 2.56;
Theta = 1.0125;
n = 1/3;
rho = .95;
sigmae = .007;
gamma = 0;
phi = 0; 
sigmaphi = .0089;


% calculate steady state equilibrium
R = 1/beta;
kn = ((R-1+delta)/alpha)^(1/(alpha-1));
yk = kn^(alpha-1);
ck = yk - delta;
xk = delta;
mc = (a/(1-a))^(-1/b)*((Theta-beta)/Theta)^(-1/b);
mk = mc*ck;

display('The steady state values of [R y/k c/k m/k n/k] =')
[R yk ck mk (1/kn)]

chi = (1+((1-a)/a)*mc^(1-b))^(-1);
Omega1 = chi*Phi + (1-chi)*b;
Omega2 = (b-Phi)*(1-chi);


% compute the rational expectations equilibrium
n1 = 4; n2 = 9; m = n1+n2; e0 = zeros(1,m);

A1 = e(1,m); A2 = e(2,m); A3 = e(3,m); A4 = e(4,m);
A5 = e0; A6 = e0; A7 = e0; A8 = alpha*yk*e(5,m);
A9 = Omega1*e(7,m) - Omega2*e(11,m); A10 = e0;
A11 = e(13,m); A12 = e0; A13 = e0;

B1 = rho*e(1,m); B2 = gamma*e(2,m)+phi*e(1,m); 
B3 = e(9,m); B4 = e(11,m);
B5 = -e(5,m)+alpha*e(3,m)+(1-alpha)*e(6,m)+e(1,m);
B6 = -yk*e(5,m)+ck*e(7,m)+e(9,m)-(1-delta)*e(3,m);
B7 = -delta*e(8,m)+e(9,m)-(1-delta)*e(3,m);
B8 = R*e(10,m)+alpha*yk*e(9,m);
B9 = e(10,m)+Omega1*e(7,m)-Omega2*e(11,m);
B10 = e(5,m)-Omega1*e(7,m)+Omega2*e(11,m)-(1+eta*(n/(1-n)))*e(6,m);
B11 = e(12,m) - e(10,m);
B12 = -e(11,m) + e(7,m) - (1/b)*e(12,m);
B13 = -e(11,m) + e(4,m) - e(13,m) + e(2,m);

A = [A1;A2;A3;A4;A5;A6;A7;A8;A9;A10;A11;A12;A13];
B = [B1;B2;B3;B4;B5;B6;B7;B8;B9;B10;B11;B12;B13];  

% % illustrate QZ decomposition
% [S,T,Q,Z] = qz(A,B)
% A - ctranspose(Q)*S*ctranspose(Z)
% B - ctranspose(Q)*T*ctranspose(Z)
% 
% [Sr,Tr,Qr,Zr] = reorder(S,T,Q,Z);
% abs(ctranspose(Qr)*Sr*ctranspose(Zr) - ctranspose(Q)*S*ctranspose(Z)) 
% abs(ctranspose(Qr)*Tr*ctranspose(Zr) - ctranspose(Q)*T*ctranspose(Z))  

[M,C,numb] = klein(A,B,n1,n2);
N = [1 0; 0 1; 0 0; 0 0];
Sigma = [sigmae^2 0; 0 sigmaphi^2];



% reproduce simulation statistics from Walsh (pp 116-117)
iter = 2000; 
for j = 1:1:iter;
    output1 = artf_data(M,N,Sigma,C,zeros(n1,1),100);           % m.file that generates artificial data
    [smooth,deviation] = hpfilter(output1.data,1600);           % m.file that computes H-P filtered series 
    output2 = sample_moments(deviation,1);                      % m.file that calculates sample moments of series    
    stdv(:,j) = sqrt(diag(output2.autocov(:,:,1)))*100;         % calculating standard deviations
    rel_stdv(:,j) = stdv(:,j)./stdv(1,j);                       % calculating relative standard deviations (y)    
    corr(:,j) = output2.autocorr(:,1,1);                        % calculating contemporaneous correlations with y    
end;

SD = mean(stdv,2);
rel_SD = mean(rel_stdv,2);
CM = mean(corr,2);

display('standard deviations for [y n c x r m i pi] = ')
[SD(1) SD(2) SD(3) SD(4) SD(6) SD(7) SD(8) SD(9)]
display('relative standard deviations for [y n c x r m i pi] = ')
[rel_SD(1) rel_SD(2) rel_SD(3) rel_SD(4) rel_SD(6) rel_SD(7) rel_SD(8) rel_SD(9)]
display('contemporaneous correlations with output for [y n c x r m i pi] = ')
[CM(1) CM(2) CM(3) CM(4) CM(6) CM(7) CM(8) CM(9)]



% % graph the impulse respone functions from Walsh (pp 78-79)
% n1 = 4; n2 = 9; m = n1+n2; e0 = zeros(1,m);
% 
% g = [.5 .8];                % two alternative values for gamma
% 
% T = 31;                     % time horizon for response functions
% IRF1 = zeros(n1,T,2);
% IRF2 = zeros(n2,T,2);
% shock = [0;sigmaphi];
% 
% for j = 1:1:2;
% 
% A1 = e(1,m); A2 = e(2,m); A3 = e(3,m); A4 = e(4,m);
% A5 = e0; A6 = e0; A7 = e0; A8 = alpha*yk*e(5,m);
% A9 = Omega1*e(7,m) - Omega2*e(11,m); A10 = e0;
% A11 = e(13,m); A12 = e0; A13 = e0;
% 
% B1 = rho*e(1,m); B2 = g(j)*e(2,m)+phi*e(1,m);
% B3 = e(9,m); B4 = e(11,m);
% B5 = -e(5,m)+alpha*e(3,m)+(1-alpha)*e(6,m)+e(1,m);
% B6 = -yk*e(5,m)+ck*e(7,m)+e(9,m)-(1-delta)*e(3,m);
% B7 = -delta*e(8,m)+e(9,m)-(1-delta)*e(3,m);
% B8 = R*e(10,m)+alpha*yk*e(9,m);
% B9 = e(10,m)+Omega1*e(7,m)-Omega2*e(11,m);
% B10 = e(5,m)-Omega1*e(7,m)+Omega2*e(11,m)-(1+eta*(n/(1-n)))*e(6,m);
% B11 = e(12,m) - e(10,m);
% B12 = -e(11,m) + e(7,m) - (1/b)*e(12,m);
% B13 = -e(11,m) + e(4,m) - e(13,m) + e(2,m);
% 
% A = [A1;A2;A3;A4;A5;A6;A7;A8;A9;A10;A11;A12;A13];
% B = [B1;B2;B3;B4;B5;B6;B7;B8;B9;B10;B11;B12;B13];
% 
% [M,C,numb] = klein(A,B,n1,n2);
% N = [1 0; 0 1; 0 0; 0 0];
% 
% for i = 1:T; 
%     IRF1(:,i,j) = M^(i-1)*N*shock*100;
%     IRF2(:,i,j) = C*IRF1(:,i,j);
% end;
% 
% end;
% 
% figure,
% plot([0:1:T-1],IRF2(1,:,1),'b-',[0:1:T-1],IRF2(2,:,1),'r-',[0:1:T-1],IRF2(1,:,2),'b--',[0:1:T-1],IRF2(2,:,2),'r--')
% title('Output and Labor')
% ylabel('percent')
% legend('y, \gamma = 0.5','n, \gamma = 0.5','y, \gamma = 0.8', 'n, \gamma = 0.8')
% 
% figure,
% plot([0:1:16],IRF2(8,1:17,1),'k-',[0:1:16],IRF2(9,1:17,1),'m-',[0:1:16],IRF2(8,1:17,2),'k--',[0:1:16],IRF2(9,1:17,2),'m--')
% title('Nominal Interest Rate and Inflation')
% ylabel('percent')
% legend('i, \gamma = 0.5', '\pi, \gamma = 0.5', 'i, \gamma = 0.8', '\pi, \gamma = 0.8')