Partial implementation of Brüggemann/Lütkepohl (2000)

[macro]

I implemented a version of one of the model selection algorithms in Brüggemann and Lütkepohl (2000) (see the "Testing Procedure" beginning on p4), and I'm hoping for feedback on its correctness and programming style and to share it with anyone who might find it helpful. There are two variations on the equation-by-equation algorithm included in my procedure; the "basic" method replicates the paper's Testing Procedure, while the "topdown" procedure uses a similar method, except that on a given iteration, the kth lag of a regressor can only be dropped if the (k+1)th lag has already been dropped.

Both use a threshold value that is a function of the current step of the procedure. The procedure also allows the user to force certain variables to stay in the regression through a second supplementary card. In my example code below, this is used to ensure the constant remains in the regression.

Code: Select all

procedure lutkephol dep_var begin_samp end_samp criterion method print
    type series dep_var
    type integer begin_samp end_samp
    option choice criterion 1 aic sbc hq
    option choice method 1 basic topdown
    option switch print 0

    local equation exog_vars
    equation exog_vars *

    local equation exog_vars_keep
    equation exog_vars_keep *

    if criterion == 1 {
        frml cT = 2.0
    }
    else if criterion == 2 {
        frml cT = 2.0 * log(t)
    }
    else if criterion == 3 {
        frml cT = 2.0 * log(log(t))
    }
    report(action = define, use = lutreport, hlabels = ||"Iteration", "Regressor", "t-stat", "threshold"||)

    compute tbl_keep = %eqntable(exog_vars_keep)
    compute rl = %rlfromeqn(exog_vars)
    cmoment(noprint) begin_samp end_samp
    # dep_var rl

    compute min_tstat = %na
    compute critical_value = %na
    compute iter = 0
    until min_tstat > critical_value {
        compute iter = iter + 1
        linreg(cmom, noprint, define = eq) dep_var begin_samp end_samp
        # rl

        compute T = %nobs
        compute K = %nreg
        compute tbl = %eqntable(eq)

        * determine which regressors can be dropped
        declare vector[integer] free_to_drop(K)

        compute idx_free_to_drop = %tablefindmiss(tbl, tbl_keep)
        dofor i = idx_free_to_drop
            compute free_to_drop(i) = 1
        end dofor idx_free_to_drop

        if method == 2 {
            * topdown
            do i = 1, K
                * if we've already determined that the variable can't be
                * dropped because it was passed in the second card, skip
                * to the next variable
                if free_to_drop(i) == 0 {
                    next
                }
                compute s = tbl(1, i)
                compute lag = tbl(2, i)
                compute idx_of_next_lag = %tablefind(tbl, s, lag + 1)

                * if the regression contains the next lag of series s,
                * then the current lag cannot be dropped
                if idx_of_next_lag > 0 {
                    compute free_to_drop(i) = 0
                }
            end do i
        }
        ewise %tstats(i) = %if(free_to_drop(i) == 0, %na, abs(%tstats(i)))
        compute idx_of_min_tstat = %minindex(%tstats)
        compute min_tstat = %tstats(idx_of_min_tstat)
        compute critical_value = sqrt((exp(cT(T)/T) - 1)*(T-K+iter-1))

        if min_tstat < critical_value {
            compute var_to_drop = %tableextract(tbl, ||idx_of_min_tstat||)
            compute new_tbl = %tableextract(tbl, %tablefindmiss(tbl, var_to_drop))
            compute rl = %rlfromtable(new_tbl)
            compute regressors = %eqnreglabels(eq)
            report(action = modify, use = lutreport, row = new, atcol = 1) iter regressors(idx_of_min_tstat) min_tstat critical_value
        }
    }
    report(action = modify, use = lutreport, row = new, atcol = 1) iter "-complete-" min_tstat critical_value
    report(action = format, use = lutreport, atcol = 1, tocol = 1, align = center)
    report(action = format, use = lutreport, atcol = 2, tocol = 2, align = left)
    report(action = format, use = lutreport, atcol = 3, tocol = 4, align = decimal, picture = "##.####")

    if print {
        report(action = show, use = lutreport)
    }

    compute %eq = eq
    compute %iterations = iter
end

This sample code demonstrates the procedure using the attached sample data.

Code: Select all

calendar(q) 1969 1

source(echo) "lutkephol.src"
compute max_lags = 8

open data "data.xls"
    data(format = xls, org = columns) /
close data

@lutkephol(criterion = aic, method = basic, print) lgdpk 1969:1 2006:4
# lgdpk{1 to max_lags} lcdk{1 to max_lags} ffedk{1 to max_lags} $
  lresinvk{1 to max_lags} lpotgdpk{1 to max_lags} constant
# constant

linreg(create, equation = %eq, print)

Does my code look correct? My main concern is that I've implemented the Testing Procedure incorrectly.

[TomDoan]

It looks fine. They wouldn't happen to have provided an example with public data? If you want to check, I would throw in some extra displays each time you have a "transaction". You can have those covered by a TRACE option (a SWITCH that defaults to 0) so they won't appear when it's run in production.

Stylistically,

procedure lutkephol dep_var begin_samp end_samp criterion method print

has (1) lutkepohl misspelled and (2) the criterion method and print are options and shouldn't be listed on the parameter line. (I'm actually surprised that there are no complaints about that for the names being used two ways).

I would also be careful about the use of T in the procedure. frml(any_integer) works fine, it doesn't have to be the global T variable.

[macro]

Thank you for the feedback; I wasn't able to find any examples using public data, although the empirical application in the paper states that it uses the data from Christiano, Eichenbaum, Evans (1996).

Thanks for catching that typo too. When you mention including a TRACE option, are you referring to something more than the report, triggered by the PRINT option, that displays the result of each iteration?

Incidentally, is there a procedure or function that would simplify code like this:

Code: Select all

if criterion == 1 {
    compute criterion_lbl = "aic"
}
else if criterion == 2 {
    compute criterion_lbl = "sbc"
}
else if criterion == 3 {
    compute criterion_lbl = "hq"
}

Since RATS translates the text passed to the option internally (presumably at compile time), it seems like it might be possible to simplify branching statements like the one above and automatically access the text that an option's number corresponds to.

[TomDoan]

Code: Select all

fixed vect[string] clabels(3)
aic
sbc
hq

later

compute criterion_lbl=clabels(criterion)

See, for instance, the use of lwindowtext in hadri.src.

[TomDoan]

Thanks for catching that typo too. When you mention including a TRACE option, are you referring to something more than the report, triggered by the PRINT option, that displays the result of each iteration?

Add

option switch trace 0

in the header

then in

if min_tstat < critical_value {
compute var_to_drop = %tableextract(tbl, ||idx_of_min_tstat||)
compute new_tbl = %tableextract(tbl, %tablefindmiss(tbl, var_to_drop))

add

if trace {
disp "Dropping" %l(var_to_drop(1,1))+"{"+var_to_drop(2,1)+"} with t-value" min_tstat
}