/*                        dickeyful.gpg
   does dickey-fuller, new gene style.
   regression is Dy(t) on rho y(t-1) + Dy(t-1) + .. + e(t)
   tabulate distribution of rho under null rho is zero, bootstrap e.
*/
output file=dickeyful.out reset;
getdata = 1;
#include getdat.g;
dataset = 1;

if dataset == 1; freq = 1; stpnt = 47; endif;
if dataset == 2; freq = 1; stpnt = 47; endif;
if dataset == 3; freq = 1; stpnt = 47; endif;
if dataset == 4; freq = 4; stpnt = 27; endif;
if dataset == 5; freq = 4; stpnt = 47; endif;

if getdata;
    call getdat(stpnt,dataset);
endif;
#include datsetup.gpg


ci = c-y;
T = rows(ci);
dci = ci[2:T]-ci[1:T-1];

proc(3) = dickful(ci,lag,prnt);
    clearg bigb,bigt,j,lhv,rhv,i,bci,eci,sigci,sebci,tci,lag;
    if prnt;
        "";
        format /RD 9,3;
        "Augmented Dickey fuller regressions of c-y";
        "D(c-y)(t) = const. + rho(c-y)(t-1) + b1 D (C-y) (t-2) + ... ";
     endif;

    j = lag[1];
    do while j <= lag[2];  /* max amount of augmenting */
             " ";
            lhv = dci[1+j:T-1];
            rhv = ones(T-1-j,1)~ci[1:T-1-j];
            if j > 0;
                i = 1;
                do while i <= j;
                    rhv = rhv~dci[1+j-i:T-i-1];
                i = i+1;
                endo;
            endif;
            bci = lhv/rhv;
            eci = lhv-rhv*bci;
            sigci = eci'eci/(T-j-1);
            sebci = (diag(inv(rhv'rhv))*diag(sigci)')^(1/2);
            tci = bci./sebci;
            if prnt;
                "coeff." bci';
                "se    " sebci';
                "tstat " tci';
            endif;
            j = j+1;
        endo;
        retp(bci,tci,eci);
endp;

proc(3) = nocivec(ci,lag,prnt);
    clearg bigb,bigt,j,lhv,rhv,i,bci,eci,sigci,sebci,tci,lag;
    if prnt;
        "";
        format /RD 9,3;
        "Augmented Dickey fuller regressions of c-y";
        "D(c-y)(t) = const. + rho(c-y)(t-1) + b1 D (C-y) (t-2) + ... ";
     endif;

    j = lag[1];
    do while j <= lag[2];  /* max amount of augmenting */
             " ";
            lhv = dci[1+j:T-1];
            rhv = ones(T-1-j,1)~ci[1:T-1-j];
            if j > 0;
                i = 1;
                do while i <= j;
                    rhv = rhv~dci[1+j-i:T-i-1];
                i = i+1;
                endo;
            endif;
            bci = lhv/rhv;
            eci = lhv-rhv*bci;
            sigci = eci'eci/(T-j-1);
            sebci = (diag(inv(rhv'rhv))*diag(sigci)')^(1/2);
            tci = bci./sebci;
            if prnt;
                "coeff." bci';
                "se    " sebci';
                "tstat " tci';
            endif;
            j = j+1;
        endo;
        retp(bci,tci,eci);
endp;



{bci,tci,eci} = dickful(ci,2|2,1);

output off;