

proc(3)=varma2SLSp(y,p,q,RR,r,nT,esd,ARmethod);
// {gamma,vc_res,res}=varma_2sls(y,p,q,RR,r,nT,ed);
// computes the two-stage least squares estimator
// of a VARMA(p,q):
//
// A_0 y_t = A_1 y_t-1 + ... + A_0 u_t + ...+ M_q u_(t-q)
//
// Restrictions have to be given as described in
// Lütkepohl 2005:
//
// vec [I-A_0 : A : M]=RR gamma + r  
//
// INPUT:
//      y           KK x TT     series
//      p           1 x 1       AR - order
//      q           1 x 1       MA - order
//      RR                      Restriction matrices
//      r                       as defined above
//      nT          1 x 1       Order of long autoregression
//      esd         1 x 1       Dummy: esd=0 - OLS
//                              esd=1 - GLS 
//
//      ARmethod    1 x 1       Method of AR estimation:
//                              1 - OLS 2 - OLS on demeaned data 
//                              3 - YW
//
// OUTPUT:
//      gamma_hut               parameter estimates
//      vc_res                  estimated covariance matrix
//      res                     estimated residuals
//
// CALLS:
//
//      reshaper
//
//  ADD COVARIANCE PARAMETERS
//

local K,Z,Nobs,T_eff,T_eff2,YY,Z_t,PInT,UU0,res,vc_res,
      cov_est,mpq,YY2,XXcnT,U_t,y_tc,gamma_hut,i;
local ar_e;  

// Dimension of Model
    K=rows(y);
    Nobs=cols(y);
    T_eff=Nobs-nT;
    mpq=maxc(p|q);
    T_eff2=T_eff-mpq;

// FIRST STAGE LEAST SQUARES:

    YY=y[.,(nT+1):Nobs];

    if ARmethod==1;
        {AR_e,cov_est,UU0}=var_ols(y,0,nT);
    elseif ARmethod==2;
        {AR_e,cov_est,UU0}=var_dm(y,nT);
    elseif ARmethod==3;
        {AR_e,cov_est,UU0}=var_yw(y,nT);
    endif; 
    

