// checkIdentfication.prn 
// 
// The program computes the Sylvester Tensor Matrix 
// for the VARMA(1,1) representation of the DSGE model
// and a "perturbed" version (setting one eigenvalue
// of the MA part equal to one of AR part) and checks
// the rank of these matrices.  
//
//
// 
// Authors: Christian Kascha and Karel Mertens 
//


new;

#include "c:\\BCAandVARMArepfiles\\GAUSS\\getmodelname.ind";



//------------------------------------------------------------
// SETTINGS:

   
model_num=4;        // 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



//------------------------------------------------------------
// GET PARAMETERS

model_name=getmodelname(model_num);

filename="c:\\BCAandVARMArepfiles\\output\\"$+"coeff_"$+model_name$+".dat";


load param[2,7]= ^filename;

const=param[.,1];
ar_mat=param[.,2:3];
ma_mat=param[.,4:5];
cov_mat=param[.,6:7];
cov_mat=cov_mat*(cov_mat');

// Compute Sylvester Tensor Matrix: 

sba=((-eye(2)).*.eye(2))~((-ma_mat).*.eye(2));

sba=sba|((eye(2)).*.eye(2))~(eye(2).*.(-ar_mat));

 


//--------------------------------------------------------------------------------
// Construct perturbed process: 

{va,ve}=eigv(ar_mat);

va=diagrv(eye(2),va);
va=real(va); 
ve=real(ve); 

{va2,ve2}=eigv(ma_mat);
va2=real(va2); 
ve2=real(ve2); 

va2=diagrv(eye(2),real(va2));
va2_w=va2;
va2_w[1,1]=-va[1,1]; 
va_w=va;

ar_mat2=ve*va_w*inv(ve); 
ma_mat2=ve2*va2_w*inv(ve2); 


// Compute New Sylvester Tensor Matrix: 

sba2=((-eye(2)).*.eye(2))~((-ma_mat2).*.eye(2));

sba2=sba2|((eye(2)).*.eye(2))~(eye(2).*.(-ar_mat2));

print "-------------------------------------------";
print "Parametrization: " model_name; 
print "Rank of Original Process:";
rank(sba);
print "Det: " det(sba); 
print ""; 
print "Rank of Pertubed Process:";
rank(sba2);

print "Det: " det(sba2);
 