// CCAprocs.ind
//
//
//      CCAprocs contains the following procedures: 
//
//
//      ss_CCAr
//      ss_CCAm
//
//
//
//
//
//
//
//
//



/*-----------------------------------------------------------------------------------------------------*/
proc(8)=ss_CCAr(y,outputt,graffs,f0,p0,n0);
//
//      Estimates a state space model in pred. error form
//      using the CCA subspace algorithm.
//      ss_CCAr uses regression to obtain the system matrices' estimates. 
//
//  
//      MODEL: 
//
//              x_t+1 = A x_t + B e_t
//              y_t   = C x_t + e_t
//
//      USAGE: 
//
//      {a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat}...
//      =ss_CCAr(y,outputt,graffs,f0,p0,n0);      
//
//
//      INPUTS: 
//
//          y               (s x T)         series 
//          outputt         (1 x 1)         1 - print output 0 - no print
//          graffs          (1 x 1)         1 - show graphs 0 - no graphs
//          f0              (1 x 1)         aux. parameter, if f0==0, f is 
//                                          estimated o/w f0 is used
//          p0              (1 x 1)         aux. parameter, if p0==0, f is 
//                                          estimated o/w p0 is used
//          n0              (1 x 1)         state dimension, if n0==0, n is 
//                                          estimated o/w n0 is used
//                                          
//      OUTPUTS: 
//
//          a_hut           (n x n)         estimate of A
//          b_hut           (n x s)         estimate of B
//          c_hut           (s x n)         estimate of C
//          e_hut           ()              estimate of residuals
//          x_ss            ()              estimate of the state vector
//          f               (1 x 1)         used f
//          p               (1 x 1)         used p
//          n_hat           (1 x 1)         used n
//
//
//
//
//      COMMENTS: 
//
//      y is assumed to have mean zero. 
//
//      if f or p are estimated, they are set to 2 * p_AIC, where
//      p_AIC is the lag-length of a VAR chosen by AIC
//
//      if the state dimension, n, is estimated, the criteria prop. by
//      .... () is used. 
//
//      checked
//
    local t,s,mean_y,y_tilde,yy,T_fp,Y_minus,Y_plus,
        xx,HH,beta_fp,Gamma_plus,Gamma_minus,W_f,W_p,
        svd_mat,U,V,SIG,n1_hat,n2_hat,UN,SIGN,VN,  
        SIGNROOT,K_P,x_ss,y_ss,ixx,c_hut,a_hut,b_hut,e_hut,
        xax,deltaa;
    local f,p,n_hat,XandE,AandB;
    local betar,betac;

    t=cols(y);
    s=rows(y); 
    
    // STEP 1 
    y_tilde=y;

    if (f0==0) or (p0==0); 
        {f,p}=getfp(y,1,round(sqrt(t)));
        f=maxc(f|n0+1);
        p=maxc(p|n0+1);
    else;
        f=f0;
        p=p0;
    endif;
    T_fp=t-f-p+1;


    // STEP 2 
    Y_minus=rev(y_tilde[.,1:(1+p-1)]');
    Y_minus=vec(Y_minus');
    for i (2,T_fp,1);
        xx=rev(y_tilde[.,i:(i+p-1)]');
        xx=vec(xx');
        Y_minus=Y_minus~xx;
    endfor;

    yy=rev(y_tilde')';

    Y_plus=rev(yy[.,1:(1+f-1)]');
    Y_plus=vec(Y_plus');
    for i (2,T_fp,1);
        xx=rev(yy[.,i:(i+f-1)]');
        xx=vec(xx');
        Y_plus=Y_plus~xx;
    endfor;
    Y_plus=rev(Y_plus')';

    HH=(Y_plus*Y_minus')/T_fp;

    //beta_fp=HH*invpd(moment(Y_minus',0)/T_fp); 
            
    beta_fp=vec(HH)/((moment(Y_minus',0)/T_fp)'.*.eye(s*f)); 
    beta_fp=reshaper(beta_fp,rows(HH),cols(HH));
        
    // STEP 3 

    Gamma_plus=moment(Y_plus',0)/T_fp;
    Gamma_minus=moment(Y_minus',0)/T_fp;

    W_f=inv(chol(Gamma_plus)');
    W_p=chol(Gamma_minus)';
    svd_mat=W_f*beta_fp*W_p;

    {U,SIG,V}=svdusv(svd_mat);

    if n0==0;
        {n_hat,n1_hat,n2_hat}=ss_dim(SIG,s,T_fp);
    else;
        n_hat=n0;
    endif;


    // STEP 4 
    UN=U[.,1:n_hat];
    SIGN=SIG[1:n_hat,1:n_hat];
    VN=V[.,1:n_hat];

    SIGNROOT=chol(SIGN)';

    K_p=SIGNROOT*VN'inv(W_p);


    // STEP 5 

    x_ss=K_p*Y_minus;

    y_ss=y_tilde[.,p+1:T-f+1];


    // SYSTEM MATRICES      

    c_hut=(y_ss')/(x_ss');
    c_hut=c_hut';

    e_hut=y_ss-c_hut*x_ss;
    
    XandE=(x_ss[.,1:cols(x_ss)-1])|e_hut[.,1:(cols(x_ss)-1)];
    
    AandB=pinv(XandE*XandE')*(XandE*x_ss[.,2:cols(x_ss)]');
    AandB=AandB';
    
    a_hut=AandB[.,1:n_hat];
    b_hut=AandB[.,n_hat+1:cols(AandB)];
    


    if outputt==1;

        print "*******************************************";
        print "Estimated State Space Model";    
        print "*******************************************";
        print "Observations :" t;
        print "Dimension y  :" s;
        print "Effective Obs:" t_fp;
        print " ";
        print "Aux. Parameters";
        print "f" f;
        print "p" p;
        print "n" n_hat;
        print " ";
        print "Estimated System Matrices:";
        print "A" a_hut;
        print "B" b_hut;
        print "C" c_hut;
        print " ";
        print "Estimated Covariance Mat.:";
        print moment(e_hut',0)/t_fp;

    endif;


    if graffs==1;
        graphset;
        _pdate="";
        _ptitle="Singular Values";
        _plctrl={-1};
        _pstype={8};
        xlabel("State Dimension");
        xax=seqa(1,1,rows(SIG));
        xy(xax,diag(SIG));
    endif;

    retp(a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat);
endp;

/*-----------------------------------------------------------------------------------------------------*/
proc(8)=ss_CCArYW(y,outputt,graffs,f0,p0,n0);
//
//      Estimates a state space model in pred. error form
//      using the CCA subspace algorithm.
//      ss_CCAr uses regression to obtain the system matrices' estimates. 
//
//  
//      MODEL: 
//
//              x_t+1 = A x_t + B e_t
//              y_t   = C x_t + e_t
//
//      USAGE: 
//
//      {a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat}...
//      =ss_CCAr(y,outputt,graffs,f0,p0,n0);      
//
//
//      INPUTS: 
//
//          y               (s x T)         series 
//          outputt         (1 x 1)         1 - print output 0 - no print
//          graffs          (1 x 1)         1 - show graphs 0 - no graphs
//          f0              (1 x 1)         aux. parameter, if f0==0, f is 
//                                          estimated o/w f0 is used
//          p0              (1 x 1)         aux. parameter, if p0==0, f is 
//                                          estimated o/w p0 is used
//          n0              (1 x 1)         state dimension, if n0==0, n is 
//                                          estimated o/w n0 is used
//                                          
//      OUTPUTS: 
//
//          a_hut           (n x n)         estimate of A
//          b_hut           (n x s)         estimate of B
//          c_hut           (s x n)         estimate of C
//          e_hut           ()              estimate of residuals
//          x_ss            ()              estimate of the state vector
//          f               (1 x 1)         used f
//          p               (1 x 1)         used p
//          n_hat           (1 x 1)         used n
//
//
//
//
//      COMMENTS: 
//
//      y is assumed to have mean zero. 
//
//      if f or p are estimated, they are set to 2 * p_AIC, where
//      p_AIC is the lag-length of a VAR chosen by AIC
//
//      if the state dimension, n, is estimated, the criteria prop. by
//      .... () is used. 
//
//      checked
//
    local t,s,mean_y,y_tilde,yy,T_fp,Y_minus,Y_plus,
        xx,HH,beta_fp,Gamma_plus,Gamma_minus,W_f,W_p,
        svd_mat,U,V,SIG,n1_hat,n2_hat,UN,SIGN,VN,  
        SIGNROOT,K_P,x_ss,y_ss,ixx,c_hut,a_hut,b_hut,e_hut,
        xax,deltaa;
    local f,p,n_hat,XandE,AandB;
    local betar,betac;
    local ae,be,ce; 

    t=cols(y);
    s=rows(y); 
    
    // STEP 1 
    y_tilde=y;

    if (f0==0) or (p0==0); 
        {f,p}=getfp(y,1,round(sqrt(t)));
        f=maxc(f|n0+1);
        p=maxc(p|n0+1);
    else;
        f=f0;
        p=p0;
    endif;
    T_fp=t-f-p+1;


    // STEP 2 
    Y_minus=rev(y_tilde[.,1:(1+p-1)]');
    Y_minus=vec(Y_minus');
    for i (2,T_fp,1);
        xx=rev(y_tilde[.,i:(i+p-1)]');
        xx=vec(xx');
        Y_minus=Y_minus~xx;
    endfor;

    yy=rev(y_tilde')';

    Y_plus=rev(yy[.,1:(1+f-1)]');
    Y_plus=vec(Y_plus');
    for i (2,T_fp,1);
        xx=rev(yy[.,i:(i+f-1)]');
        xx=vec(xx');
        Y_plus=Y_plus~xx;
    endfor;
    Y_plus=rev(Y_plus')';

    HH=(Y_plus*Y_minus')/T_fp;

    //beta_fp=HH*inv(moment(Y_minus',0)/T_fp); 
            
    beta_fp=vec(HH)/((moment(Y_minus',0)/T_fp)'.*.eye(s*f)); 
    beta_fp=reshaper(beta_fp,rows(HH),cols(HH));
        
    // STEP 3 

    Gamma_plus=moment(Y_plus',0)/T_fp;
    Gamma_minus=moment(Y_minus',0)/T_fp;

    W_f=inv(chol(Gamma_plus)');
    W_p=chol(Gamma_minus)';
    svd_mat=W_f*beta_fp*W_p;

    {U,SIG,V}=svdusv(svd_mat);

    if n0==0;
        {n_hat,n1_hat,n2_hat}=ss_dim(SIG,s,T_fp);
    else;
        n_hat=n0;
    endif;


    // STEP 4 
    UN=U[.,1:n_hat];
    SIGN=SIG[1:n_hat,1:n_hat];
    VN=V[.,1:n_hat];

    SIGNROOT=chol(SIGN)';

    K_p=SIGNROOT*VN'inv(W_p);


    // STEP 5 

    x_ss=K_p*Y_minus;

    y_ss=y_tilde[.,p+1:T-f+1];


    // SYSTEM MATRICES      

    c_hut=(y_ss')/(x_ss');
    c_hut=c_hut';

    e_hut=y_ss-c_hut*x_ss;
    
    XandE=(x_ss[.,1:cols(x_ss)-1])|e_hut[.,1:(cols(x_ss)-1)];
    

    
    {a_hut,be,ce}=var_yw(x_ss,1);
  
    b_hut=(x_ss[.,2:cols(x_ss)]')/(e_hut[.,1:(cols(x_ss)-1)]');
    b_hut=b_hut';
  


    if outputt==1;

        print "*******************************************";
        print "Estimated State Space Model";    
        print "*******************************************";
        print "Observations :" t;
        print "Dimension y  :" s;
        print "Effective Obs:" t_fp;
        print " ";
        print "Aux. Parameters";
        print "f" f;
        print "p" p;
        print "n" n_hat;
        print " ";
        print "Estimated System Matrices:";
        print "A" a_hut;
        print "B" b_hut;
        print "C" c_hut;
        print " ";
        print "Estimated Covariance Mat.:";
        print moment(e_hut',0)/t_fp;

    endif;


    if graffs==1;
        graphset;
        _pdate="";
        _ptitle="Singular Values";
        _plctrl={-1};
        _pstype={8};
        xlabel("State Dimension");
        xax=seqa(1,1,rows(SIG));
        xy(xax,diag(SIG));
    endif;

    retp(a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat);
endp;


/*-----------------------------------------------------------------------------------------------------*/
proc(8)=ss_CCAm(y,outputt,graffs,f0,p0,n0);
//
//      Estimates a state space model in pred. error form
//      using the CCA subspace algorithm.
//      ss_CCAm uses covar. matrices to obtain the system matrices' estimates. 
//
//  
//      MODEL: 
//
//              x_t+1 = A x_t + B e_t
//              y_t   = C x_t + e_t
//
//      USAGE: 

//      {a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat}...
//      =ss_CCAm(y,outputt,graffs,f0,p0,n0);      
//
//
//      INPUTS: 
//
//          y               (s x T)         series 
//          outputt         (1 x 1)         1 - print output 0 - no print
//          graffs          (1 x 1)         1 - show graphs 0 - no graphs
//          f0              (1 x 1)         aux. parameter, if f0==0, f is 
//                                          estimated o/w f0 is used
//          p0              (1 x 1)         aux. parameter, if p0==0, f is 
//                                          estimated o/w p0 is used
//          n0              (1 x 1)         state dimension, if n0==0, n is 
//                                          estimated o/w n0 is used
//                                          
//      OUTPUTS: 
//
//          a_hut           (n x n)         estimate of A
//          b_hut           (n x s)         estimate of B
//          c_hut           (s x n)         estimate of C
//          e_hut           (s x T_fp)      estimate of residuals
//                                          from p+1 ... T-f
//          x_ss            (n x T_fp)      estimate of the state vector
//                                          from p+1 ... T-f
//          f               (1 x 1)         used f
//          p               (1 x 1)         used p
//          n_hat           (1 x 1)         used n
//
//
//      COMMENTS: 
//
//      y is assumed to have mean zero. 
//
//      if f or p are estimated, they are set to 2 * p_AIC, where
//      p_AIC is the lag-length of a VAR chosen by AIC
//
//      if the state dimension, n, is estimated, the criteria prop. by
//      .... () is used. 
//
//      checked!
//
    local t,s,mean_y,y_tilde,yy,T_fp,Y_minus,Y_plus,
        xx,HH,beta_fp,Gamma_plus,Gamma_minus,W_f,W_p,
        svd_mat,U,V,SIG,n1_hat,n2_hat,UN,SIGN,VN,  
        SIGNROOT,K_P,x_ss,y_ss,ixx,c_hut,a_hut,b_hut,e_hut,
        xax,deltaa;
    local f,p,n_hat;
    local betar,betac;

    t=cols(y);
    s=rows(y); 
    
    // STEP 1 
    y_tilde=y;

    if (f0==0) or (p0==0); 
        {f,p}=getfp(y,1,round(sqrt(t)));
        f=maxc(f|n0+1);
        p=maxc(p|n0+1);
    else;
        f=f0;
        p=p0;
    endif;
    T_fp=t-f-p+1;


    // STEP 2 
    Y_minus=rev(y_tilde[.,1:(1+p-1)]');
    Y_minus=vec(Y_minus');
    for i (2,T_fp,1);
        xx=rev(y_tilde[.,i:(i+p-1)]');
        xx=vec(xx');
        Y_minus=Y_minus~xx;
    endfor;

    yy=rev(y_tilde')';

    Y_plus=rev(yy[.,1:(1+f-1)]');
    Y_plus=vec(Y_plus');
    for i (2,T_fp,1);
        xx=rev(yy[.,i:(i+f-1)]');
        xx=vec(xx');
        Y_plus=Y_plus~xx;
    endfor;
    Y_plus=rev(Y_plus')';

    HH=(Y_plus*Y_minus')/T_fp;

    beta_fp=HH*invpd(moment(Y_minus',0)/T_fp);
    
    //beta_fp=vec(HH)/((moment(Y_minus',0)/T_fp)'.*.eye(s*f)); 
    //beta_fp=reshaper(beta_fp,rows(HH),cols(HH));
   
   
    // STEP 3 
    Gamma_plus=moment(Y_plus',0)/T_fp; 
    Gamma_minus=moment(Y_minus',0)/T_fp;
     
 
    W_f=inv(chol(Gamma_plus)');
    W_p=chol(Gamma_minus)';
    svd_mat=W_f*beta_fp*W_p;

    {U,SIG,V}=svdusv(svd_mat);


    if n0==0;
        {n_hat,n1_hat,n2_hat}=ss_dim(SIG,s,T);
    else;
        n_hat=n0;
    endif;
    

    // STEP 4 
    UN=U[.,1:n_hat];
    SIGN=SIG[1:n_hat,1:n_hat];
    VN=V[.,1:n_hat];

    SIGNROOT=chol(SIGN)';

    K_p=SIGNROOT*VN'inv(W_p);


    // STEP 5 

    x_ss=K_p*Y_minus;

    y_ss=y_tilde[.,p+1:T-f+1];


    // SYSTEM MATRICES      
    xx=K_p*Gamma_minus*(K_p'); 
    ixx=inv(xx);

    c_hut=HH[1:s,.]*K_p'ixx;

    e_hut=y_ss-c_hut*x_ss;


    a_hut=K_p*(HH[1:s,.]|Gamma_minus[1:(p-1)*s,.])*K_p'ixx;

    deltaa=Gamma_plus[1:s,1:s]-c_hut*(K_p*Gamma_minus*K_p')*c_hut';

    b_hut=(K_p*Gamma_minus[.,1:s]-a_hut*K_p*(HH[1:s,.]'))'/deltaa;
    b_hut=b_hut';


    if outputt==1;

        print "*******************************************";
        print "Estimated State Space Model";    
        print "*******************************************";
        print "Observations :" t;
        print "Dimension y  :" s;
        print "Effective Obs:" t_fp;
        print " ";
        print "Aux. Parameters";
        print "f" f;
        print "p" p;
        print "n" n_hat;
        print " ";
        print "Estimated System Matrices:";
        print "A" a_hut;
        print "B" b_hut;
        print "C" c_hut;
        print " ";
        print "Estimated Covariance Mat.:";
        print moment(e_hut',0)/t_fp;

    endif;


    if graffs==1;
        graphset;
        _pdate="";
        _ptitle="Singular Values";
        _plctrl={-1};
        _pstype={8};
        xlabel("State Dimension");
        xax=seqa(1,1,rows(SIG));
        xy(xax,diag(SIG));
    endif;

    retp(a_hut,b_hut,c_hut,e_hut,x_ss,f,p,n_hat);
endp;

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


proc(1)=ss_frc(y,statev,res,A,B,C,f,p,horizon);
    // computes forecasts from a SS model estimated
    // by CCA. 
    // checked 
    local sspred,x_old,x_new,L,ABC,t0,KK;
    local Aj,yp,e_new,e_old,TT;

    TT=cols(y);
    KK=rows(C);
    ABC=(A-B*C);
     
    t0=TT-f+1;
    
    x_old=statev[.,cols(statev)];
   
    format 8,8; 
    for j (1,f,1);
        
        x_new=ABC*x_old+B*y[.,t0+j-1];
        x_old=x_new; 
         
    endfor;

    sspred={};
    Aj=eye(cols(x_old));

    for j (1,horizon,1);

        yp=C*Aj*x_old;
        Aj=Aj*A;  
        sspred=sspred~yp;

    endfor;

retp(sspred);
endp;


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

proc(3)=ss_dim(SIG,dim_y,T_eff);
// seems allright!; 
    local HT,sig_vec,svc_vec,ba_vec,n_hat,n1_hat,n2_hat,xax;

    HT=ln(T_eff);
    sig_vec=diag(SIG);

    
    svc_vec=sig_vec[2:rows(SIG)].^2+seqa(1,1,rows(SIG)-1)*2*dim_y*HT/T_eff;
 

    ba_vec=-ln(ones(rows(SIG)-1,1)-sig_vec[1:(rows(SIG)-1)].^2)+seqa(1,1,rows(SIG)-1)*2*dim_y*HT/T_eff;



    n1_hat=minindc(svc_vec);
    n2_hat=minindc(ba_vec);

    n_hat=n2_hat;


retp(n_hat,n1_hat,n2_hat);
endp;


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


proc(2)=getfp(y,const,pmax_AIC);
    // gets auxialiry parameters f,p in 
    // the CCA subspace algorithm
    // ok 
    local f,p,method;
    local aic,aic_p,hq,hq_p,sc,sc_p; 
    method=1;
 
    {aic,aic_p,hq,hq_p,sc,sc_p}=ord_VAR(y,const,1,pmax_aic,method);

    f=2*aic_p;
    p=2*aic_p;

retp(f,p);
endp;


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;




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

proc(5)=ss_BQ(A,C,K,covmat,horizon,impuls,scalemat);
// {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  = ss_BQ(A,C,K,covmat,horizon,impuls);
// computes the structural impulse responses of an estimated
// state-space model.
//
//
//      OUTPUT: 
//
//
//          y_ir        impulse responses
//          y_ira       accumulated impulse responses
//
// checked
local longrunmat,longrunmat_struc,cont_impact,
      y_ir,y_ira,impuls_struc;


longrunmat=eye(2)+c*inv(eye(cols(a))-a)*k;



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

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

if rows(scalemat) ne 0;
cont_impact=inv(scalemat)*cont_impact;
endif;


impuls_struc=cont_impact*impuls;

y_ir=ss_IR(A,K,C,horizon+1,impuls_struc);

y_ira=zeros(2,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)=ss_BQ2(A,C,K,covmat,horizon,impuls,scalemat);
// {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  = ss_BQ(A,C,K,covmat,horizon,impuls);
// computes the structural impulse responses of an estimated
// state-space model.
//
// checked
local longrunmat,longrunmat_struc,cont_impact,
      y_ir,y_ira,impuls_struc;


longrunmat=eye(2)+c*inv(eye(cols(a))-a)*k;


cont_impact=getD0(longrunmat,covmat);

longrunmat_struc=longrunmat*cont_impact; 


if rows(scalemat) ne 0;
cont_impact=inv(scalemat)*cont_impact;
endif;


impuls_struc=cont_impact*impuls;

y_ir=ss_IR(A,K,C,horizon+1,impuls_struc);

y_ira=zeros(2,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)=ss_BQ3(A,C,K,covmat,horizon,impuls,scalemat);
// {longrunmat,longrunmat_struc,cont_impact,y_ir,y_ira}
//  = ss_BQ(A,C,K,covmat,horizon,impuls);
// computes the structural impulse responses of an estimated
// state-space model.
//
// checked
local longrunmat,longrunmat_struc,cont_impact,
      y_ir,y_ira,impuls_struc,A_1,iA1;


A_1=(eye(cols(a))-a);
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=eye(2)+c*iA1*k;


cont_impact=getD04(longrunmat,covmat);

longrunmat_struc=longrunmat*cont_impact; 


if rows(scalemat) ne 0;
cont_impact=inv(scalemat)*cont_impact;
endif;


impuls_struc=cont_impact*impuls;

y_ir=ss_IR(A,K,C,horizon+1,impuls_struc);

y_ira=zeros(2,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)=ss_IR(A,B,C,horizon,impuls);
// checked 
    local y_ir,A_old,KK;

    KK=rows(C);
    y_ir=zeros(KK,horizon);

    y_ir[.,1]=impuls;

    A_old=eye(rows(A));

    for i (2,horizon,1);
        y_ir[.,i]=C*A_old*B*impuls;
        A_old=A_old*A;
    endfor;

    retp(y_ir);
endp;

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

proc(1)=ss_IRmats(A,B,C,horizon);
// not checked 
    local y_irmat,A_old,KK,y_ir;

    KK=rows(C);

    A_old=eye(rows(A));

    y_irmat=vec(eye(KK));


    for i (2,horizon,1);
        y_ir=C*A_old*B;
        A_old=A_old*A;
        y_irmat=y_irmat~(vec(y_ir));
    endfor;

    retp(y_irmat);
endp;

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

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

    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); 

    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; 

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

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(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;

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