% ------------------------------------------------------
% This program solves the model and computes 
%  (1) the VARMA representation
%  (2) the State Space represenation
% from the policy rules for the alternative series:
%
%   y_t = [\Delta \log Y_t, \Delta \log L_t]; 
%
%
%
% CK modified 24.10.08
% ------------------------------------------------------
clear;

global model_num;
global model_name;
%--------------------------------------------------------------------------
% SPECIFY THE MODEL: 

model_num=1;        % specifies the model: 
                    % 1 KP benchmark model  
                    % 2 KP w indivisible labor 
                    % 3 KP w Frisch E. 0.63
                    % 4 CKM benchmark model
                    % 5 CKM w indivisble labor
                    % 6 CKM w Frisch E. 0.63
                    % 7 CKM w sigma_l/2
                    % 8 CKM w sigma_l/3
                    
 save_d=0;          % 1 - save VARMA coefficients 0 - dont 
%--------------------------------------------------------------------------
getmodelname; 
model_run;




%--------------------------------------------------------------------------
% DEFINE  PARAMETERS AS IN THE PAPER:


% Get parameters for shocks:
muu=PAR_mu;
rho=PAR_rho;
tau_bar=PAR_taulbar;
sigma_x=PAR_sigmaz;
sigma_l=PAR_sigmatau/tau_bar;


% First Equation, capital

phi_11=hx(1,1);
phi_12=hx(1,3);
phi_1=(1-phi_11)*VAR_K+phi_11*muu-phi_12*tau_bar;


% Second equation, productivity
phi_21=gx(4,1);
phi_22=gx(4,3);
phi_2=VAR_PROD-phi_21*VAR_K+phi_21*muu-phi_22*tau_bar;

% Third Equation, labor
phi_31=gx(3,1);
phi_32=gx(3,3);
phi_3=VAR_L-phi_31*VAR_K+phi_31*muu-phi_32*tau_bar;

% Fourth Equation, detrended output
phi_4=phi_2+phi_3;
phi_41=phi_21+phi_31;
phi_42=phi_42+phi_32;


%--------------------------------------------------------------------------
% STATE SPACE REPRESENTATION 
% ... using the notation in Fernandez-Villaverde, Rubio-Ramirez, Sargent
% and Watson (2006): A,B,C'S (AND D)'S FOR UNDERSTANDING VARS
%


% State transition equation:
% x_t+1 = K_1 + A x_t + B e_t 

K_1=[phi_1-phi_11*muu+phi_12*(1-rho)*tau_bar
    (1-rho)*tau_bar];

A=[phi_11 , phi_12*rho
    0 rho];

B=[-phi_11*sigma_x, phi_12*sigma_l
    0 sigma_l];



% Observation equation:
% y_t = K_2 + C x_t + D e_t 

K_2=[phi_41*inv(phi_11)*phi_1+(1-phi_41)*muu+phi_42*(1-rho)*tau_bar
    phi_3-phi_31*muu+phi_32*(1-rho)*tau_bar];

C=[phi_41*(1-inv(phi_11)), phi_41*inv(phi_11)*phi_12-(1-rho)*phi_42
    phi_31, phi_32*rho];

D=[(1-phi_41)*sigma_x, phi_42*sigma_l
    -phi_31*sigma_x, phi_32*sigma_l];

%--------------------------------------------------------------------------
% GET REDUCED FORM VARMA REPRESENTATION:

    K=B*inv(D);

    AR1=C*A*inv(C);
    MA1=C*(K*C-A)*inv(C);
 
    K_3=K_2+C*K_1-C*A*inv(C)*K_2;

    covmat=D*D';


% ( REDUCED FORM as in Ferndanez-Villaverde et al:
% AR1S=C*((A-K*C)*inv(C)+K);
% MA1S=-C*(A-K*C)*inv(C);)  

%--------------------------------------------------------------------------
% PRINT SOME RESULTS:
clc;
format short; 

disp('Model:'); 
disp(model_name);
dips(' - alternative series - ');

AR1 
MA1
covmat
disp('Impact');
D(:,1)

disp('Eigevalues of AR part:');
eig(AR1)
disp('Eigenvalues of MA part:');
eig(MA1)



%--------------------------------------------------------------------------
% SAVE PARAMETERS:

if save_d==1;

   coeff=[K_3 AR1 MA1 D];
   
    name=['c:\BCAandVARMArepfiles\output\Acoeff_' model_name '.dat'];

    save(name,'coeff','-ascii');

end;

%--------------------------------------------------------------------------

