// Computation of Hankel Matrix using impulse responses

new;

prcsn 80;

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

//------------------------------------------------------------------
// Specifiy model: 

model_num=8;        // 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 SERIES and TRUE PARAMETERS

model_name=getmodelname(model_num);

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

load param[2,7]= ^filename2;

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

aa0=eye(2);
aa1=ar_mat;
mm1=ma_mat; 

//------------------------------------------------------------------
// Compute Hankel matrix of impulse responses and Kronecker indices:  


HH=getIR_Hankelmat(aa0,aa1,mm1,4); 

KK=2;


{pKvec}=getIndices(HH,KK);

print "------------------------------";
print "Parametrization : " model_name;
print "Kronecker Indices: " pKvec;





//--------------------------------------------------------------------------------
//--------------------------------------------------------------------------------

proc(1)=getIR_Hankelmat(AA0,AA,MM,dim_Hankelmat); 
// create Hankel Matrix from ARMA 
// Let A=[A_1,A_2,...,A_p], M=[M_1 M_2, ... M_q]
// collecting the parameters in  
//
// A_0 y_t=A_1 y_t-1 + ... + A_p y_t-p + u_t + M_1 u_t-1 + ... + M_q u_t-q
//
//  INPUT:
//      AA0,AA,MM       coefficient matrices
//      dim_Hankelmat   max. lag computed in the Hankel matrix, 
//                      dim_Hankelmat * dimen. of y gives the dimension
//                      of the Hankelmat 
//
//  change _svdtol=0.001 to alter the threshold for the rank procedure
//
//
    local HH,j,KK,mat; 

    KK=rows(AA0); 

    HH={};

    for j (1,dim_Hankelmat,1);

        mat=varma_irfun_m(AA0,AA,MM,dim_Hankelmat+(j-1));
        mat=mat[.,(j*KK+1):cols(mat)];

        HH=HH|mat; 

    endfor; 

    retp(HH);
endp; 


proc(1)=getIndices(HH,KK);
// 
//
//  INPUT:
//      HH          Hankelmatrix
//      KK          dimension of series
//
//  OUTPUT:
//      pKvec       vector of Kronecker Indices
//
local pKvec,HHn,n,k,nrank,orank,m;

    // Vector of Missing Values
    pKvec=zeros(1,KK); pKvec=miss(pKvec,0);

   
    HHn={};
    m=rows(HH)/KK; 
    
    for n (1,(m),1); // loop over LAGS in time 
    

        for k (1,KK,1); // loop over elements in y_t

            if isinfnanmiss(pKvec[k]); // check wether Kronecker Index is shorter 
                                       // than time index!
                
                if rows(HHn) ne 0;     // compute old rank
                    orank=rank(HHn);
                else; orank=0;
                endif;            

                HHn=HHn|HH[(KK*(n-1)+k),.]; // Add another row
                
                nrank=rank(HHn);            // compute new rank
                
                if (nrank-orank)==0;        // Check whether the rank increased.    
                    if (rows(HHn)-1)>0;                    
                    HHn=HHn[1:(rows(HHn)-1),.];
                    else; HHn={};           // if it did not take away the row
                    endif;
                    pKvec[k]=n-1;           
                endif;                      
           
            endif;


        endfor;

    endfor;


retp(pKvec);
endp;


proc(1)=varma_irfun_m(AA0,AA,MM,h);
// irfun=varma_ifun(AA0,AA,MM,h);
// computes the impulse response function
// of a VARMA (p,q) model as
// AA0 y_t = AA1 y_(t-1) + ... AA0 u_t + ... + MMq u_(t-q);
//
// INPUTS:
//          AA0     KK x KK 
//          AA      KK x KKp    AR coefficients
//                              AA=[AA1: ... : AAp]
//          MM      KK x KKq    MA coefficients
//                              MM=[MM1: ... : MMq]
//          h       1 x 1       horizon
//
//          Remark: If p or q=0 then AA={}, MM={}, respectively
//
// OUTPUTS:
//          irfun   KK^2 x (h+1)Matrix of vectorized ir
//                              functions for j=0,...,h 
//
// See Luetkepohl (2005) pp. 426 for Documentation.
local AAA,A11,A12,A21,A22,p,q,KK,JJ,HH,AAAJ,irfun;

KK=rows(AA0);

p=cols(AA)/KK;
q=cols(MM)/KK;

// Construction of AAA:

A11=zeros(KK,KK);
A12=zeros(KK,KK);
A21=zeros(KK,KK);
A22=zeros(KK,KK);


if p ne 0;
A11=AA;
endif;

if p>1;
A11=A11|(eye(KK*(p-1))~zeros((p-1)*KK,KK));
endif;

if q ne 0;
A12=MM;
    if p>1;
        A12=A12|(zeros(KK*(p-1),KK*q));//
    endif;
endif;

if (q ne 0);
    A21=zeros(KK*q,KK);
endif;

if (q ne 0) and (p>0);
    A21=zeros(KK*q,KK*p);
endif;

if q ne 0;
A22=zeros(KK,KK*q);
endif;

if q>1;
A22=A22|(eye(KK*(q-1))~zeros(KK*(q-1),KK));
endif;

if q==0;
A12={};
A21={};
A22={};
endif;

AAA=(A11~A12)|(A21~A22);




JJ=eye(KK);
if p>1;
    JJ=JJ~(zeros(KK,KK*(p-1)));
endif;

if q>0;
    JJ=JJ~zeros(KK,KK);
    if q>1;
    JJ=JJ~(zeros(KK,KK*(q-1)));
    endif;
endif;



HH=eye(KK);

if p>1;
    HH=HH|zeros(KK*(p-1),KK);
endif;

if q > 0;
    HH=HH|eye(KK);
    if q>1;
    HH=HH|(zeros(KK*(q-1),KK));
    endif;
endif;


irfun={};
AAAj=eye(cols(AAA));

for j (0,h,1);

    irfun=irfun~(JJ*AAAJ*HH);
    AAAJ=AAAJ*AAA;

endfor;


retp(irfun);
endp;