// SECOND STAGE LEAST SQUARES:
// Setting up matrices:

    YY2=YY[.,(mpq+1):T_eff];

    XXcnT={};

    for i (1+mpq,T_eff,1);
        Z_t=(rev(YY[.,(i-p):(i-1)]'))';
        Z_t=vec(Z_t);
        U_t=(rev(UU0[.,(i-q):(i-1)]'))';
        U_t=vec(U_t);
        y_tc=YY[.,i]-UU0[.,i];
        XXcnT=XXcnT~(y_tc|Z_t|U_t); 
    endfor;

    z=vec(YY2)-(XXcnT'.*.eye(K))*r;

    if esd == 0;
        cov_est=eye(K);
    endif;

    gamma_hut=inv(RR'((XXcnT*(XXcnT')).*.inv(cov_est))*RR)*RR'(XXcnT.*.inv(cov_est))*z;

    res=z-((XXcnT').*.eye(K))*RR*gamma_hut; 
    res=reshaper(res,K,T_eff2);
    vc_res=moment(res',0)/Nobs;

retp(gamma_hut,vc_res,res);
endp;


//-------------------------------------------------
// AUXILIARY PROCEDURES: 
//
//  var procedures
//  procedures to check the eigenvalues
//  of the estimated polynomials



proc(3)=var_yw(y,p);
    /*
    ESTIMATION     VAR by Yule-Walker

    USAGE:
    {AR_e,COV_e,RES_e}=var_yw(y,p);
    
    INPUTS
    y       k x T   series
    p       1 x 1   lag order

    OUTPUT
    a       k x k*p+1   parameters
    covm    k x k       covariance matrix
    u                   residuals

    */
local a,T,i,covm,u,ydm,k,yd;
local G0,GP,ys,YB,meany,xx,xx0,Te;

    T=cols(y);  // sample observations
    k=rows(y);
    Te=T-p; 

    meany=sumc(y')/(T);
    
    ydm=y-meany;
    
    ys=trimr(ydm',p,0)';

    if p ne 0; 

    YB=rev(ydm[.,1:(1+p-1)]');
    YB=vec(YB');
    for i (2,T-p,1);
        xx=rev(ydm[.,i:(i+p-1)]');
        xx=vec(xx');
        YB=YB~xx;
    endfor;

    // Making covariance matrices: 

    
    xx0={};
    for h (0,(p-1),1); 

        yd=(ydm[.,(1+h):T])*((ydm[.,(1):(T-h)])');
        xx=(yd)/(T-h);         
        xx0=xx0~xx;
 
   endfor; 

   G0=xx0;
   xx=xx0;
   for h (1,(p-1),1); 
        xx=(xx0[.,(k*h+1):(k*(h+1))]')~xx[.,1:(cols(xx)-k)]; 
        G0=G0|xx;
   endfor; 


    GP={};
    for h (1,p,1); 
    
        yd=(ydm[.,(1+h):T])*((ydm[.,(1):(T-h)])');
        xx=(yd)/(T-h);         
        GP=GP~xx;
    endfor;  

    a=GP*inv(G0);
    
    u=(ys-a*YB);
    else;

    u=ydm;
    a=meanc(y[.,p+1:T]');   

    endif; 



    covm=moment(u',0)/(Te-k*p-1);

retp(a,covm,u);
endp;


proc(3)=var_dm(y,p);
    /*
    ESTIMATION     VAR by substracting the sample mean
                   prior to estimation

    USAGE:
    {AR_e,COV_e,RES_e}=var_ols(y,p);
    
    INPUTS
    y       k x T   series
    p       1 x 1   lag order

    OUTPUT
    a       k x k*p+1   parameters
    covm    k x k       covariance matrix
    u                   residuals

    */
local a,ys,yb,T,xx,i,covm,u,ydm,k,Te;

    T=cols(y);
    k=rows(y);
    Te=T-p; 
    ydm=y-meanc(y[.,p+1:T]');
    
    ys=trimr(ydm',p,0)';

    if p ne 0; 
  
    YB=rev(ydm[.,1:(1+p-1)]');
    YB=vec(YB');
    for i (2,T-p,1);
        xx=rev(ydm[.,i:(i+p-1)]');
        xx=vec(xx');
        YB=YB~xx;
    endfor;

    a=(ys')/(YB');
    
    a=a';
    u=(ys-a*YB);
    else; 

    u=ydm;
    a=meanc(y[.,p+1:T]');   

    endif; 


    covm=moment(u',0)/(Te-k*p-1);

retp(a,covm,u);
endp;




proc(3)=var_ols(y,mud,p);
/*
    ESTIMATION     VAR by OLS

    USAGE:
    {AR_e,COV_e,RES_e}=var_ols(y,mud,p);
    
    INPUTS
    y       k x T   series
    mud     1 x 1   dummy: 1 -> constant, 0 -> no constant
    p       1 x 1   lag order

    OUTPUT
    a       k x k*p+1   parameters
    covm    k x k       covariance matrix
    u                   residuals



*/
local a,cvm,ys,yb,T,xx,i,covm,u,k,Q,R,Te;


T=cols(y);
ys=trimr(y',p,0)';
k=rows(y);
Te=T-p; 

if p ne 0; 
    YB=rev(y[.,1:(1+p-1)]');
    YB=vec(YB');
    for i (2,T-p,1);
        xx=rev(y[.,i:(i+p-1)]');
        xx=vec(xx');
        YB=YB~xx;
    endfor;


if mud==1;
yb=ones(1,(cols(y)-p))|yb;
else;
yb=yb;
endif;

a=(ys')/(YB');
a=a';
u=(ys-a*YB);

else; 
    if mud==1;
        yb=ones(k,cols(ys));     
        a=meanc(ys');
        u=(ys-a);
    else; 
        a=0;
        u=ys;
    endif; 
endif; 



covm=moment(u',0)/(Te-k*p-mud);


retp(a,covm,u);
endp;

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;





proc(1)=check_maxroot(MM);
// {mroot}=check_maxroot calculates
// the maximum eigenvalue of the companion matrix
// to the matrix polynomial
// M(z) = I + M_1 z + .. + M_q z 
// 
// INPUTS:
//      MM = [M_1 : M_2 : ... : M_q]
//      (K x Kq)
//
// OUTPUTS:
//      mroot
//
// CALLS:
//      - 
// COMMENT: The Eigenvalues of the companion
//          matrix correspond to the negative, 
//          inverse roots of M(z), lambda= - 1/z
//
local mroot,compmat,K,q;

K=rows(MM);
q=cols(MM)/K;

if q > 1;
compmat=MM|(-eye(K*(q-1))~zeros(K*(q-1),K));
else;
compmat=MM;
endif;

mroot=maxc(abs(eig(compmat)));

retp(mroot);
endp;


proc(1)=check_maxroot2(AA);
// {mroot}=check_maxroot calculates
// the maximum eigenvalue of the companion matrix
// to the matrix polynomial
// M(z) = I - M_1 z - .. - M_q z 
// 
// INPUTS:
//      MM = [M_1 : M_2 : ... : M_q]
//      (K x Kq)
//
// OUTPUTS:
//      mroot
//
// CALLS:
//      - 
// COMMENT: The Eigenvalues of the companion
//          matrix correspond to the negative, 
//          inverse roots of M(z), lambda= - 1/z
//
local mroot,compmat,K,p;

K=rows(AA);
p=cols(AA)/K;

if p > 1;
compmat=AA|(eye(K*(p-1))~zeros(K*(p-1),K));
else;
compmat=AA;
endif;


mroot=maxc(abs(eig(compmat)));

retp(mroot);
endp;







