           Unit Root Tests (Isaac and Rapach,May96)

/********************************************************************
*********************************************************************
UNIT ROOT TESTING: a collection of procedures
Authors: Alan G. Isaac and David Rapach
==> For public, non-commercial use only
==> Please include proper attribution should this code be used in
    a research project or in other code
==> Authors make no performance guarantees
==> Program written for GAUSS Version 3.1.13
    (vendor: Aptech Systems, Maple Valley, WA)
*********************************************************************
********************************************************************/

/* Procedure Definition Area */

/********************************************************************
Proc            : ADF
Authors         : Alan G. Isaac and David Rapach
Last revised    : 27 May 1996
Input           : y -- (N x 1) vector, the time series of interest
                  p -- scalar, maximum lag (max autocorr order)
Output          : Prints augmented Dickey-Fuller t-stats
Globals         : None
Notes           : Requires procs `varlags()' and `adf_ls()'(provided)
References      : Russell Davidson and James G. MacKinnon,
                      _Estimation and Inference in Econometrics_
                      (New York: Oxford, 1993)
                      See references therein for original sources
                      See p.708, table 20.1 for critical values
                  Harris, 1992, Economics Letters
********************************************************************/
proc(0)=adf(y,p);
    local obs,i,b,ty,temp,oldcv,oldnv;
    obs=rows(y);
    format 8,4;
    "One-sided test of H0: Unit root vs. H1: Stationary";?;
    "Approximate asymptotic critical values (t-ratio):";
    "------------------------------------------------------------";
    "  1%      5%      10%      Model";
    "------------------------------------------------------------";
    "-2.56   -1.94   -1.62";;
    "     Simple ADF (no constant or trend)";
    "-3.43   -2.86   -2.57";;
    "     ADF with constant (no trend)";
    "-3.96   -3.41   -3.13";;
    "     ADF with constant & trend";
    "------------------------------------------------------------";
"You began with a series of "$+ftocv(obs,1,0)$+" observations, in which case";
"Harris (1992) recommends using "$+ftocv(int(12*(obs/100)^(1/4)),1,0)$+" lags";
"
------------------------------------------------------------
                     ADF coef on lagged level
                             (t-ratio)
                For model including the following:
                ----------------------------------
         Lags   No Constant  Constant    Trend           Obs
                No Trend     No Trend    Constant
------------------------------------------------------------";
    i=0;
        oldcv=__fmtcv;oldnv=__fmtnv;
    do while i<=p;
        {b,ty,obs}=adf_ls(y,i,"off");
        __fmtcv="*.*s"~12~8;__fmtnv="*.*lg"~12~3;
        temp=printfmt(i~b[1,.]~obs,1);print;
        temp=printfmt("  "~(0$+"("$+ftocv(ty[1,.],1,2)$+")")~"  ",0);print;
        i=i+1;
    endo;
"------------------------------------------------------------";
        __fmtcv=oldcv;__fmtnv=oldnv;
retp;
endp;

/********************************************************************
Proc            : ADF_LS
Authors         : Alan G. Isaac and David Rapach
Last revised    : 27 May 1996
Input           : y   -- (N x 1) vector, the time series of interest
                  p   -- scalar, maximum lag (max autocorr order)
                  prt -- string, prt="on" to print results
Output          : Prints augmented Dickey-Fuller t-stats
                  b   -- coefs from ADF regression
                  t   -- t-stats from ADF regression
                  obs -- scalar, rows(y)-(p+1)
Globals         : None
Notes           : Requires procs `varlags()' (provided)
References      : Russell Davidson and James G. MacKinnon,
                      _Estimation and Inference in Econometrics_
                      (New York: Oxford, 1993)
                      See references therein for original sources
********************************************************************/
proc(3)=adf_ls(y,p,prt);
    local labelx,addx,y_1,dy,dylags,temp,report,rhs,gap,oldcv,oldnv,
          T,b,e,SSE,sig2,covb,seb,ty,i;
    T=rows(y);                     @ number of observations @
    labelx="x(-1)";
    addx={};
    {y,y_1}=varlags(y,1);          @ generating lags @
    dy=y-y_1;                      @ differencing @
    if p /= 0;
      {dy,dylags}=varlags(dy,p);   @ generating lags @
      addx=addx~dylags;
      labelx=labelx|(0$+"dx(-"$+ftocv(seqa(1,1,p),1,0)$+")");
    endif;
    labelx=labelx|"constant"|"trend"|"RSS";
    T=rows(dy);                    @ # of usable obs @
    y_1=y_1[p+1:p+T];              @ usable lagged y @
    addx=y_1~addx;
    gosub regress(addx,miss(zeros(2,1),0));
      pop temp;
      report=temp;
    addx=addx~ones(T,1);           @ add constant to regressors @
    gosub regress(addx,miss(0,0)); @ add constant to regressors @
      pop temp;
      report=report~temp;
    addx=addx~seqa(-int(T/2),1,T); @ add trend to regressors @
    temp={};
    gosub regress(addx,temp);
      pop temp;
      report=report~temp;
    if prt $== "on";
      print;
      "ADF results for "$+ftocv(p,1,0)$+" lags";
