nulldata 100

# seed so that experiment can be repeated exactly
set seed 89675430

# params of error process: can set some to zero, but
# here we do AR(3)
scalar rho1 = 0.5
scalar rho2 = -0.3
scalar rho3 = 0.2

# y params: d = drift, a = autoregressive parameter
scalar d = 0.02
scalar a = 0.9

# test down from
scalar maxlag = 6

# test-recoding scalars
scalar zval, pz, zlag, pzstar

# replications
scalar replics = 1000

# matrix to record results across relications
matrix lsel = zeros(replics, 4)

# outer loop: replications
loop j=1..replics -q

    # build stationary AR(p) error process for u
    series u = 0
    series e = normal()
    u = rho1*u(-1) + rho2*u(-2) + rho3*u(-3) + e

    # construct y series for "ADF testing"
    series y = 0
    series y = d + a*y(-1) + u

    # difference for DF regression
    dy = diff(y)
    lags maxlag ; dy
    # comprehensive list of lagged differences
    list dyL = dy_*

    # Information Criterion matrix
    matrix IC = zeros(3, maxlag+1)
    zlag = 0

    # inner loop: test for lag length
    loop for (i=maxlag; i>=0; i--) -q
        ols dy 0 y(-1) dyL -q
        IC[1,i+1] = $aic
        IC[2,i+1] = $bic
        IC[3,i+1] = $hqc
        if (i > 0 && zlag == 0) 
            # we haven't yet found a significant "final lag"
            zval = abs($coeff(dy_$i)/$stderr(dy_$i))
            pz = 2 * pvalue(z, zval)
            # printf "pz = %g\n", pz
            if (pz <= 0.10) 
                zlag = i
                pzstar = pz
            endif
        endif
        if i > 0
            # discard the last lag
            dyL -= dy_$i
        endif
    endloop

    # just for checking...
    # print IC

    matrix iCmin = iminr(IC) - 1

    # Verbose printing to check results (set replics = 1)
    # matrix Cmin = minr(IC)
    # printf "AIC min = %8.3f at lag %d\n", Cmin[1], iCmin[1]
    # printf "BIC min = %8.3f at lag %d\n", Cmin[2], iCmin[2]
    # printf "HQC min = %8.3f at lag %d\n", Cmin[3], iCmin[3]
    # if zlag > 0
    #   printf "z p-val = %8.3f at lag %d\n", pzstar, zlag
    # else
    #   print "z-test selects no lags"
    # endif

    lsel[j,1:3] = iCmin'
    lsel[j,4] = zlag

endloop

# Uncomment to print the full results
# print lsel
matrix stats = meanc(lsel) | sdc(lsel)
colnames(stats, "AIC BIC HQC z-test")
printf "\n%d replications, test down from %d lags\n", replics, maxlag
print "Lag order selected: mean on row 1, s.d. on row 2"
printf "\n%8.2f\n", stats