/*********************** PROC JOHANTST *****************************
AUTHOR: Avinash S. Bhati

        -- The author makes no performance guarantees.
        -- For public non-commercial use only.

LAST REVISED: 19 December, 1997.

NOTES:  PROC JOHANTEST tests for the cointegration rank for a VAR
        process, estimates the TRACE and LMAX stats, the eigen
        values, and the eigenvectors. Assuming r = 1, it computes
        the long-run equilibrium coefficients (BETA), the adjustment
        coefficients (ALPHA), the long-run (PI) coefficient, the
        covariance matrix of the errors (OMEGA), the R-squares
        for each of the equations in the VECM. In addition,
        assuming r=1, it will also test for "a" linear restriction
        on the long-run equilibrium coefficients (BETAS). The program
        is based on the methodology detailed in:

        Johansen, Soren, "Cointegration and Hypothesis Testing of
           Cointegration Vectors in Gaussian Vector Autoregressive
           Models", Econometrica, Vol.59, No.6 (Nov 1991) 1551-1580.

        More specifically, the testing of cointegrating rank is
        explained and derived on pp. 1551-1556, and the testing
        of linear restrictions on the BETA matrix (assuming r=1) is
        defined as Corollary 5.3 on pp. 1564.

        The method developed in Johansen (1991) is general in terms
        of the deterministic terms. In the proc below, the costant
        term in the VECM is constrained so that the Cointegrating
        relationship has a constant term.

FORMAT: {alm,atr,aeval,abet,aalp,alrpi,aomeg,arest3,arestp3,arsqs}
         = johantst(data,lg,H);

INPUT:  data = T x P matrix of input variables in levels
        lj   = maximum lag length in the underlying VAR model
        H    = (P+1) column vector for a linear constraint on the
               BETA matrix. Restrictions are of the form H'B=0.
               The H matrix can ONLY be a vector. i.e. joint
               hypothesis CANNOT be tested within the PROC.

            Example 1: H={1,0,0,0,0}
               Tests the hypothesis that b_1=0.
            Example 2: H={1,1,0,0,0}
               Tests the hypothesis that b_1+b_2=0.

               note: if no restrictions are to be tested, define
               H = (P+1) vector of zeros.

OUTPUT: alm    = P+1 x 1 matrix of MAXIMAL EIGENVALUE statistics
                 (last element = 0 because of the constant)
        alm2   = Small sample version of alm
        atr    = P+1 x 1 matrix of TRACE statistics
                 (last element = 0 because of the constant)
        atr2   = Small sample version of atr
        aeval  = P+1 x 1 matrix of EIGEN VALUES
                 (last element = 0 because of the constant)
        abet   = P+1 x 1 matrix of un-normalized long run
                 BETA coefficients. The last element is the constant
                 and the order of the remaining coefficients is the
                 same as that of the incoming data matrix.
        aalp   = P x 1 matrix of adjustment coefficients. The order
                 of the coefficients is the same as that of the
                 incoming data matrix.
        alrpi  = P x P+1 matrix of the PI matrix. PI=ALPHA*BETA'
        aomeg  = P x P covariance matrix of the error terms in the VECM
        arest  = scalar. The Chi-Square Statstic for the restriction
                 defined by the H vector. If the H matrix is a vector
                 of zeros, then this is a missing value (".").
        arestp = scalar. The marginal significance level of the Chi-
                 Square Statistic (arest) defined above. Again, if there
                 is no restriction, this is a missing value (".").
        arsqs  = P x 1 matrix of the R-Squares for the various equations
                 in the VECM. The order is the same as that of the incoming
                 data matrix.

GLOBALS: None

EXTERNAL PROCS: VARLAGS (provided below)
*******************************************************************/
proc(12)=johantst(levels,k,h);
local x,dx,Z0T,Z1T,Z2T,T,R0T,R1T,S10,S01,S00,S11,S00I,S11I,
      G,A,EVALS,EV,EVECS,BETA,LMAX,LMAX2,TRAC,TRAC2,ALPHA,LRPI,OMEGA,
      L_1,V_H,A2,RESTR,RESTR_P,EE,RSS,TSS,RSQ;
