tolerance=1.e-10;
_svdtol=1.e-10;
cls;
"                                MINIMAL (Version 2)";
" finds the irreducible cointegrating vectors in a set of I(1) variables"
"                    by tests on the Johansen beta matrix";
"                         (James Davidson, 1997)";

/* The data file must be an ASCII file with the data in free format,
   entered either by observation or by variable.
   Substitute the path to your data files in the following line. */

path="c:\\documents\\ox\\japi\\";

/* (Note that GAUSS requires a double backslash to stand for the
   DOS single backslash.)

   The first line of the file must contain two numbers, separated by
   a space: the number of variables in the data set, followed by
   0 if the data are listed by observation, or 1 if listed by variable.

   The second line of the file must contain names for the variables
   (up to 8 characters) separated by spaces.

   The data must be listed beginning at the third line of the file.

   If no file type is given for the data file, ".DAT" is assumed.
   The output file is assigned the same DOS path as the data file.

   If no data file name is entered, the program runs with a simulated
   data set, five I(1) variables with cointegrating rank 2. The output
   is written to file TEST.OUT.

*/

"";
"Enter the name of the data file: ";
filename=CONS;
if STRLEN(filename) GT 0;
  if STRINDX(filename,".",1) EQ 0;
    xfile=path $+ filename $+ ".DAT";
  else;
    xfile=path $+ filename;
  endif;
  "";
  "Enter name and path for output file (just RETURN for <datafile>.OUT)";
  ofname=CONS;
  if STRLEN(ofname) EQ 0;
    dotspot=STRINDX(xfile,".",1);
    if dotspot>0;
      ofname=STRSECT(xfile,1,dotspot-1) $+ ".OUT";
    else;
      ofname=xfile $+ ".OUT";
    endif;
  endif;
  output file=^ofname reset;  
  xfile = "fk.mat";
  load x[]=^xfile;
  m=x[1];
  ov=x[2];
  names=x[3:m+2];
  x=x[m+3:rows(x)];
  nk=rows(x);
  n=int(nk/m);
/* ov = 0 for data by observation, 1 for data by variable*/
  if ov EQ 0;
    "data arranged by observation";
    x=reshape(x,n,m);
  else;
    "data arranged by variable";
    x=reshape(x,m,n);
    x=x';
  endif;
  a=zeros(m,1);
  anames=zeros(m,1);
else;
  ofname="TEST.OUT";
  output file=^ofname reset;
  m=5;
  n=100;
  let names = "X1" "X2" "X3" "X4" "X5" ;
  sigma4=1;
  sigma5=1;
  anvar=5;
  let a = 1 2 3 4 5 ;
  x=RNDN(n,m);
  i=1;
  do while i<=3;
    it=2;
    do while it<=n;
      x[it,i]=x[it,i]+x[it-1,i];
      it=it+1;
    endo;
    i=i+1;
  endo;
  it=1;
  do while it<=n;

    x[it,4]=x[it,4]*sigma4+x[it,1]+x[it,2]+x[it,3];
    x[it,5]=x[it,5]*sigma5+x[it,1]-x[it,2];

    it=it+1;
  endo;
  anvar=5;
  anames=names;
endif;

osltr=zeros(45,4);
call gentr;
oslmx=zeros(45,4);
call genmx;
chisq=zeros(15,4);
call genchi;
kpcrit=zeros(2,4);
call genkpss;
dfcrit=zeros(1,6);
call gendfcrit;

FORMAT /LD 4,0;
 "number of observations: ";;n;;", number of variables: ";;m;
FORMAT /LD 10,4;

intcep=ones(n,1);
trendv=seqa(1,1,n);
nsq=n*n;

/* KPSS and PP tests of I(0) versus I(1) */

i=1;
anvar=0;
FORMAT /RD 12,4;

  "       Tests for unit roots in the data series:";
  "    Rejections at: *=10%, **=5%, ***=2.5%, ****=1%.";
  "-Variable-|----KPSS-[H0 = I(0)]--------|-Phillips-Perron--[H0 = I(1)]--|-";
  "----------|---Trend-------NoTrend------|---Trend-------NoTrend---------|-";

  i=1;
  let oput[1,10] = 0 0 0 0 0 0 0 0 0 0 ;
  let mask[1,10] = 0 1 0 1 0 0 1 0 1 0 ;
  do while i<=m;
    c0=names[i];
    o1=kpss(1,i,1,1);
    c1=doi1test(o1,1);
    o2=kpss(1,i,1,0);
    c2=doi1test(o2,2);
    o3=pptest(1,i,1,1,1);
    c3=doi1test(o3,3);
    o4=pptest(1,i,1,0,1);
    c4=doi1test(o4,4);
    rw=c0~o1~c1~o2~c2~"    "~o3~c3~o4~c4;
    oput=oput|rw;
    i=i+1;
  endo;
oput=oput[2:rows(oput),.];
let fs[1,3]= "*.*lf  " 10 10 ;
let fd[1,3]= "*.*lf" 8 3 ;
let fc[1,3]= "*.*lf " 4 4 ;
ft=fs|fd|fc|fd|fc|fc|fd|fc|fd|fc;
y = PRINTFM(oput,mask,ft);
"If you wish to delete one or more variables from the CVAR, type the";
"names one at a time, followed by Enter. To finish, type Enter again.";

delvar=getnames;
a=zeros(m,1);
anames=zeros(m,1);
i=1;
anvar=0;
do while i<=m;
  j=2;
  delit=0;
  do while j<=rows(delvar);
    if upper(delvar[j]) $== upper(names[i]);
      delit=1;
    endif;
    j=j+1;
  endo;
  if delit==0;
    anvar=anvar+1;
    a[anvar]=i;
    anames[anvar]=names[i];
  endif;
  i=i+1;
