// varmaAUXprocs.ind contains the following procedures: 

//      varma_simu
//      far_restrictions
//      fma_restrictions
//      dma_restrictions
//      stand_restrictions 
//      getvarma_mat
//      varma_irfun



proc(2)=varma_simu(constant,ar_mat,ma_mat,icovm,T,inival,uerr,ustate);
// {y,yma}=varma_simu(constant,ar_mat,ma_mat,icovm,T,inival,uerr,ustate);
// simulates the following Varma model for T periods:
//
//  y_t = constant + ar_mat *(y'_t-1 ...y'_t-p)'+ ...
//       u_t + ma_mat * (u'_t-1 ...u_t-q')'
//
//  y_t is assumed to be stationary.
//
//  INPUT:
//           constant     K x 1
//           ar_mat       K x pK
//           ma_mat       K x qK
//           T            1 x 1
//           icovm         K x K covariance matrix
//           inival       1 x 1 number of initial values
//           uerr         K x T user-specified shocks
//           ustate             state for the rndKMn procedure
// 
//  OUTPUT:
// 
//           y            K x T  series
//           yma          K x T  series adjusted by 
//                               theoretical mean
//
//  CALLS: 
//           - - - 
//
//  Errors are assumed to be normal.
//  Version: 1.0
//
//  CHECKED for BIVARIATE VAR WITH TWO LAGS
//  


local K,p,q,longrunmat,mean_yt,maxlag,y,yma,erreg,yreg,
      err,newstate;

// dimensions of model

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


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


mean_yt=inv(longrunmat)*constant;

maxlag=maxc(p|q);

if (rows(uerr))==0;

    y=ones(1,T+inival+maxlag).*.mean_yt;
    {err,newstate}=rndKMn(K,T+inival+maxlag,ustate);
    err=chol(icovm)'err;

else;

    inival=0;
    T=cols(uerr);
    y=ones(1,T+inival+maxlag).*.mean_yt;
    err=zeros(K,maxlag)~uerr;    

endif;

