// VAR estimation methods, programmed according to Luetkepohl 2005,
// New Introduction to Multivariate Time Series Analysis


proc(6)=ord_VAR(y,mud,minp,maxp,method);
/*
    Lag selection in a VAR with with/without intercept
    computes the AIC,HQ and Schwarz/BIC criterion. 

    USAGE:
    {aic,aic_p,hq,hq_p,sc,sc_p}=ord_VAR(y,mud,maxp,method);
    
    INPUTS
    y           k x T       series
    mud         1 x 1       dummy: 0 - no intercept, 1 - intercept
    minp        1 x 1       minimal lag order
    maxp        1 x 1       maximal lag order
    method      1 x 1       estimation method dummy
                            1 - OLS
                            2 - YW
                            3 - OLS on demeaned data

    OUTPUT
    aic         maxp+1 x 2  aic for lags 0 - maxp
    aic_p       1 x 1       optimal lag length acc. to aic      
    hq          maxp+1 x 2  hq for lags 0 - maxp          
    hq_p        1 x 1       optimal lag length acc. to hq
    sc          maxp+1 x 2  sc for lags 0 - maxp
    sc_p        1 x 1       optimal lag length acc. to sc

    

*/
local aic,hq,sc,ae,covmat;
local detcovm,KK,u,TT;
local p,aic_p,hq_p,sc_p;

KK=rows(y);

TT=cols(y);


aic=zeros(maxp+1,1);
hq=zeros(maxp+1,1);
sc=zeros(maxp+1,1);


for p (0,maxp,1);
    
    if method==1; 
        {ae,covmat,u}=var_ols(y[.,(maxp+1-p):cols(y)],mud,p);
    elseif method==2;
        {ae,covmat,u}=var_yw(y[.,(maxp+1-p):cols(y)],p);
    elseif method==3;
        {ae,covmat,u}=var_dm(y[.,(maxp+1-p):cols(y)],p);
    endif; 
        
    detcovm=det((u*u')/(TT-maxp));
    
    AIC[p+1]=ln(detCovM)+2*(p*(KK^2))/(TT-maxp);
    HQ[p+1]=ln(detCovM)+2*ln(ln(TT-maxp))*(p*(KK^2))/(TT-maxp);
    SC[p+1]=ln(detCovM)+ln(TT-maxp)*(p*(KK^2))/(TT-maxp);

endfor;


aic=seqa(0,1,maxp+1)~aic;
hq=seqa(0,1,maxp+1)~hq;
sc=seqa(0,1,maxp+1)~sc;

aic=aic[minp+1:rows(aic),.]; 
hq=hq[minp+1:rows(hq),.]; 
sc=sc[minp+1:rows(sc),.]; 

aic_p=aic[minindc(aic[.,2]),1];
hq_p=hq[minindc(hq[.,2]),1];
sc_p=sc[minindc(sc[.,2]),1];



retp(aic,aic_p,hq,hq_p,sc,sc_p);
endp;



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



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;


