/* 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 MGJOHA: Johansen's ML Testing of Cointegrated Systems    */
/* (c) 1996 Michele Gambera and Kristen Strellec                              */
proc (0) = MGJOHA(x,ritardi,cirels);
  /* Maximum Likelihood Cointegration Analysis: MODEL SELECTION              */
  /* this procedure follows Hamilton, s. 20.3 and Johansen's book (1995)     */

  /* x is a matrix with the series to be tested for cointegration;           */
  /* ritardi is the number of lagged deltaX's included in VAR equations,     */
  /*   optimally selected by SC in MGJOH;                                    */
  /* cirels is the number of CI vectors resulting from MGJOH.                */

local regres,suu,svv,suv,u,v,ave,avet,eva,eve,obse,leva,test1,test2,test3,
      conto,conta,conte,u1,v1,DXt,DXt_1,a,big_s,pro,w1,caso,Xt_1;
/* step 1: auxiliary regressions */
DXt=x-lag1(x);
DXt_1=lag1(DXt);
/* since the regressors are the same for all equations, there */
/* is no gain in efficiency with simultaneous estimation      */
_olsres=1;
__altnam=0;
__output=1;
"";"*** PROCEDURE MGJOHA ***";
"Model Selection Testing (Johansen 1995, chapter 11)";"";
conto=1;
regres=DXt_1;
do while conto<ritardi;
    regres=regres~(lagn(x,conto+1)-lagn(x,conto+2));
    conto=conto+1;
endo;
regres=trimr(regres,ritardi+1,0);
DXt=trimr(DXt,ritardi+1,0);
Xt_1=trimr(lag1(x),ritardi+1,0);
pro=seqa(1,1,rows(regres));
 output off;
 conte=1;
 u=zeros(rows(regres),1);
 v=u;
@ no trend in levels @
/* step 1: auxiliary regressions */
__con=0;
   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,Xt_1[.,conte],regres);
     conte=conte+1;
     u=u~u1;
     v=v~v1;
   endo;
 u=u[.,2:cols(u)];
 v=v[.,2:cols(v)];

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

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

 /* calculate likelihood ratio tests (Hamilton, p.649) */
 leva=ln(1-eva);
 test1=obse*leva;
 test2=sumc(test1[1:cirels]);
 test3=sumc(test1[cirels+1:rows(test1)]);

@ constant in CI vectors @
/* step 1: auxiliary regressions */
 conte=1;
 u=zeros(rows(regres),1);
 v=u;
   __con=0;
   {a,v,a}=olsqr2(ones(rows(regres),1),regres);
   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,Xt_1[.,conte],regres);
     conte=conte+1;
     u=u~u1;
     v=v~v1;
   endo;
 v=v[.,2:cols(v)]~v[.,1];
 u=u[.,2:cols(u)];

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

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

 /* calculate likelihood ratio tests (Hamilton, p.649)      */
 /* note: Johansen proves that the p+1th eigenvalue is zero */
 leva=ln(1-eva);
 test1=obse*leva;
 test2=test2|sumc(test1[1:cirels]);
 test3=test3|sumc(test1[cirels+1:(rows(test1)-1)]);

@ linear trend in levels @
/* step 1: auxiliary regressions */
__con=1;
 conte=1;
 u=zeros(rows(regres),1);
 v=u;
   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,Xt_1[.,conte],regres);
     conte=conte+1;
     u=u~u1;
     v=v~v1;
   endo;
 u=u[.,2:cols(u)];
 v=v[.,2:cols(v)];

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

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

 /* calculate likelihood ratio tests (Hamilton, p.649) */
 leva=ln(1-eva);
 test1=obse*leva;
 test2=test2|sumc(test1[1:cirels]);
 test3=test3|sumc(test1[cirels+1:rows(test1)]);

@ trend in CI vectors @
/* step 1: auxiliary regressions */
 conte=1;
 u=zeros(rows(regres),1);
 v=u;
   __con=1;
   {a,v,a}=olsqr2(pro,regres);
   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,Xt_1[.,conte],regres);
     conte=conte+1;
     u=u~u1;
     v=v~v1;
   endo;
v=v[.,2:cols(v)]~v[.,1];
u=u[.,2:cols(u)];

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

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

 /* calculate likelihood ratio tests (Hamilton, p.649)      */
 /* note: Johansen proves that the p+1th eigenvalue is zero */
 leva=ln(1-eva);
 test1=obse*leva;
 test2=test2|sumc(test1[1:cirels]);
 test3=test3|sumc(test1[cirels+1:(rows(test1)-1)]);

@ quadratic trend in levels @
/* step 1: auxiliary regressions */
__con=1;
 conte=1;
 u=zeros(rows(regres),1);
 v=u;
   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,Xt_1[.,conte],pro~regres);
     conte=conte+1;
     u=u~u1;
     v=v~v1;
   endo;
 u=u[.,2:cols(u)];
 v=v[.,2:cols(v)];

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

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

 /* calculate likelihood ratio tests (Hamilton, p.649) */
 leva=ln(1-eva);
 test1=obse*leva;
 test2=test2|sumc(test1[1:cirels]);
 test3=test3|sumc(test1[cirels+1:rows(test1)]);
"";
 /* Print output */
output on;
"LAGS:" ritardi;"CI Rank assumed:" cirels;
"Testing for model specification:";
"Deterministic parts in CI vectors vs. constant and trends in model";
"For case identification, see Osterwald-Lenum, Oxford Bull. 1992";
"For asymptotic theory, see Johansen's book (1995)";
$"Ho CASE"~"TEST"~"P-value";
caso={"0","1*","1","2*"};
test1=test2;
conto=2;
do while conto<=5;
    test1[conto]=test2[conto-1]-test2[conto];
    conto=conto+1;
endo;
test2=test3;
conto=2;
do while conto<=5;
    test2[conto]=test3[conto]-test3[conto-1];
    conto=conto+1;
endo;

test1=trimr(test1,1,0);
test2=trimr(test2,1,0);
test1[2]=test2[2];
test1[4]=test2[4];
test2=cdfchic(test1,cirels|(cols(x)-cirels)|cirels|(cols(x)-cirels));
call printfmt(rev(caso~test1~test2),0~1~1);

"";

endp;