endo;

if anvar LE 1;
  "Nothing to do, program terminating.";
  end;
endif;
a=a[1:anvar];
anames=anames[1:anvar];

FORMAT /RD 4,0;
anvar;;"variables retained for cointegration tests.";
FORMAT /LZ 9,0;
$anames';

"";
"Akaike & Schwarz Criteria for choice of lags:";
klag=0;
limit=(n-1)/(3*anvar);
if limit>12;
  limit=12;
endif;
do while limit>1;
  if limit*n*anvar>140000;
    limit=limit-1;
  else;
    break;
  endif;
endo;
do while (klag+1) LE limit;
  klag1=klag+1;
  nk1=n-klag1;
  like=loglik(a,anvar,1);
  akaike=like-klag1*anvar^2-anvar;
  schwarz=like-(klag1*anvar^2+anvar)*ln(nk1);
FORMAT /LZ 2,0;
  klag;;"lags: AIC = ";;
FORMAT /LZ 10,4;
akaike;;"  SIC = ";;schwarz;
  klag=klag+1;
endo;
"How many lagged differences do you want to include in the CVAR?";
ans=CONS;
if ans $/= "";
  klag=STOF(ans);
else;
  klag=0;
endif;
if klag GT n/2 OR klag LT 0;
  klag=0;
endif;
klag;

"";
FORMAT /LZ 8,4;
"The eigenvalues in the reduced rank regression in ascending order are,";
  { tracet,maxlamt,lamda }=coranktest(a,anvar,1);
"   1) assuming a deterministic trend:";
lamda';
  { trace0,maxlam0,lamda }=coranktest(a,anvar,0);
"   2) assuming no deterministic trend:";
lamda';
 "";
  "    Johansen tests of H0: rank = r vs. H1: rank > r,";
  "                       + Chisq(p-r) test for trend, given rank = r.";
  "    Rejections at: *=10%, **=5%, ***=2.5%, ****=1%.";
  "--r--|--------Lamdamax----------|-------Trace-------------|--Trend test----";
  "-----|---Trend-----NoTrend------|--Trend-----NoTrend------|----------------";

  i=1;
  let oput[1,11]= 0 0 0 0 0 0 0 0 0 0 0 ;
  let mask[1,11] = 1 1 0 1 0 1 0 1 0 1 0;
  do while i<=anvar;
    o0=anvar-i;
    o1=maxlamt[i];
    c1=dotest(o1,i,1);
    o2=maxlam0[i];
    c2=dotest(o2,i,2);
    o3=tracet[i];
    c3=dotest(o3,i,3);
    o4=trace0[i];
    c4=dotest(o4,i,4);
    o5=o4-o3;
    c5=dotest(o5,i,5);
    rw=o0~o1~c1~o2~c2~o3~c3~o4~c4~o5~c5;

    oput=oput|rw;
    i=i+1;
  endo;
  oput=oput[2:rows(oput),.];
  let fr[1,3]= "*.*lf  " 3 0 ;
  let fd[1,3]= "*.*lf" 8 3 ;
  let fc[1,3]= "*.*lf " 4 4 ;
  ft=fr|fd|fc|fd|fc|fd|fc|fd|fc|fd|fc;
  y = PRINTFM(oput,mask,ft);
   "";
"Enter the cointegrating rank";
ans=CONS;
ans;
corank=STOF(ans);

"Assume a trend in the model? (Type [Y] or N):";
ans=UPPER(CONS);
trend=1;
if ans $== "N";
  trend=0;
else;
  ans="Y";
endif;
$ans;
  { lamdas,nus,s00,s10,s11,s11inv,dx,x1 } = joheigen(a,anvar,trend);
  { aa,bigb,beta,beta0 } =
           testcomps(lamdas,nus,s00,s10,s11,s11inv,corank,anvar);
  FORMAT /LZ 9,3;
  "The beta' matrix is";
  if trend==1;
     $anames';
     beta';
  else;
     $anames'~"Const";
     beta'~beta0';
  endif;

"";
testype=1;
"Significance level (%) of the cointegration tests? Enter 1, 2.5, [5] or 10.";
"Alternatively, enter T to perform a single cointegration test,";
"                  or A to print ALL the tests.";
ans=upper(CONS);
if ans $== "1";
  sigra=4;
elseif ans $== "2.5";
  sigra=3;
elseif ans $== "10";
  sigra=1;
elseif ans $== "T";
  testype=2;
elseif ans $== "A";
  testype=3;
else;
  sigra=2;
  ans="5";
endif;
$ans;
if testype==1;
  "Use rule of thumb correction to critical values? (Enter Y or [N])";
  ans=upper(CONS);
  if ans $== "Y";
    corrfact=1;
  else;
    corrfact=0;
  ans="N";
  endif;
  $ans;
  "";
  ncoint=0;
  coints=zeros(10*anvar,m+1);
  yes=elim(a,anames,anvar,1);
  if yes EQ 0;
     "All restrictions on Beta rejected";
  endif;
elseif testype==2;
  ncoint=1;
  b=zeros(m,1);
  bnames=zeros(m,1);
  "Type the names of the variables you want to include in your vector";
  "one at a time, followed by Enter. To finish, type Enter again.";

  delvar=getnames;
  i=1;
  bnvar=0;
  do while i<=m;
    j=2;
    incl=0;
    do while j<=rows(delvar);
      if upper(delvar[j]) $== upper(names[i]);
        incl=1;
      endif;
      j=j+1;
    endo;
    if incl==1;
      bnvar=bnvar+1;
      b[bnvar]=i;
      bnames[bnvar]=names[i];
    endif;
    i=i+1;
  endo;
  b=b[1:bnvar];
  bnames=bnames[1:bnvar];
  b';
  $bnames';
  { stat,foundvec,foundse,degf }=ranktest(b,bnvar);
  call writeres(stat,b,bnames,bnvar,
                    degf,"-",foundvec,foundse);
