/* 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 MGJOHB: Johansen's ML Estimation of Coefficients Systems */
/* (c) 1996 Michele Gambera and Kristen Strellec                             */
proc (3) = MGJOHB(x,ritardi,caso,cirels);
  /* Maximum Likelihood Cointegration Analysis                               */
  /* 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;                                    */
  /* caso is the case number, optimally selected by MGJOHA: 0 no const,      */
  /*   1 const, 2 trend, 3 const in CI vector, 4 trend in CI vector          */
  /* cirels is the number of CI vectors resulting from MGJOH.                */
/* The procedure returns the coefficients of the ECT, the CI vectors, and a  */
/*   vector containing -T*ln(1-eigenvalue)                                   */

local regres,suu,svv,suv,u,v,ave,avet,eva,eve,obse,leva,alfahat,betahat,
      conto,conta,conte,u1,v1,DXt,DXt_1,a,big_s,pro,w1,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 MGJOHB ***";
"MLE of Cointegrated Systems (Johansen)";"";
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;
obse=rows(DXt);

if caso==0;       @ 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)]';

elseif caso==3;       @ 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)]';

elseif caso==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)]';

elseif caso==4;       @ 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)]';

elseif caso==2;       @ 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)]';

endif;
/* logs of (1-eigenvalue) */
leva=-obse*ln(1-eva);

/* Obtain MLE of a: CI implies that the coefficient of the ECT is ab'  */
/*   and b is given by the eigenvectors                                */
/* normalization */
betahat=orth(eve[.,1:cirels]);
if cols(betahat)<cirels;
  cls;
  "The matrix of the CI vectors is not full rank; terminating procedure.";
  "Returning non-normalized eigenvectors.";
  retp(0,eve[.,1:cirels],leva);
endif;
conto=1;
do while conto<=cirels;
betahat[.,conto]=betahat[.,conto]./sqrt(betahat[.,conto]'svv*betahat[.,conto]);
conto=conto+1;
endo;
/* estimate a */
alfahat=-suv*betahat;
"";
 /* Print output */
output on;
"LAGS:" ritardi;
"alfahat: coefficients of b'Xt in VECM";
alfahat;
"betahat: cointegration vectors, normalized as in Johansen's papers";
betahat;
"-T*ln(1-eigenvalue), to be used in testing";
leva;
"";
retp(alfahat,betahat,leva);
endp;