"------------------------------------------------------------
                      ADF regression coefficients
                             (t-ratio)
                For model including the following:
                ----------------------------------
                No Constant  Constant    Trend
     RHS var    No Trend     No Trend    Constant
--------------------------------------------------";
    i=1;
    oldcv=__fmtcv;oldnv=__fmtnv;
    __fmtcv="*.*s"~12~8;__fmtnv="*.*lg"~12~3;
    temp=printfmt(labelx[1]~report[1,1|3|5],0~1~1~1);print;
    do while i<=p+3;
      temp=printfmt("  "~(0$+"("$+ftocv(report[i,2|4|6],1,2)$+")")~"  ",0);
      print;
      i=i+1;
      temp=printfmt(labelx[i]~report[i,1|3|5],0~1~1~1);print;
    endo;
"------------------------------------------------------------";
"Note: Final line is sum of squared residuals";
/* Could use these for Dickey Fuller (1981) Phi tests */
        __fmtcv=oldcv;__fmtnv=oldnv;
    endif;
retp(report[.,1|3|5],report[.,2|4|6],T); @ coef, t-stat, usable obs @
regress:
    pop gap;
    pop rhs;
    b=dy/rhs;                     @ OLS estimates @
    e=dy-rhs*b;                   @ resids @
    SSE=e'e;                      @ sum of squared resids @
    sig2=SSE/(T-cols(rhs));       @ e2 estimate @
    covb=sig2*invpd(rhs'rhs);     @ var-cov(b) @
    seb=sqrt(diag(covb));         @ std errs @
    return((b|gap|SSE)~((b./seb)|gap|miss(0,0)));
endp;

/**********************  PROC VARLAGS  *****************************
**   author: Alan G. Isaac
**   last update: 5 Dec 95      previous: 15 June 94
**   FORMAT
**          { x,xlags } = varlags(var,lags)
**   INPUT
**        var  - T x K matrix
**        lags - scalar, number of lags of var (a positive integer)
**   OUTPUT
**          x -     (T - lags) x K matrix, the last T-lags rows of var
**          xlags - (T - lags) x lags*cols(var) matrix,
**                  being the 1st through lags-th
**                  values of var corresponding to the values in x
**                  i.e, the appropriate rows of x(-1)~x(-2)~etc.
**   GLOBAL VARIABLES: none
**********************************************************************/
proc(2)=varlags(var,lags);
    local xlags;
    xlags = shiftr((ones(1,lags) .*. var)',seqa(1-lags,1,lags)
                                            .*. ones(cols(var),1),miss(0,0))';
    retp(trimr(var,lags,0),trimr(xlags,0,lags));
endp;

/***************************************************************************
Proc            : KPSS
Author          : David Rapach
Last revised    : 27 May 1996
Input           : y -- (T x 1) vector of obs on variable of interest
                : p -- specified max order of autocorrelation
Output          : Prints KPSS statistics
Globals         : None
Reference       : Denis Kwiatowski, Peter C.B. Phillips, Peter Schmidt,
                      and Tongcheol Shin, "Testing the null hypothesis of
                      stationarity against the alternative of a unit root:
                      How sure are we that economic series have a unit
                      root?" _Journal of Econometrics_ 54 (Oct./Dec.
                      1992): 159-178
Notes           : KPSS(trend) based on y(t)=a+bt+e(t) -- trend stationarity
                  KPSS(level) based on y(t)=a+e(t)    -- level stationarity
***************************************************************************/
proc(0)=kpss(y,l);
    local T,trend,X,b,e,g0,p,g,i,w,s2,S,ei,etathat,etamhat,lag,
          eta,oldnv;

    /* Level-stationary test */

    T=rows(y);              @ # of obs @
    X=ones(T,1);            @ RHS var @
    b=y/X;                  @ OLS estimates @
    e=y-X*b;                @ OLS resids @
    g0=e'e/T;               @ gamma_0=(1/T)sum[e(t)^2] @
    etamhat=zeros(l+1,1);
    p=1;
    do until p>l;
        g=zeros(p,1);       @ autocovs @
        i=1;
        do until i>p;
            g[i,1]=(e[1+i:T,1]'e[1:T-i,1])/T;@ gamma_i=(1/T)sum[e(t)e(t-i)] @
            i=i+1;
        endo;
        w=zeros(p,1);
        i=1;
        do until i>p;
            w[i,1]=(p+1-i)/(p+1);   @ Bartlett window weight @
            i=i+1;
        endo;
        s2=g0+2*w'g;                @ consistent error variance estimate @
        S=zeros(T,1);               @ resid partial sum process @
        i=1;
        do until i>T;
            ei=e[1:i,.];
            S[i,.]=sumc(ei);        @ S(i)=sum[e(i)] @
            i=i+1;
        endo;
        etamhat[1,.]=(1/(g0*T^2))*S'S;      @ KPSS eqn (13), l=0 @
        etamhat[p+1,.]=(1/(s2*T^2))*S'S;    @ KPSS eqn (13) @
        p=p+1;
    endo;

    /* Trend-stationary test */

    trend=seqa(1,1,T);      @ linear time trend @
    X=ones(T,1)~trend;      @ regressor matrix @
    b=y/X;                  @ OLS estimates @
    e=y-X*b;                @ OLS resids @
    g0=e'e/T;               @ gamma_0=(1/T)sum[e(t)^2] @
    etathat=zeros(l+1,1);
    p=1;
    do until p>l;
        g=zeros(p,1);       @ autocovs @
        i=1;
        do until i>p;
            g[i,1]=(e[1+i:T,1]'e[1:T-i,1])/T;@ gamma_i=(1/T)sum[e(t)e(t-i)] @
            i=i+1;
        endo;
        w=zeros(p,1);
        i=1;
        do until i>p;
            w[i,1]=(p+1-i)/(p+1);   @ Bartlett window weight @
            i=i+1;
        endo;
        s2=g0+2*w'g;                @ consistent error variance estimate @
        S=zeros(T,1);               @ resid partial sum process @
        i=1;
        do until i>T;
            ei=e[1:i,.];
            S[i,.]=sumc(ei);        @ S(i)=sum[e(i)] @
            i=i+1;
        endo;
        etathat[1,.]=(1/(g0*T^2))*S'S;      @ KPSS eqn (17), l=0 @
        etathat[p+1,.]=(1/(s2*T^2))*S'S;    @ KPSS eqn (17) @
        p=p+1;
    endo;
    format 8,4;
    "One-sided test of H0: Stationary vs. H1: Unit root";?;
    "Approximate asymptotic critical values";
    "--------------------------------------------------";
    "              10%     5%      1%";
    "--------------------------------------------------";
    "Level        0.347   0.463   0.739";
    "Trend        0.119   0.146   0.216";
    "--------------------------------------------------";?;
    lag=seqa(0,1,l+1);
    eta=etamhat~etathat;
    oldnv=__fmtnv;
    __fmtnv="*.*lg"~8~3;
    "KPSS stats";
    "--------------------------------------------------";
    "      Lag     Level    Trend";
    "--------------------------------------------------";
    i=1;
    do until i>l+1;
        call printfmt(lag[i,.],1);"    " eta[i,.];
        i=i+1;
    endo;
    "--------------------------------------------------";
    "Obs";;call printfmt(T,1);
    __fmtnv=oldnv;
endp;