else;
  ncoint=0;
  call alltest(a,anames,anvar,1);
endif;
"Tests complete";
end;

proc(1)=getnames;
local delvar,ans;
  delvar="";
  nameloop:
  ans=CONS;
  "";
  if (ans $== "") OR (rows(delvar)>=m);
    goto contlabel;
  else;
    delvar=delvar|ans;
  endif;
  goto nameloop;
  contlabel:
retp(delvar);
endp;

proc(1) = doi1test(stat,k);
local i,j,code,result,table;
  let code= "*" "**" "***" "****";
  result = "";
  if k<=2;
    j=1;
    do while j<=4;
      table=kpcrit[2-k,j];
      if stat>table;
        result=code[j];
      endif;
      j=j+1;
    endo;
  else;
    j=1;
    do while j<=3;
      table=dfcrit[1,3*(4-k)+(4-j)];
      if stat<table;
        if j==3;
          result=code[j+1];
        else;
          result=code[j];
        endif;
      endif;
      j=j+1;
    endo;
  endif;
retp(result);
endp;

proc(1) = dotest(stat,i,k);
local j,code,result,table;
  let code= "*" "**" "***" "****";
  result = "";
  j=1;
  do while j<=4;
    if k EQ 1;
       table=oslmx[30+i,j];
    elseif k EQ 2;
       table=oslmx[15+i,j];
    elseif k EQ 3;
       table=osltr[30+i,j];
    elseif k EQ 4;
       table=osltr[15+i,j];
    else;
       table=chisq[i,j];
    endif;
    if stat>table;
      result=code[j];
    endif;
    j=j+1;
  endo;
retp(result);
endp;


/* This is a self-calling subroutine which eliminates one
variable at a time from the set, and tests for cointegration
in each subset.     */

proc elim(a,anames,anvar,ia);
local ib,b,bnames,bnvar,coivec,result;
  coivec=0;
  ib=ia;
  bnvar=anvar-1;
/* eliminate one variable at a time from a */

  do while ib<=anvar;
    result=0;
    { b,bnames }=getvec(a,anames,anvar,ib);
    if bnvar > 2;
      if elim(b,bnames,bnvar,ib) EQ 0;
        result=testb(b,bnames,bnvar);
      else;
        result=1;
      endif;
    else;
      result=testb(b,bnames,bnvar);
    endif;
    if result EQ 1 and testype NE 3;
      coivec=1;
    endif;
      /* ... and move to the next deletion. */
    ib=ib+1;
  endo;
retp(coivec);
endp;

proc testb(b,bnames,bnvar);
local stat,crit,coivec,nvar,foundvec,foundse,degf;
local coincorr,kpsstat,kpsscrit;
  coivec=0;
  if checkco(b,bnvar) EQ 0;
    { stat,foundvec,foundse,degf }=ranktest(b,bnvar);
    if corrfact==1;
      if bnvar==2;
        coincorr=4;
      elseif bnvar==3;
        coincorr=2;
      else;
        coincorr=1;
      endif;
    else;
      coincorr=1;
    endif;
    if degf>0;
      crit=chisq[degf,sigra]*coincorr;
    else;
       crit=0;
    endif;
    if stat<=crit;
      call storeco(b,bnvar);
      call writeres(stat,b,bnames,bnvar,
                  degf,crit,foundvec,foundse);
      coivec=1;
    endif;
  else;
    coivec=1;
  endif;
retp(coivec);
endp;


proc(0)=writeres(stat,b,bnames,bnvar,
                    degf,crit,foundvec,foundse);
local dropone,cnvar,ib,c,cnames,ppstat,ppcrit;
local temp1,temp2,temp3;
  print "";
  format /RZ 4,4;
  if testype/=2;
     ncoint;;":";;
  endif;
  " ChisqStat(";;degf;;") = ";;
  format /LZ 6,3;
  stat;;"(";;crit;;")";
  if bnvar GT 2;
    dropone=zeros(bnvar,1);
    cnvar=bnvar-1;
    ib=1;
    do while ib LE bnvar;
      { c,cnames }=getvec(b,bnames,bnvar,ib);
      { stat,temp1,temp2,degf }=ranktest(c,cnvar);
      dropone[ib]=stat;
      ib=ib+1;
    endo;
  endif;
  FORMAT /LZ 9,0;
  if trend EQ 0;
    print   "Variables:      ";; $bnames'~" Const";
  else;
    print   "Variables:     ";; $bnames';
  endif;
  FORMAT /LZ 9,3;
  print     "Coefficients:  ";; foundvec';
  print     "Std. errors:   ";; foundse';
  if bnvar GT 2;
    FORMAT /LZ 9,4;
    print   "Omission tests:";; dropone';
  endif;
  ppstat=pptest(foundvec,b,bnvar,trend,1);
  ppcrit=dfcrit[1,2+trend*3];
  print "Phillips-Perron statistic= ";;ppstat;;
              "(5% critical value =";;ppcrit;;")";
@    kpsstat=kpss(foundvec,b,bnvar,trend);
     kpsscrit=kpcrit[trend+1,2];
     print "           KPSS statistic= ";;kpsstat;;
              "(5% critical value = ";;kpsscrit;;")";
@
endp;

proc(2) = getvec(b,bnames,bnvar,ic);
local c,cnames;
  if ic EQ 1;
    c=b[2:bnvar];
    cnames=bnames[2:bnvar];
  elseif ic EQ bnvar;
    c=b[1:bnvar-1];
    cnames=bnames[1:bnvar-1];
  else;
    c=b[1:ic-1]|b[ic+1:bnvar];
    cnames=bnames[1:ic-1]|bnames[ic+1:bnvar];
  endif;
