                  Johansen Test (Isaac & Rapach)

/* ----JOHANSEN.SRC---- */
/* ----Procedures for the Johansen (1988) Cointegration Tests---- */

/***********************************************************************
----Written by Alan G. Isaac and David E. Rapach
----Last revised 23 August 1994
----For public, non-commercial use only
----Please include proper attribution should this code be used in a
    research project
----The authors make no performance guarantees
----This program was written using the The GAUSS System, Version 3.1
      (vendor: Aptech Systems, Inc., Maple Valley, WA)
----Consult Hamilton (1994), tables B.10 & B.11 for critical values (see
    references therein for sources)
*************************************************************************/

/* ----PROCEDURE DEFINITION AREA---- */

/*************************************************************************
----PROC JOHANSEN
----Written by Alan G. Isaac and David E. Rapach
----Last revised 23 August 1994
----FORMAT: {ev,eta0,zeta0}=johansen(x,k)
----INPUT:  x - data matrix (first row is first observation)
            k - maximum lag for VAR model
----OUTPUT: ev - eigenvalues
            eta0 - trace statistic
            zet0 - maximal-eigenvalue statistic
----GLOBAL VARIABLES: none
----EXTERNAL PROCEDURES: varlags
----This follows the presentation on pp. 636-37 in
    Hamilton, James D. _Time Series Analysis_ (Princeton, 1993).

    The original source is
    Johansen, Soren.  'Statistical Analysis of Cointegration Vectors,'
    _Journal of Economic Dynamics and Control_ (v. 12, 1988), pp. 231-54.

----NB: Constant included in regression
************************************************************************/

proc(3)=johansen(x,k);
    local x_lag1,xs,x_1,dx,dxlags,T,bu,u,bv,v,Suu,Suv,Svu,Svv,SuuI,SvvI,A,
    l,z0,eta0,lmax,zeta0;
    x_lag1=x[k+1:rows(x)-1,.];  /* --levels, lag 1-- */
    {xs,x_1}=varlags(x,1);      /* --creating matrices for differencing-- */
    dx=xs-x_1;                  /* --differencing-- */
    {dx,dxlags}=varlags(dx,k);  /* --constructing data matrices-- */
    T=rows(dxlags);             /* --# of observations-- */
    dxlags=ones(T,1)~dxlags;    /* --regressor matrix-- */
    bu=dx/dxlags;               /* --OLS for difference regression-- */
    u=dx-dxlags*bu;             /* --residuals for difference regression-- */
    bv=x_lag1/dxlags;           /* --OLS for level regression-- */
    v=x_lag1-dxlags*bv;         /* --residuals for level regression-- */
    Suu=(u'u)/T;                /* --correlations for diff. residuals-- */
    Suv=(u'v)/T;                /* --cross correlations-- */
    Svu=Suv';                   /* --cross correlations'-- */
    Svv=(v'v)/T;                /* --correlations for level residuals-- */
    SuuI=invpd(Suu);            /* --inverse of duu-- */
    SvvI=invpd(Svv);            /* --inverse of dvv-- */
    A=SvvI*Svu*SuuI*Suv;        /* --matrix A for eigenvalue problem-- */
    l=eig(A);                   /* --eigenvalues for A-- */
    l=rev(sortc(l,1));          /* --ordering eigenvalues-- */
    "eigenvalues:" l';
    z0=ln(1-l);                 /* --obtaining log(1-eigenvalues)-- */
    eta0=-T*sumc(z0);           /* --trace statistic-- */
    "trace statistic:" eta0;
    lmax=l[1,1];                /* --largest eigenvalue-- */
    zeta0=-T*ln(1-lmax);        /* --maximal eigenvalue statistic-- */
    "maximal-eigenvalue statistic:" zeta0;
    retp(l,eta0,zeta0);
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;
