// LRidentification.ind
//
//  ... contains procedures for identifying structural shocks given
//      the output of a bivariate VARMA or SS Model 
//
//



proc(5)=varma_BQ(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  =varma_BQ(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  estimates the structural matrices of a bivariate ARMA
//  using long-run restrictions and the Cholesky Decomposition 
//  and computes impulse responses
//  and accumulated impulse responses. 
//

local s,p,q,A_1,M_1,longrunmat,longrunmat_struc,mean_yt,
      cont_impact,y_ir,y_ira,yreg,erreg,err,maxlag;

s=rows(ar_mat);
p=cols(ar_mat);
p=p/s;
q=cols(ma_mat);
q=q/s;



A_1=eye(s);
for i (1,p,1);
A_1=A_1-ar_mat[.,((s*(i-1))+1):(s*i)];
endfor;


M_1=eye(s);
for i (1,q,1);
M_1=M_1+ma_mat[.,((s*(i-1))+1):(s*i)];
endfor;

longrunmat=inv(A_1)*M_1;

mean_yt=inv(A_1)*const;

/*********************************************/
/*  Get Structural Parameters                */
/*********************************************/

longrunmat_struc=chol(longrunmat*cov_mat'longrunmat');

cont_impact=inv(longrunmat)*(longrunmat_struc');

longrunmat_struc=longrunmat_struc';

/*********************************************/
/*  Make structural impulse responses        
*/
/*********************************************/

maxlag=maxc(p|q);

y_ir=ones(1,horizon+maxlag+1).*.mean_yt;
err=zeros(s,horizon+maxlag+1);
err[.,maxlag+1]=cont_impact*impuls;

for i (maxlag+1,horizon+maxlag+1,1);
    
    yreg=rev(y_ir[.,(i-p):(i-1)]')';
    yreg=vec(yreg);
    
    erreg=rev(err[.,(i-q):(i-1)]')';
    erreg=vec(erreg);
    
    y_ir[.,i]=const+ar_mat*yreg+err[.,i]+ma_mat*erreg;
    
endfor;

y_ir=y_ir-ones(1,horizon+maxlag+1).*.mean_yt;

y_ir=y_ir[.,maxlag+1:cols(y_ir)];

y_ira=zeros(s,horizon+1);
y_ira[.,1]=y_ir[.,1];
for i (2,horizon+1,1);

    y_ira[.,i]=y_ira[.,(i-1)]+y_ir[.,i];

endfor; 

retp(longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira);
endp;




proc(5)=varma_BQ2(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  =varma_BQ2(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  estimates the structural matrices of a bivariate ARMA
//  using long-run restrictions and the BQ method 
//  and computes impulse responses
//  and accumulated impulse responses. 
//

local s,p,q,A_1,M_1,longrunmat,longrunmat_struc,mean_yt,
      cont_impact,y_ir,y_ira,yreg,erreg,err,maxlag;

s=rows(ar_mat);
p=cols(ar_mat);
p=p/s;
q=cols(ma_mat);
q=q/s;



A_1=eye(s);
for i (1,p,1);
A_1=A_1-ar_mat[.,((s*(i-1))+1):(s*i)];
endfor;


M_1=eye(s);
for i (1,q,1);
M_1=M_1+ma_mat[.,((s*(i-1))+1):(s*i)];
endfor;


longrunmat=pinv(A_1)*M_1;

mean_yt=pinv(A_1)*const;

/*********************************************/
/*  Get Structural Parameters                */
/*********************************************/

cont_impact=getD0(longrunmat,cov_mat);

longrunmat_struc=longrunmat*cont_impact; 

/*********************************************/
/*  Make structural impulse responses        
*/
/*********************************************/

maxlag=maxc(p|q);

y_ir=ones(1,horizon+maxlag+1).*.mean_yt;
err=zeros(s,horizon+maxlag+1);
err[.,maxlag+1]=cont_impact*impuls;

for i (maxlag+1,horizon+maxlag+1,1);
    
    yreg=rev(y_ir[.,(i-p):(i-1)]')';
    yreg=vec(yreg);
    
    erreg=rev(err[.,(i-q):(i-1)]')';
    erreg=vec(erreg);
    
    y_ir[.,i]=const+ar_mat*yreg+err[.,i]+ma_mat*erreg;
    
endfor;

y_ir=y_ir-ones(1,horizon+maxlag+1).*.mean_yt;

y_ir=y_ir[.,maxlag+1:cols(y_ir)];

y_ira=zeros(s,horizon+1);
y_ira[.,1]=y_ir[.,1];
for i (2,horizon+1,1);

    y_ira[.,i]=y_ira[.,(i-1)]+y_ir[.,i];

endfor; 

retp(longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira);
endp;


proc(5)=varma_BQ3(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  =varma_BQ2(const,ar_mat,ma_mat,cov_mat,impuls,horizon);
//
//  estimates the structural matrices of a bivariate ARMA
//  using long-run restrictions and the BQ method 
//  and computes impulse responses
//  and accumulated impulse responses. 
//

local s,p,q,A_1,M_1,longrunmat,longrunmat_struc,mean_yt,
      cont_impact,y_ir,y_ira,yreg,erreg,err,maxlag;
local iA1; 

s=rows(ar_mat);
p=cols(ar_mat);
p=p/s;
q=cols(ma_mat);
q=q/s;



A_1=eye(s);
for i (1,p,1);
A_1=A_1-ar_mat[.,((s*(i-1))+1):(s*i)];
endfor;


M_1=eye(s);
for i (1,q,1);
M_1=M_1+ma_mat[.,((s*(i-1))+1):(s*i)];
endfor;

iA1=(A_1[2,2]~(-A_1[1,2]))|((-A_1[2,1])~A_1[1,1]);
iA1=iA1/(A_1[1,1]*A_1[2,2]-A_1[1,2]*A_1[2,1]); 

longrunmat=iA1*M_1;

mean_yt=iA1*const;

/*********************************************/
/*  Get Structural Parameters                */
/*********************************************/

cont_impact=getD04(longrunmat,cov_mat);

longrunmat_struc=longrunmat*cont_impact; 

/*********************************************/
/*  Make structural impulse responses        
*/
/*********************************************/

maxlag=maxc(p|q);

y_ir=ones(1,horizon+maxlag+1).*.mean_yt;
err=zeros(s,horizon+maxlag+1);
err[.,maxlag+1]=cont_impact*impuls;

for i (maxlag+1,horizon+maxlag+1,1);
    
    yreg=rev(y_ir[.,(i-p):(i-1)]')';
    yreg=vec(yreg);
    
    erreg=rev(err[.,(i-q):(i-1)]')';
    erreg=vec(erreg);
    
    y_ir[.,i]=const+ar_mat*yreg+err[.,i]+ma_mat*erreg;
    
endfor;

y_ir=y_ir-ones(1,horizon+maxlag+1).*.mean_yt;

y_ir=y_ir[.,maxlag+1:cols(y_ir)];

y_ira=zeros(s,horizon+1);
y_ira[.,1]=y_ir[.,1];
for i (2,horizon+1,1);

    y_ira[.,i]=y_ira[.,(i-1)]+y_ir[.,i];

endfor; 

retp(longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira);
endp;



proc(1)=getD0(PhiU,SigU); 
    local a,b,d11,d12,d21,d22,D0;
    local PhuD;     
  

    a=PhiU[1,1]/PhiU[1,2];
    b=2*SigU[1,2]*a+(a^2)*SigU[1,1]+SigU[2,2];

    d12=sqrt((SigU[2,2]*SigU[1,1]-(SigU[1,2]^2))/b);
    d22=-a*d12;
    d11=sqrt(SigU[1,1]-d12^2);

    d21=(SigU[1,2]-d12*d22)/d11;

    D0=(d11~d12)|(d21~d22); 

    // imposing the sign restriction: 
    PhUD=PhiU*D0;
    if PhUD[1,1]<=0;
        D0=-1*D0; 
    endif;
    

    retp(D0); 
endp; 


proc(1)=getD04(PhiU,SigU); 
    local amma,PhiE_11,PhiE_21,PhiE_22,PHiE,D0;

    amma=PhiU*SigU*(PhiU');
    
    PhiE_11=sqrt(amma[1,1]); 
    PhiE_21=amma[1,2]/PhiE_11; 
    PhiE_22=sqrt(amma[2,2]-(PhiE_21^2));
    PhiE=(PhiE_11~0)|(PhiE_21~PhiE_22); 

    D0=inv(PhiU)*PhiE; 

    retp(D0); 
endp; 

// auxiliary procedure: 

proc(1)=reshaper(y,r,c);
// {yt}=reshaper(y,r,c)
// transforms the matrix y to a matrix
// yt with r rows and c columns. It is
// basically a inv(vec(y)) operator since
// it is columnwise. 
local newmat,r_y,c_y;

r_y=rows(y);
c_y=cols(y);

newmat=zeros(r,1);

for i (1,r_y/r,1);

    newmat=newmat~y[(r*(i-1)+1):(i*r),.];

endfor;

newmat=newmat[.,2:cols(newmat)];

retp(newmat);
endp;