for i ((maxlag+1),(T+inival+maxlag),1);
    yreg=rev(y[.,(i-p):(i-1)]')';
    yreg=vec(yreg);
    erreg=rev(err[.,(i-q):(i-1)]')';
    erreg=vec(erreg);
    y[.,i]=constant+ar_mat*yreg+err[.,i]+ma_mat*erreg;
endfor;


y=y[.,(maxlag+inival+1):cols(y)];
yma=y-mean_yt;

retp(y,yma);
endp;

proc(2)=far_restrictions(K,p,q);
// {RR,r}=far_restrictions(K,p,q)
// computes the restriction matrices for 
// a VARMA(p,q) in Final Equations form:
//
//  y_t = alfa_1 y_(t-1)+ ... + ...+u_t + ... M_q u_(t-q) 
//
//  Restriction matrices are constructed as
//  in Lüktepohl 2005, i.e.
//
//  (alfa_0     
//   beta )     = vec [I-A_0 : A : M]
//
//   (alfa_0
//    beta )    = RR gamma + r
//
//   INPUT:
//          K   1 x 1   Dimension of Series
//          p   1 x 1   AR - Order
//          q   1 x 1   MA - Order
//
//   OUTPUT:
//          RR  ((p+1+q)*K^2) x (p + q*K^2) 
//          r   ((p+1+q)*K^2) x 1 
//                      ... as above
//   see documentation!
//
//   CALLS:
//           - - - 
//

local numb_ma_para,numb_totalpara,RR,r;

numb_ma_para=q*(K^2);
numb_totalpara=(p+1)*(K^2)+numb_ma_para;

RR=(eye(p)~zeros(p,numb_ma_para)).*.vec(eye(K));
RR=zeros((K^2),(p+q*(K^2)))|RR|(zeros(numb_ma_para,p)~eye(numb_ma_para));

r=zeros((numb_totalpara),1);

retp(RR,r);
endp;


proc(2)=fma_restrictions(K,p,q);
// {RR,r}=fma_restrictions(K,p,q)
// computes the restriction matrices for 
// a VARMA(p,q) in Final Moving Average form:
//
//  y_t = A_1 y_(t-1)+ ... + u_t - m_1 u_(t-1)- ... - m_q u_(t-q) 
//
//  Restriction matrices are constructed as
//  in Lüktepohl 2005, i.e.
//
//  (alfa_0     
//   beta )     = vec [I-A_0 : A : M]
//
//   (alfa_0
//    beta )    = RR gamma + r
//
//   INPUT:
//          K   1 x 1   Dimension of Series
//          p   1 x 1   AR - Order
//          q   1 x 1   MA - Order
//
//   OUTPUT:
//          RR  ((p+1+q)*K^2) x (p*K^2+q) = # parameters x # free parameters  
//          r   ((p+1+q)*K^2) x 1 
//                      ... as above
//   see documentation!
//
//   CALLS:
//           - - - 
//

local numb_ma_para,numb_totalpara,RR,r;

numb_ma_para=q*(K^2);
numb_totalpara=(p+1)*(K^2)+numb_ma_para;

RR=(zeros(q,p*(K^2))~eye(q)).*.vec(eye(K));
RR=zeros((K^2),(q+p*(K^2)))|(eye(p*(K^2))~zeros(p*(K^2),q))|RR;

r=zeros((numb_totalpara),1);

retp(RR,r);
endp;


proc(4)=EC_restrictions(p);
// {RRc,rc,RR,r}=EC_restrictions(p)
// computes the restriction matrices for 
// a K-dimensional VARMA(p1,p2,..pK) in Echelon form:
//
//  A_0 y_t = A_1 y_(t-1)+ ... + ...+ A_0 u_t + ... M_q u_(t-q) 
//
//  Restriction matrices are constructed as
//  in Lüktepohl 2005, i.e.
//
//  (alfa_0     
//   beta )     = vec [I-A_0 : A : M]
//
//   (alfa_0
//    beta )    = RRc gamma + rc
//
//   and
//
//  (alfa_0     
//   beta )     = vec [A_0 : A : M]
//
//   (alfa_0
//    beta )    = RR gamma + r
//
//
//   INPUT:
//          p   K x 1   Vector of Kronecker Indices
//
//
//   OUTPUT:
//          RRc         Restriction matrices  
//          rc          as defined above.
//          RR  
//          r   
//
//   CALLS:
//           - - - 
//

local KK,maxp,pkimat,k,i,j,AA,MM,AAMM,
      AAp,MMq,numbfreeparam,numbparam,
      lauf,RRc,rc,RR,r;


KK=rows(p);
maxp=maxc(p);

// computation of pki_mat:

pkimat=zeros(KK,KK);

for k (1,KK,1);
    for i (1,KK,1);

        if k==i;
            pkimat[k,i]=p[k]; 
        elseif k > i;
            pkimat[k,i]=minc((p[k]+1)|p[i]);
        else;
            pkimat[k,i]=minc((p[k])|p[i]);  
        endif;

    endfor;
endfor;

// Computation of 
// I-A_0 A_1 A_2 .... A_p M_1 .... M_q
// by putting 0 and 1 s for restricted
// and unrestricted parameters, respectively

// AR - Part
AA={};

for j (0,maxp,1);

    AAp=zeros(KK,KK);
    
    for k (1,KK,1);
        for i (1,KK,1);
        
        // Diagonal Elements:
        if k==i;
            if j==0;
            AAp[k,i]=0;
            elseif (j ne 0) and (j le p[k]);            
            AAp[k,i]=1;
            endif;
        // Off-Diagonals:
        else;
            if (j ge (p[k]-pkimat[k,i]+1)) and (j le p[k]);
                AAp[k,i]=1;
            endif;
        endif;

        endfor;
    endfor;

    AA=AA~AAp;

endfor;


MM={};
// MA - Part
for j (1,maxp,1);
    MMq=zeros(KK,KK);
    for k (1,KK,1);
        for i (1,KK,1);

        if  (j le p[k]);
        MMq[k,i]=1;
        endif;

        endfor;
    endfor;
    MM=MM~MMq;
endfor;

// Setting up [I-A_0 : A : M]:

AAMM=AA~MM;

AAMM=vec(AAMM);

// Making RRc:
numbfreeparam=sumc(AAMM);
numbparam=2*maxp*(KK^2)+KK^2;

lauf=1;
RRc=zeros(numbparam,numbfreeparam);

for i (1,numbparam,1);
    if AAMM[i] ne 0;
    RRc[i,lauf]=1;
    lauf=lauf+1;
    endif;
endfor;

// Making rc:

rc=zeros((numbparam),1);

// Making RR and r:


RR=RRc;
RR[1:(KK^2),.]=-RR[1:(KK^2),.];
r=vec(eye(KK))|zeros((numbparam-(KK^2)),1);


retp(RRc,rc,RR,r);
endp;



proc(2)=dma_restrictions(s,p,q);
//  proc(2)=dma_rest(p,q,s)  
//  computes the restriction matrices for 
//  a K-dimensional VARMA in diagonal ma form:
//
//  A_0 y_t = A_1 y_(t-1)+ ... + ...+ A_0 u_t + ... M_q u_(t-q) 
//
//  Restriction matrices are constructed as
//  in Lüktepohl 2005, i.e.
//
//  (alfa_0     
//   beta )     = vec [I-A_0 : A : M]
//
//   (alfa_0
//    beta )    = RRc gamma + rc
//
//   and
//
//  (alfa_0     
//   beta )     = vec [A_0 : A : M]
//
//   (alfa_0
//    beta )    = RR gamma + r
//
//  INPUTS:
//          p       AR - Order
//          q       MA - Order
//          s       Dimension of Series
//
//  OUTPUTS:
//          RR   Restriction Matricees 
//          r   
//
//  CALLS:
//          none
//
    local RR,r,J,numP,mmat,i;


//number of unrestricted parameters:
    J=(s^2)*p+s*q;

// total number of parameters:
    numP=(s^2)*(p+q+1);


    RR=zeros(s^2,J);

    RR=RR|(eye((s^2)*p)~zeros(((s^2)*p),J-((s^2)*p)));

    mmat=zeros(s^2,s);

    for i (1,s,1); 
        mmat[(i-1)*(s+1)+1,i]=1; 
    endfor; 


    mmat=eye(q).*.mmat; 


    RR=RR|(zeros((s^2)*q,(s^2)*p)~mmat); 

    r=zeros(numP,1);

    retp(RR,r);
endp;



proc(2)=stand_restrictions(K,p,q);
// {RR,r}=stand_restrictions(K,p,q)
// computes the restriction matrices for 
// a VARMA(p,q) in standard form:
//
//  y_t = A_1 y_(t-1)+ ... + ...+u_t + ... M_q u_(t-q) 
//
//  Restriction matrices are constructed as
//  in Lüktepohl 2005, i.e.
//
//  (alfa_0     
//   beta )     = vec [I-A_0 : A : M]
//
//   (alfa_0
//    beta )    = RR gamma + r
//
//   INPUT:
//          K   1 x 1   Dimension of Series
//          p   1 x 1   AR - Order
//          q   1 x 1   MA - Order
//
//   OUTPUT:
//          RR  ((p+1+q)*K^2) x (p*K^2 + q*K^2) 
//          r   ((p+1+q)*K^2) x 1 
//                      ... as above
//   see documentation!
//
//   CALLS:
//           - - - 
//

local numb_freepara,numb_totalpara,RR,r;

numb_totalpara=(p+1+q)*(K^2);
numb_freepara=(p+q)*(K^2);

RR=zeros(K^2,numb_freepara);
RR=RR|(eye(numb_freepara));

r=zeros((numb_totalpara),1);

retp(RR,r);
endp;




proc(3)=getvarma_mat(gammma,RR,r,K,p,q);
// {AA0,AA,MM}=getvarma_mat(gammma,RR,r,K,p,q) transforms
// the parameter vector gamma in the equation
//
// beta=vec [I-A_0 : A : M]=RR gamma + r  
//
// to the matrices A_0, A and M. 
//
// INPUTS: 
//          gamma  (n1 x 1)     Matrix of free parameters
//          RR     (n2 x n1)    Restriction matrices
//          r      (n2 x n1)    
//          K       1 x 1       Dimension of the series
//          p       1 x 1       AR - order
//          q       1 x 1       MA - order
//
//  OUTPUTS:
//
//          AA0     K x K       
//          AA      K x pK      
//          MM      K x qK      
//
//  CALLS:
//          reshaper

local beta,nAR1,nAR2,AA0,AA,MM;

beta=RR*gammma+r;

nAR1=((K^2)+1);
nAR2=nAR1+p*(K^2)-1;


AA0=beta[1:(K^2),.];

AA0=-reshaper(AA0,K,K)+eye(K);



AA=beta[nAR1:nAR2,.];
AA=reshaper(AA,K,p*(K));

MM=beta[nAR2+1:rows(beta),.];

MM=reshaper(MM,K,q*(K));

retp(AA0,AA,MM);
endp;




proc(1)=varma_irfun(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 
// CALLS:
//          -   
// See Luetkepohl (2005) pp. 426 for Documentation.
//
// NOT CHECKED
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~vec(JJ*AAAJ*HH);
    AAAJ=AAAJ*AAA;

endfor;


retp(irfun);
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;
