/* 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.              */

/*****************************************************************************/
/***  Procedure MGKS to perform ADF test on individual series to test for   **/
/***     unit root. Specifies lags optimally according to AIC, SC.          **/
/***  (C) 1996 By Michele Gambera and Kristen Strellec                      **/
/***                                                                        **/
/***  FORMAT:  DXt=a+bt+rXt-1+g1DXt-1+....+ gkDXt-k +et                     **/
/***     x is a column vector, the series to be tested for unit root        **/
/***     mxlg is the maximum number of lags included in the equation        **/
/***     costante=0 for no constant, 1 for constant, 2 for constant & trend **/
/***  USE: {akaike,schwarz}=mgks(x,mxlg,costante)                           **/
/*****************************************************************************/

Proc (2) = MGKS(x,mxlg,costante);

local deltax,aka,sc,conto,ttest,ma,bta,sebta,siga,a,littler,ftest,pro;
local lin,test,RT,Maik,Msc,Maikr,Mscr,testlags,z,dep,arpar,k;
test={};
testlags={};
_olsres=1;
__output=1;
__altnam=0;

deltax=trimr(x-lag1(x),1,0);
/**  Creates Xt-Xt-1; returns (n-1) observations  **/
conto=1;
if costante==0;
        __con=0;           /** no constant term added **/
        else;
        __con=1;           /** constant term added **/
endif;

Do while conto<=mxlg;
        output off;

/***  Create matrix of Xt-1, deltax(t-1),..., deltax(t-k)   ***/
/***  Adjust for missing observations                       ***/
        k=0;
        z=trimr(lag1(x),1,0);
                if (mxlg>0);
                        do until k>=conto;
                        k=k+1;
                        z=z~lagn(deltax,k);
                        endo;
                endif;
        dep=trimr(deltax,k,0);
        z=trimr(z,k,0);
        pro=seqa(1,1,rows(dep));
        if costante<2;
          {a,ma,bta,a,a,sebta,siga,a,a,a,a}=OLS(0,dep,z);
        else;
          {a,ma,bta,a,a,sebta,siga,a,a,a,a}=OLS(0,dep,z[.,1]~pro~z[.,2:cols(z)])
;
        endif;
        output on;
        if costante==0;
          ttest=bta[1]/sebta[1];
          arpar=bta[1];
        else;
          ttest=bta[2]/sebta[2];
          arpar=bta[2];
/*****  ADF F-test (Wald)  ****/
          if costante==2;
            littler={0,0};
            bta=bta[2:3];
            ma=invpd(ma[1:rows(ma)-1,1:rows(ma)-1]);
            ma=ma[2:3,2:3];
            ftest=(bta-littler)'invpd((siga^2)*ma)*(bta-littler)/2;
          endif;
        endif;

/****  AKAIKE & SCHWARZ CRITERIA FOR LAG SPECIFICATION                  ****/
/****  AIC(k)=ln(sigma^2) + 2k/T  where sigma^2 is the ML estimator     ****/
/****                             k=lag length and T is #obs            ****/
/****  SC(k)=ln(sigma^2) + klnT/T  where sigma^2 is the ML estimator    ****/
/****                             k=lag length and T is #obs            ****/
/****    Adjustments made for degree of freedom for  OLS estimator of   ****/
/****                   sigma^2                                         ****/

aka=ln((siga^2)/rows(dep)*(rows(dep)-(conto+costante)))+2*conto/rows(dep);
sc=ln((siga^2)/rows(dep)*(rows(dep)-(conto+costante)))+conto*ln(rows(dep))/
rows(dep);
         if costante==2;
          lin=conto~aka~sc~ttest~ftest;
        else;
          lin=conto~aka~sc~arpar~ttest;
        endif;
test=test|lin;
conto=conto+1;
endo;
testlags=testlags|test;

RT=MININDC(test);
 @ returns column vector containing index of smallest element in each column @
Maik=RT[2];
Msc=RT[3];
MAIKR=test[Maik,.];
MSCR=test[Msc,.];
"";
"Number of observations" rows(dep);
if costante==2;
    "Schwarz Optimal ADF";"Lags, AIC, SC, t-test, F test";
else;
    "Schwarz Optimal ADF";"Lags, AIC, SC, (1-rho), t-test";
endif;
mscr;"";
retp(maikr,mscr);

endp;