x=levels;
dx=trimr(x-lagn(x,1),1,0);
{Z0T,Z2T}=VARLAGS(dx,k-1);
T=rows(Z0T);
/******************************************************************/
/* The Z, R, and S notations below follow identically from        */
/* Johansen (1991)                                                */
/******************************************************************/
Z1T=trimr(lagn(x,k),k,0)~ones(rows(x)-k,1);
R0T=Z0T-Z2T*(Z0T/Z2T);
R1T=Z1T-Z2T*(Z1T/Z2T);
S01=(R0T'R1T)/T;
S10=S01';
S00=(R0T'R0T)/T ; S00I=INVPD(S00);
S11=(R1T'R1T)/T ; S11I=INVPD(S11);
G=inv(chol(S11)');                      /* Cholesky Decomposition */
A=G*S10*S00I*S01*G';             /* Matrix for Eigenvalue Problem */
{EVALS,EVECS}=eigrs2(A);    /* Extracting Eigen- values & vectors */
EV=(EVALS~EVECS);
EVALS=EV[.,1];
LMAX=-T*ln(1-EVALS);                            /* Lmax statistic */
LMAX2=-(T-cols(x)*k)*ln(1-EVALS);            /* Small Sample Lmax */
TRAC=-T*cumsumc(ln(1-EVALS));                  /* Trace Statistic */
TRAC2=-(T-cols(x)*k)*cumsumc(ln(1-EVALS));  /* Small Sample Trace */
EVECS=rev((G'*ev[.,2:cols(EV)])')';
BETA=EVECS[.,1];          /* Single Long Run Cointegrating Vector */
ALPHA=S01*BETA*inv(BETA'*S11*BETA);    /* Adjustment Coefficients */
LRPI=ALPHA*BETA';       /* Long Run Pi matrix: (alph)(beta)'=(Pi) */
EE=Z0T - (Z2T~Z1T)*(Z0T/(Z2T~Z1T));    /* Residuals from the VECM */
RSS=diag((EE-meanc(EE)')'(EE-meanc(EE)')); /* Res. Sum of Squares */
TSS=diag((Z0T-meanc(Z0T)')'(Z0t-meanc(Z0T)'));        /* Total SS */
RSQ=1-RSS./TSS;     /* R squares for each of the equation in VECM */
OMEGA=S00-ALPHA*BETA'*S10; /* COVARIANCE matrix of errors of VECM */
L_1=EVALS[rows(EVALS)];
V_H=EVECS[.,2:cols(EVECS)];
IF sumc(sumc(abs(h))) .== 0;
   RESTR=MISS(0,0);
   RESTR_P=MISS(0,0);
ELSE;
   A2=(inv((1/L_1-1)*(h'*V_H*V_H'*h)));
   RESTR=T*((h'*BETA).^2)*A2;          /* LR test of restrictions */
   RESTR_P=CDFCHIC(RESTR,1);       /* Marginal Significance level */
ENDIF;
RETP(rev(LMAX),rev(LMAX2),rev(TRAC),rev(TRAC2),rev(EVALS),BETA,ALPHA,
     LRPI,OMEGA,RESTR,RESTR_P,RSQ);
endp;


/**********************  PROC VARLAGS  *****************************
**   author: Alan G. Isaac
**   last update: 5 Dec 95      previous: 15 June 94
**   FORMAT
**          { x,xlags } = varlags(var,lags)
**   INPUT
**        var  - T x K matrix
**        lags - scalar, number of lags of var (a positive integer)
**   OUTPUT
**          x -     (T - lags) x K matrix, the last T-lags rows of var
**          xlags - (T - lags) x lags*cols(var) matrix,
**                  being the 1st through lags-th
**                  values of var corresponding to the values in x
**                  i.e, the appropriate rows of x(-1)~x(-2)~etc.
**   GLOBAL VARIABLES: none
**********************************************************************/
proc(2)=varlags(var,lags);
    local xlags;
    xlags = shiftr((ones(1,lags) .*. var)',seqa(1-lags,1,lags)
                                     .*. ones(cols(var),1),miss(0,0))';
    retp(trimr(var,lags,0),trimr(xlags,0,lags));
endp;