retp(c,cnames);
endp;

proc(0)= storeco(b,bnvar);
local i;
  ncoint=ncoint+1;
  i=1;
  do while i LE bnvar;
    coints[ncoint,i]=b[i];
    i=i+1;
  endo;
  coints[ncoint,m+1]=bnvar;
endp;

proc checkco(b,bnvar);
local i,j,k,check;
  check=0;
  k=1;

  /* look at each detected CIV */

  do while k LE ncoint AND check EQ 0;
    j=1;

     /* work through the non-zero elements.  */

    do while j LE coints[k,m+1];

      /* Compare each element with the elements of b.
         If there are as many matches as elements, we know kth CIV
         is in b */
      i=1;
      do while i LE bnvar;
        if b[i] EQ coints[k,j];
          check=check+1;
        endif;
        i=i+1;
      endo;
      j=j+1;
    endo;
    if check LT coints[k,m+1];
      check=0;
    endif;
    k=k+1;
  endo;
  retp(check);
endp;


proc(0)=alltest(a,anames,anvar,ia);
local ib,b,bnames,bnvar,coivec,result;
local stat,foundvec,foundse,degf;
  coivec=0;
  ib=ia;
  bnvar=anvar-1;
/* eliminate one variable at a time from a */
  do while ib<=anvar;
    result=0;
    { b,bnames }=getvec(a,anames,anvar,ib);
    if bnvar > 2;
      call alltest(b,bnames,bnvar,ib);
    endif;
    { stat,foundvec,foundse,degf }=ranktest(b,bnvar);
    ncoint=ncoint+1;
    call writeres(stat,b,bnames,bnvar,
                    degf,"-",foundvec,foundse);
    ib=ib+1;
  endo;
endp;

