/*   PROC JOHANSEN
**
**   Pedro J. F. de Lima, May 1993
**   Department of Economics, The Johns Hopkins University
**      This code is written and submitted for public, non-commercial use.
**      There are no performance guarantees.
**      Please acknowledge this code (and its author)
**      if you find it useful in your own work.
**
**   FORMAT
**          {tr, l_max, z} = johansen(levels,nlags);
**
**   INPUT
**
**         levels -  levels of N series to test for cointegration
**          nlags -  number of lags in the VAR specification
**
**   OUTPUT
**
**            tr  - Trace Statistic, for the hypothesis that there are at
**                  most q=0, 1, 2, ...,n-1 cointegrating vectors against
**                  the completly unrestricted VAR(p) model.
**                  That is Ho: r<=q H1: 0<=q<n.
**
**         l-max  - lmax statistic for the null hypothesis that there are
**                  at most r<=q cointegrating vectors (q=0,1,...,n-1),
**                  against the alternative that only one additional
**                  cointegrating vector exists.
**
**              z -  N x N*3 matrix with the following three N x N matrices:
**                      - E = NxN matrix of normalized eigenvectors (E'Spp E=I),
**                        the first column of E giving the the eigenvector
**                        associated with the largest eigenvalue and so on.
**                      - S00 = NxN matrix with the product moment matrix S00
**                      - S0p = NxN matrix with the product moment matrix S0p
**
**              These matrices are required to compute the cointegrating
**              vectors á, the parameters à, the long run matrix ã=àá',
**              and the variance matrix of the residuals, once r, the # of
**              cointegrating vectors, is defined.
**              See the procedure COINTEGR
*/

proc (3) = johansen(levels,p);
local N,T,diff,x,i,y0t,ypt,q,r0t,rpt,s00,s0p,sp0,spp,c,cc,
lambda,vlambda,lr;

N=cols(levels);
T=rows(levels);
diff=(levels-lag1(levels));
x=diff;
i=1;do while i<=p-1;
x=x~lagn(diff,i);
i=i+1;
endo;
x=packr(x);


y0t=x[.,1:N];
ypt=packr(lagn(levels,p));
x=x[.,N+1:p*N];
q=invpd(x'x);

R0t=y0t-x*q*x'y0t;
Rpt=ypt-x*q*x'ypt;

s00=r0t'r0t/T;
s0p=r0t'rpt/T;
sp0=s0p';
spp=rpt'rpt/T;

c=chol(spp)';
cc=inv(c);
{ lambda, vlambda}=eigrs2(cc*sp0*invpd(s00)*s0p*cc');
retp(rev(-T*cumsumc(ln(1-lambda))), rev(-T*ln(1-lambda)),
(rev((cc'vlambda)')')~s00~s0p);
endp;


/*   PROC CI_VEC
**
**   Pedro J. F. de Lima, May 1993
**   Department of Economics, The Johns Hopkins University
**
**   FORMAT
**          { beta, alpha, lrun, omega} = ci_vec(mat,r);
**
**   INPUT
**          mat - n x 3n  matrix  = to various as defined in proc johansen
**           r  -  number of cointegrating vectors
**
**   OUTPUT
**         beta - estimate of the cointegrating vector á.
**        alpha - estimate of à.
**         lrun - estimate of the long run matrix ã=àá'
**        omega - estimate of ê, the variance matix of the VAR(p) innovations
*/

proc (4) = ci_vec(mat,r);
local e,n,beta,alpha;
beta=mat[.,1:r];
n=rows(mat);
alpha=-mat[.,2*n+1:3*n]*beta;
retp(beta,alpha,-mat[.,2*n+1:3*n]*beta*beta',mat[.,n+1:2*n]-alpha*alpha');
endp;