/* December 1997 */

/* Programs: mgks.g, mgjoh.g, mgjoha.g, mgjohb.g         */
/* Authors: Michele Gambera with Kristin Strellec        */
/* The code is written and submitted for public, non-commercial use.  */
/* There are no performance guarantees.                  */
/* Please report bugs and suggestions to gambera@psu.edu */
/* Please acknowledge this code (and its authors)        */
/*  if you find it useful in your own work.              */

/* DEFINE PROCEDURE MGJOH: Johansen's ML LAG SPECIFICATION and CI RANK        */
/*   Testing of Cointegrated Systems                                          */
/* (c) 1996 Michele Gambera and Kristen Strellec                              */
proc (0) = MGJOH(x,ritardi,costante);
  /* Maximum Likelihood Cointegration Analysis                                */
  /* this procedure follows Hamilton, section 20.2 and 20.3                   */

  /* x is a matrix with the series to be tested for cointegration             */
  /* ritardi is the max number of lagged deltaX's included in VAR equations   */
  /* costante=0 for no constant in VAR, 1 for constant, 2 for const and trend */
  /* note: 2 implies trend in VAR(deltaX), ie, quadratic trend in levels of X */

local regres,suu,svv,suv,u,v,ave,avet,eva,eve,obse,leva,test1,test2,pro,w1,
      conto,conta,conte,u1,v1,DXt,DXt_1,a,big_s,aka,sik,useful,rt;
"";"*** PROCEDURE MGJOH --- COINTEGRATION ANALYSIS ***";
"Testing for number of lags in VAR (VECM) and CI rank";
"";
/* step 1: auxiliary regressions */
DXt=x-lag1(x);
DXt_1=lag1(x)-lagn(x,2);
/* since the regressors are the same for all equations, there */
/* is no gain in efficiency with simultaneous estimation      */
_olsres=1;
__altnam=0;
__output=1;
__con=0;
pro=seqa(1,1,rows(x));
if costante>=1;
   __con=1;
endif;
conto=1;
regres=DXt_1;
useful=ones(1,3+cols(x));
useful=useful*10000;
do while conto<=ritardi;
 output off;

 conte=1;
 u=zeros(rows(x)-1,1);
 v=u;
 if costante<2;         @ linear or no trend in levels @
   do while conte<=cols(x);
     {a,a,a,a,a,a,a,a,a,u1,a}=ols(0,DXt[.,conte],regres);
     {a,a,a,a,a,a,a,a,a,v1,a}=ols(0,lag1(x[.,conte]),regres);
     conte=conte+1;
     u=u~(zeros(rows(u)-rows(u1),1)|u1);
     v=v~(zeros(rows(v)-rows(v1),1)|v1);
   endo;
   v=v[conto+1:rows(v),2:cols(v)];
 else;                  @ quadratic trend in levels @
   do while conte<=cols(x);
     {a,a,a,a,a,a,a,a,a,u1,a}=ols(0,DXt[.,conte],pro~regres);
     {a,a,a,a,a,a,a,a,a,v1,a}=ols(0,lag1(x[.,conte]),pro~regres);
     conte=conte+1;
     u=u~(zeros(rows(u)-rows(u1),1)|u1);
     v=v~(zeros(rows(v)-rows(v1),1)|v1);
   endo;
 v=v[conto+1:rows(v),2:cols(v)];
 endif;
 u=u[conto+1:rows(u),2:cols(u)];

 /* step 2: canonical correlations */
 obse=rows(u1);
 suu=obse^(-1)*u'*u;
 svv=obse^(-1)*v'*v;
 suv=obse^(-1)*u'*v;
 big_s=invpd(svv)*suv'*invpd(suu)*suv;
 {eva,eve}=eigv(big_s);

 /* Information criteria for lag specification */
 aka=ln(det(big_s))+2*((rows(big_s)^2)*conto/obse);
 sik=ln(det(big_s))+(ln(obse)*(conto*rows(big_s)^2)/obse);

 /* sort eigenvalues in descending order */
 ave=eva'|eve;
 avet=rev(sortc(ave',1));
 eva=avet[.,1];
 eve=avet[.,2:cols(avet)]';

 useful=useful|aka~sik~conto~eva';
conto=conto+1;
endo;
output on;
conto=1;
rt=minindc(useful);

do while conto<3;
 if conto==1;
     eva=useful[rt[1],.]';
     "";"AKAIKE best lag specification of VECM";"Lags:" eva[3];"";
 else;
     eva=useful[rt[2],.]';
     "";"SCHWARZ best lag specification of VECM";"Lags:" eva[3];"";
 endif;
 eva=trimr(eva,3,0);

 /* calculate likelihood ratio tests (Hamilton, p.645) */
 /* test 1; Ho: h-1 CI restrictions exist, Ha: h restrictions exist */
 leva=ln(1-eva);
 test1=-obse*leva;
"";
"'Lambda Max' TEST -- Ho: h-1 CI Vectors, Ha: h CI Vectors (and n-h random
walks):";
"Critical values:";
 if costante==0;
    "see Hamilton, table B11, case 1, or Osterwald-Lenum, table 0";
 elseif costante==1;
    "see Hamilton, table B11, case 3, or Osterwald-Lenum, table 1";
 else;
    "see Osterwald-Lenum, table 2";
 endif;
"Ho: 1 CI vector (and 1 random walk):";
"Test results:";
"Ho:CI Relations                 test";
 seqa(0,1,rows(big_s))~test1;
 "";


 /* test 2; Ho: 0 CI restrictions exist, Ha: h restrictions exist */
 test2=zeros(cols(big_s),1);
 conta=cols(big_s)-1;
 test2[conta]=test1[conta];
 "";
 "'Trace' TEST -- Ho: 0 CI Vectors, Ha: h CI Vectors";
 do while conta>0;
       test2[conta]=test2[conta+1]+test1[conta];
       conta=conta-1;
 endo;
 "Critical values: ";
 "see Hamilton, table B10, same case as above, or Osterwald-Lenum, same table";
 "Ho: No CI relations.";
 "Ha: 1 CI Relation:";
 "Ha: 2 CI relations, etc (n relationships imply all stationary)";
 "Test results:";
 "Ha:CI Relations                 test";
 seqa(1,1,rows(big_s))~test2;
 "";"";

regres=regres~(lagn(x,conto+1)-lagn(x,conto+2));
conto=conto+1;
endo;


endp;