proc(4) = ranktest(b,bnvar);
local i,aitch,stat,gee,rankcov,invcov,foundvec,foundse,foundcov;
local oldcov,cov,testvec,testmat,testnus,testlams,testev,nrestr;
local pee,rankpee,queue,pbetah,qa,qpbetah,iota;
  { gee,aitch } = restric(b,bnvar);
  nrestr=rows(aitch);
  iota=ones(nrestr,1);
  testmat=beta'aitch'aitch*beta;
  { testlams,testnus }=EIGHV(testmat);
  testev=testnus[.,1];
  testvec=aitch*beta*testev;
  if corank>1;
    pee=eye(corank)*testlams[1]-testmat;
    pee=pinv(pee);
    queue=aitch'aitch*beta;
    pbetah=pee*beta'aitch';
    qa=queue*testev;
    qpbetah=queue*pbetah;
    cov=gencov(testev,aitch,pbetah,qa,qpbetah,nrestr);
  else;
    cov=(testev'aa*testev)*aitch*bigb*aitch';
  endif;
  if iota'cov*iota>1e-12;
    { invcov,rankcov }=geninv(cov);
    stat=nsq*testvec'invcov*testvec;
  else;
    stat=0;
    rankcov=0;
  endif;
  if bnvar>0;
    foundvec=gee*beta*testev;
    if trend EQ 0;
      foundvec=foundvec|(beta0*testev);
    endif;
    if corank>1;
      pbetah=pee*beta'gee';
      qpbetah=queue*pbetah;
      foundse=diag(gencov(testev,gee,pbetah,qa,qpbetah,bnvar));
    else;
      foundse=diag((testev'aa*testev)*gee*bigb*gee');
    endif;
    foundse=foundse^(.5)/n;
  else;
    foundvec=0;
  endif;
retp(stat,foundvec,foundse,rankcov);
endp;

proc(1)=gencov(testev,aitch,pbetah,qa,qpbetah,nrestr);
local i,j,cov;
local a11,a12,a13,a22,a23,a33,b11,b12,b13,b22,b23,b33;
  a11=testev'aa*testev;
  a12=testev'aa*pbetah;
  a13=a11;
  a22=pbetah'aa*pbetah;
  a23=a12;
  a33=a11;
  b11=aitch*bigb*aitch';
  b12=aitch*bigb*qa;
  b13=aitch*bigb*qpbetah;
  b22=qa'bigb*qa;
  b23=qa'bigb*qpbetah;
  b33=qpbetah'bigb*qpbetah;
  cov=zeros(nrestr,nrestr);
  i=1;
  do while i<=nrestr;
    j=1;
    do while j<=i;
      cov[i,j]=a11*b11[i,j]+a12[i]*b12[j]+a13*b13[i,j]
              +a12[j]*b12[i]+a22[i,j]*b22+a23[i]*b23[j]
              +a13*b13[j,i]+a23[j]*b23[i]+a33*b33[i,j];
      cov[j,i]=cov[i,j];
      j=j+1;
    endo;
    i=i+1;
  endo;
retp(cov);
endp;

proc(2) = geninv(a);
local lamdas,nus,i,rank,dimn,ainv,laminv;
  a=a*1;
  { lamdas,nus } = EIGHV(a);
  i=1;
  dimn=cols(a);
  rank=dimn;
  do while i<=dimn;
    if lamdas[i] GT tolerance;
      break;
    else;
      rank=rank-1;
    endif;
    i=i+1;
  endo;
  if rank>0;
    nus=nus[.,dimn-rank+1:dimn];
    laminv=DIAGRV(ZEROS(rank,rank),1/lamdas[dimn-rank+1:dimn]);
    ainv=nus*laminv*nus';
  else;
    ainv=ones(dimn,dimn);
  endif;
retp(ainv,rank);
endp;

proc(8)=joheigen(b,bnvar,trend);
local q,i,j,k,nk,nk2,n1,x1,dx,lags,temp,temp2,l11,lamdas,nus;
local s00,s10,s11,s11inv;
  nk=n-klag;
  nk2=n-klag-2;
  n1=n-1;
  x1=zeros(n-klag-1,1);
  dx=zeros(n1,1);
  lags=ones(n-klag-1,1);
  j=1;
  do while j<=bnvar;
    i=b[j];
    x1=x1~x[1:nk-1,i];
    dx=dx~(x[2:n,i]-x[1:n1,i]);
    k=1;
    do while k<=klag;
      lags=lags~dx[k:nk2+k,j+1];
      k=k+1;
    endo;
    j=j+1;
  endo;
  if trend EQ 0;
    x1=x1[.,2:cols(x1)]~ones(n-klag-1,1);
    if klag>0;
     lags=lags[.,2:cols(lags)];
    endif;

    /* this gives us QbarH */

  else;
    x1=x1[.,2:cols(x1)];

   /* this gives us QbarG and QbarW */

  endif;
  dx=dx[klag+1:n1,2:cols(dx)];
  if (klag>0) OR (trend>0);
    temp=lags*invpd(lags'lags);
    x1=x1-temp*lags'x1;
    dx=dx-temp*lags'dx;
  endif;
  s00=dx'dx/n;
  s11=x1'x1/n;
  s10=x1'dx/n;
  s11inv=invpd(s11);
  l11=chol(s11inv);
  temp2=l11*s10;
  q=temp2*invpd(s00)*temp2';
  { lamdas,nus } = EIGHV(q);
  nus=l11'nus;
retp(lamdas,nus,s00,s10,s11,s11inv,dx,x1);
endp;

proc(4) = testcomps(lamdas,nus,s00,s10,s11,s11inv,corank,anvar);
local aa,lam,lams;
local beta0,c,hbh,itemp,nuhats,cov,testvec,testev,beta,alpha,biglam,biglaminv;
  FORMAT /LZ 10,6;
  beta=nus[.,anvar-corank+2-trend:anvar+1-trend];
  alpha=s10'beta;
  biglam=s00-alpha*beta's10;
  biglaminv=invpd(biglam);
  if trend EQ 0;
    beta0=beta[anvar+1,.];
    beta=beta[1:anvar,.];
    s11=s11[1:anvar,1:anvar];
    s11inv=s11inv[1:anvar,1:anvar];
    nuhats=nus[1:anvar,2:anvar-corank+1];
  else;
    beta0=0;
    nuhats=nus[.,1:anvar-corank];
  endif;
  c=beta*invpd(beta'beta);
  bigb=n*(eye(anvar)-beta*c')*S11inv*(eye(anvar)-c*beta');
  lams=lamdas[anvar-corank+2-trend:anvar+1-trend];
  lam=zeros(corank,corank);
  lam=DIAGRV(lam,lams);
  aa=invpd(lam)-eye(corank);
retp(aa,bigb,beta,beta0);
endp;

proc(2) = restric(b,bnvar);
local gee,omits,i,j,k,aitch;
  omits=anvar-bnvar;
  aitch=zeros(omits,anvar);
  gee=zeros(bnvar,anvar);
  j=1;
  k=1;
  i=1;
  do while j<=anvar;
    if b[k] NE a[j];
      aitch[i,j]=1;
      i=i+1;
    else;
      gee[k,j]=1;
      if k<bnvar;
        k=k+1;
      endif;
    endif;
    j=j+1;
  endo;
retp(gee,aitch);
endp;

proc(3) = coranktest(a,anvar,trend);
local lamdas,nus,s00,s10,s11inv,sx,x1,i,traceq,maxlam;
  { lamdas,nus,s00,s10,s11,s11inv,dx,x1 } = joheigen(a,anvar,trend);
  traceq=zeros(anvar,1);
  maxlam=zeros(anvar,1);
  maxlam[1]=-(n-klag-1)*ln(1-lamdas[2-trend]);
  traceq[1]=maxlam[1];
  i=2;
  do while i<=anvar;
    maxlam[i]=-(n-klag-1)*ln(1-lamdas[i-trend+1]);
    traceq[i]=traceq[i-1]+maxlam[i];
    i=i+1;
  endo;
retp(traceq,maxlam,lamdas);
endp;
end;

proc(1) = loglik(a,anvar,trend);
local biglam,loglik,lamdas,nus,s00,s10,s11inv,sx,x1,i;
local traceq,maxlam,logl2;
  { lamdas,nus,s00,s10,s11,s11inv,dx,x1 } = joheigen(a,anvar,trend);
  biglam=s00-s10'*s11inv*s10;
  loglik=-(n-klag-1)*ln(det(biglam));
retp(loglik);
endp;

proc(1)=kpss(foundvec,b,bnvar,trend);
local e,eta,bigs,etahat,t,sigma,lamda,omega;
  { e } = getres(foundvec,b,bnvar,trend);
  eta=0;
  t=1;
  bigs=0;
  do while t<=n;
    bigs=bigs+e[t];
    eta=eta+bigs^2;
    t=t+1;
  endo;
  { sigma,lamda,omega }=lrcov(e);
  etahat=eta/(n^2*omega);
retp(etahat);
endp;

proc(3)= lrcov(x);
local n,m,el,i,q0,qi,qqi,sigma,lamda,omega;
  n=rows(x);
  m=cols(x);
  sigma=x'x/n;
  el=int(4*(n/100)^(1/4));
  i=1;
  lamda=zeros(m,m);
  do while i <= el;
    q0=x[i+1:n,.];
    qi=x[1:n-i,.];
    qqi=(qi'q0)*(1-i/(el+1))/n;
    lamda=lamda+qqi;
    i=i+1;
  endo;
  omega = sigma+lamda+lamda';
  retp(sigma,lamda,omega);
endp;

proc(1)=pptest(foundvec,b,bnvar,trend,whichz);
local i,j,m,u,alfa,ulag,ucur,uu,stat;
local v,sk2,stl2,stl,lamda,squu,ztee;

{ u } = getres(foundvec,b,bnvar,trend);

/* Phillips-Perron Z-hat-t and z-hat-alfa statistics */

  ulag=u[1:n-1];
  ucur=u[2:n];
  alfa=ulag'ucur/(ulag'ulag);
  v=ucur-alfa*ulag;
  { sk2,lamda,stl2 } = lrcov(v);
  uu=u'u;
  squu=uu^0.5;
  stl=stl2^0.5;
  if whichz EQ 0;
      stat=n*(alfa-1)-(stl2-(sk2*n^2)/uu)/2;
  else;
    stat=squu*(alfa-1)/stl-lamda/(stl*squu/n);
  endif;
  retp(stat);
endp;

proc(1)=getres(foundvec,b,bnvar,trend);
local i,j,z,res;
  z=zeros(n,1);
  j=1;
  do while j<=bnvar;
    i=b[j];
    z=z~x[.,i];
    j=j+1;
  endo;
  z=z[.,2:cols(z)];
  if trend EQ 0;
    foundvec=foundvec[1:bnvar];
  endif;
  res=z*foundvec;
  z=intcep;
  if trend EQ 1;
    z=z~trendv;
  endif;
  res=res-z*invpd(z'z)*z'res;
retp(res);
endp;

proc(0)=genmx;
/* cases 12-15 have been generated by extrapolation using
   cubic trend regressions */
/* QbarW */
oslmx[1,1]=6.50;   oslmx[1,2]=8.18;   oslmx[1,3]=9.72;   oslmx[1,4]=11.65;
oslmx[2,1]=12.91;  oslmx[2,2]=14.90;  oslmx[2,3]=17.07;  oslmx[2,4]=19.19;
oslmx[3,1]=18.90;  oslmx[3,2]=21.07;  oslmx[3,3]=22.89;  oslmx[3,4]=25.75;
oslmx[4,1]=24.78;  oslmx[4,2]=27.14;  oslmx[4,3]=29.16;  oslmx[4,4]=32.14;
oslmx[5,1]=30.84;  oslmx[5,2]=33.32;  oslmx[5,3]=35.80;  oslmx[5,4]=38.78;
oslmx[6,1]=36.35;  oslmx[6,2]=39.43;  oslmx[6,3]=41.86;  oslmx[6,4]=44.59;
oslmx[7,1]=42.06;  oslmx[7,2]=44.91;  oslmx[7,3]=47.59;  oslmx[7,4]=51.30;
oslmx[8,1]=48.43;  oslmx[8,2]=51.07;  oslmx[8,3]=53.85;  oslmx[8,4]=57.07;
oslmx[9,1]=54.01;  oslmx[9,2]=57.00;  oslmx[9,3]=59.80;  oslmx[9,4]=63.37;
oslmx[10,1]=59.19; oslmx[10,2]=62.42; oslmx[10,3]=64.98; oslmx[10,4]=68.61;
oslmx[11,1]=65.07; oslmx[11,2]=68.27; oslmx[11,3]=70.69; oslmx[11,4]=74.36;
oslmx[12,1]=71.00; oslmx[12,2]=74.00; oslmx[12,3]=76.00; oslmx[12,4]=80.00;
oslmx[13,1]=76.00; oslmx[13,2]=80.00; oslmx[13,3]=81.00; oslmx[13,4]=85.00;
oslmx[14,1]=82.00; oslmx[14,2]=86.00; oslmx[14,3]=86.00; oslmx[14,4]=91.00;
oslmx[15,1]=87.00; oslmx[15,2]=91.00; oslmx[15,3]=91.00; oslmx[15,4]=96.00;

/* QbarH - add 15 */
oslmx[16,1]=7.52;  oslmx[16,2]=9.24;  oslmx[16,3]=10.80; oslmx[16,4]=12.97;
oslmx[17,1]=13.75; oslmx[17,2]=15.67; oslmx[17,3]=17.63; oslmx[17,4]=20.20;
oslmx[18,1]=19.77; oslmx[18,2]=22.00; oslmx[18,3]=24.07; oslmx[18,4]=26.81;
oslmx[19,1]=25.56; oslmx[19,2]=28.14; oslmx[19,3]=30.32; oslmx[19,4]=33.24;
oslmx[20,1]=31.66; oslmx[20,2]=34.40; oslmx[20,3]=36.90; oslmx[20,4]=39.79;
oslmx[21,1]=37.45; oslmx[21,2]=40.30; oslmx[21,3]=43.22; oslmx[21,4]=46.82;
oslmx[22,1]=43.25; oslmx[22,2]=46.45; oslmx[22,3]=48.99; oslmx[22,4]=51.91;
oslmx[23,1]=48.91; oslmx[23,2]=52.00; oslmx[23,3]=54.71; oslmx[23,4]=57.95;
oslmx[24,1]=54.35; oslmx[24,2]=57.42; oslmx[24,3]=60.50; oslmx[24,4]=63.71;
oslmx[25,1]=60.25; oslmx[25,2]=63.57; oslmx[25,3]=66.24; oslmx[25,4]=69.94;
oslmx[26,1]=66.02; oslmx[26,2]=69.74; oslmx[26,3]=72.64; oslmx[26,4]=76.73;
oslmx[27,1]=72.00; oslmx[27,2]=76.00; oslmx[27,3]=78.00; oslmx[27,4]=83.00;
oslmx[28,1]=78.00; oslmx[28,2]=82.00; oslmx[28,3]=85.00; oslmx[28,4]=90.00;
oslmx[29,1]=83.00; oslmx[29,2]=88.00; oslmx[29,3]=91.00; oslmx[29,4]=97.00;
oslmx[30,1]=89.00; oslmx[30,2]=95.00; oslmx[30,3]=97.00; oslmx[30,4]=105.00;
/* QbarG - add 30 */
oslmx[31,1]= 2.69; oslmx[31,2]= 3.76; oslmx[31,3]= 4.95; oslmx[31,4]= 6.65;
oslmx[32,1]=12.07; oslmx[32,2]=14.07; oslmx[32,3]=16.05; oslmx[32,4]=18.63;
oslmx[33,1]=18.60; oslmx[33,2]=20.97; oslmx[33,3]=23.09; oslmx[33,4]=25.52;
oslmx[34,1]=24.73; oslmx[34,2]=27.07; oslmx[34,3]=28.98; oslmx[34,4]=32.24;
oslmx[35,1]=30.90; oslmx[35,2]=33.46; oslmx[35,3]=35.71; oslmx[35,4]=38.77;
oslmx[36,1]=36.76; oslmx[36,2]=39.37; oslmx[36,3]=41.86; oslmx[36,4]=45.10;
oslmx[37,1]=42.32; oslmx[37,2]=45.28; oslmx[37,3]=47.96; oslmx[37,4]=51.57;
oslmx[38,1]=48.33; oslmx[38,2]=51.42; oslmx[38,3]=54.29; oslmx[38,4]=57.69;
oslmx[39,1]=53.98; oslmx[39,2]=57.12; oslmx[39,3]=59.33; oslmx[39,4]=62.80;
oslmx[40,1]=59.62; oslmx[40,2]=62.61; oslmx[40,3]=65.44; oslmx[40,4]=69.09;
oslmx[41,1]=65.38; oslmx[41,2]=68.83; oslmx[41,3]=72.11; oslmx[41,4]=75.95;
oslmx[42,1]=72.00; oslmx[42,2]=76.00; oslmx[42,3]=80.00; oslmx[42,4]=84.00;
oslmx[43,1]=79.00; oslmx[43,2]=84.00; oslmx[43,3]=88.00; oslmx[43,4]=92.00;
oslmx[44,1]=87.00; oslmx[44,2]=92.00; oslmx[44,3]=97.00; oslmx[44,4]=102.00;
oslmx[45,1]=96.00; oslmx[45,2]=102.00;oslmx[45,3]=102.00;oslmx[45,4]=114.00;
endp;

proc(0)=gentr;
/* cases 12-15 have been generated by extrapolation using
   cubic trend regressions */
/* QbarW */
osltr[1,1]=6.50;   osltr[1,2]=8.18;   osltr[1,3]=9.72;   osltr[1,4]=11.65;
osltr[2,1]=15.66;  osltr[2,2]=17.95;  osltr[2,3]=20.08;  osltr[2,4]=23.52;
osltr[3,1]=28.71;  osltr[3,2]=31.52;  osltr[3,3]=34.48;  osltr[3,4]=37.22;
osltr[4,1]=45.23;  osltr[4,2]=48.28;  osltr[4,3]=51.54;  osltr[4,4]=55.43;
osltr[5,1]=66.49;  osltr[5,2]=70.60;  osltr[5,3]=74.04;  osltr[5,4]=78.87;
osltr[6,1]=90.39;  osltr[6,2]=95.18;  osltr[6,3]=99.32;  osltr[6,4]=104.20;
osltr[7,1]=118.99; osltr[7,2]=124.25; osltr[7,3]=129.75; osltr[7,4]=136.06;
osltr[8,1]=151.38; osltr[8,2]=157.11; osltr[8,3]=162.75; osltr[8,4]=168.92;
osltr[9,1]=186.54; osltr[9,2]=192.84; osltr[9,3]=198.06; osltr[9,4]=204.79;
osltr[10,1]=226.34;osltr[10,2]=232.49;osltr[10,3]=238.26;osltr[10,4]=246.27;
osltr[11,1]=269.53;osltr[11,2]=277.39;osltr[11,3]=283.84;osltr[11,4]=292.65;
osltr[12,1]=316.00;osltr[12,2]=325.00;osltr[12,3]=331.00;osltr[12,4]=341.00;
osltr[13,1]=367.00;osltr[13,2]=376.00;osltr[13,3]=382.00;osltr[13,4]=393.00;
osltr[14,1]=421.00;osltr[14,2]=431.00;osltr[14,3]=437.00;osltr[14,4]=449.00;
osltr[15,1]=479.00;osltr[15,2]=489.00;osltr[15,3]=495.00;osltr[15,4]=509.00;

/* QbarH - add 11 */
osltr[16,1]=7.52;  osltr[16,2]=9.24;  osltr[16,3]=10.80; osltr[16,4]=12.97;
osltr[17,1]=17.85; osltr[17,2]=19.96; osltr[17,3]=22.05; osltr[17,4]=24.60;
osltr[18,1]=32.00; osltr[18,2]=34.91; osltr[18,3]=37.61; osltr[18,4]=41.07;
osltr[19,1]=49.65; osltr[19,2]=53.12; osltr[19,3]=56.06; osltr[19,4]=60.16;
osltr[20,1]=71.86; osltr[20,2]=76.07; osltr[20,3]=80.06; osltr[20,4]=84.45;
osltr[21,1]=97.18; osltr[21,2]=102.14;osltr[21,3]=106.74;osltr[21,4]=111.01;
osltr[22,1]=126.58;osltr[22,2]=131.70;osltr[22,3]=136.49;osltr[22,4]=143.09;
osltr[23,1]=159.48;osltr[23,2]=165.58;osltr[23,3]=171.28;osltr[23,4]=177.20;
osltr[24,1]=196.37;osltr[24,2]=202.92;osltr[24,3]=208.81;osltr[24,4]=215.54;
osltr[25,1]=236.54;osltr[25,2]=244.15;osltr[25,3]=251.30;osltr[25,4]=257.68;
osltr[26,1]=282.45;osltr[26,2]=291.40;osltr[26,3]=298.31;osltr[26,4]=307.64;
osltr[27,1]=332.00;osltr[27,2]=341.00;osltr[27,3]=348.00;osltr[27,4]=359.00;
osltr[28,1]=384.00;osltr[28,2]=395.00;osltr[28,3]=403.00;osltr[28,4]=415.00;
osltr[29,1]=440.00;osltr[29,2]=454.00;osltr[29,3]=462.00;osltr[29,4]=476.00;
osltr[30,1]=501.00;osltr[30,2]=517.00;osltr[30,3]=525.00;osltr[30,4]=542.00;

/* QbarG - add 30 */
osltr[31,1]= 2.69; osltr[31,2]= 3.76; osltr[31,3]= 4.95; osltr[31,4]= 6.65;
osltr[32,1]=13.33; osltr[32,2]=15.41; osltr[32,3]=17.52; osltr[32,4]=20.04;
osltr[33,1]=26.79; osltr[33,2]=29.68; osltr[33,3]=32.56; osltr[33,4]=35.65;
osltr[34,1]=43.95; osltr[34,2]=47.21; osltr[34,3]=50.35; osltr[34,4]=54.46;
osltr[35,1]=64.84; osltr[35,2]=68.52; osltr[35,3]=71.80; osltr[35,4]=76.07;
osltr[36,1]=89.48; osltr[36,2]=94.15; osltr[36,3]=98.33; osltr[36,4]=103.18;
osltr[37,1]=118.50;osltr[37,2]=124.24;osltr[37,3]=128.45;osltr[37,4]=133.57;
osltr[38,1]=150.53;osltr[38,2]=156.00;osltr[38,3]=161.32;osltr[38,4]=168.36;
osltr[39,1]=186.39;osltr[39,2]=192.89;osltr[39,3]=198.82;osltr[39,4]=204.95;
osltr[40,1]=225.85;osltr[40,2]=233.13;osltr[40,3]=239.46;osltr[40,4]=247.18;
osltr[41,1]=269.96;osltr[41,2]=277.71;osltr[41,3]=284.87;osltr[41,4]=293.44;
osltr[42,1]=317.00;osltr[42,2]=326.00;osltr[42,3]=334.00;osltr[42,4]=344.00;
osltr[43,1]=369.00;osltr[43,2]=378.00;osltr[43,3]=388.00;osltr[43,4]=399.00;
osltr[44,1]=424.00;osltr[44,2]=435.00;osltr[44,3]=445.00;osltr[44,4]=458.00;
osltr[45,1]=484.00;osltr[45,2]=495.00;osltr[45,3]=507.00;osltr[45,4]=522.00;
endp;

proc(0)=genchi;
chisq[1,1]=2.71; chisq[1,2]=3.84; chisq[1,3]=5.02; chisq[1,4]=6.63;
chisq[2,1]=4.61; chisq[2,2]=5.99; chisq[2,3]=7.38; chisq[2,4]=9.21;
chisq[3,1]=6.25; chisq[3,2]=7.81; chisq[3,3]=9.35; chisq[3,4]=11.3;
chisq[4,1]=7.78; chisq[4,2]=9.49; chisq[4,3]=11.1; chisq[4,4]=13.3;
chisq[5,1]=9.24; chisq[5,2]=11.1; chisq[5,3]=12.8; chisq[5,4]=15.1;
chisq[6,1]=10.6; chisq[6,2]=12.6; chisq[6,3]=14.4; chisq[6,4]=16.8;
chisq[7,1]=12.0; chisq[7,2]=14.1; chisq[7,3]=16.0; chisq[7,4]=18.5;
chisq[8,1]=13.4; chisq[8,2]=15.5; chisq[8,3]=17.5; chisq[8,4]=20.1;
chisq[9,1]=14.7; chisq[9,2]=16.9; chisq[9,3]=19.0; chisq[9,4]=21.7;
chisq[10,1]=16.0; chisq[10,2]=18.3; chisq[10,3]=20.5; chisq[10,4]=23.2;
chisq[11,1]=17.3; chisq[11,2]=19.7; chisq[11,3]=21.9; chisq[11,4]=24.7;
chisq[12,1]=18.5; chisq[12,2]=21.0; chisq[12,3]=23.3; chisq[12,4]=26.2;
chisq[13,1]=19.8; chisq[13,2]=22.4; chisq[13,3]=24.7; chisq[13,4]=27.7;
chisq[14,1]=21.1; chisq[14,2]=23.7; chisq[14,3]=26.1; chisq[14,4]=29.1;
chisq[15,1]=22.3; chisq[15,2]=25.0; chisq[15,3]=27.5; chisq[15,4]=30.6;
endp;

proc(0)=genkpss;
kpcrit[1,1]=0.347; kpcrit[1,2]=0.463; kpcrit[1,3]=0.574; kpcrit[1,4]=0.739;
kpcrit[2,1]=0.119; kpcrit[2,2]=0.146; kpcrit[2,3]=0.176; kpcrit[2,4]=0.216;
endp;

proc(0)=gendfcrit;
dfcrit[1,1]=-3.4335-5.999/n-29.25/n^2;
dfcrit[1,2]=-2.8621-2.738/n-8.36/n^2;
dfcrit[1,3]=-2.5671-1.438/n-4.48/n^2;
dfcrit[1,4]=-3.9638-8.353/n-47.44/n^2;
dfcrit[1,5]=-3.4126-4.039/n-17.83/n^2;
dfcrit[1,6]=-3.1279-2.418/n-7.58/n^2;
endp;

end;

