```
# This is the top-level user-interface function.
# It calls BCtest_fromto()
# which in turn calls caus_freq() to calculate the actual test for a given freq,
# and then puts all test results in a matrix.
# At the end this function graphs the results.
# some checks:
if sum(missing(target)) > 0
funcerr "Missing values are not supported (found in target)."
endif
if sum(missing(cause)) > 0
funcerr "Missing values are not supported (found in cause)."
endif
if nelem(condvars) > 0
if sum(ok(condvars)-1) < 0
# missing() doesn't work for lists; the -1 seems to broadcast fine
# the lower-than-term is correct because everything is inverted here
funcerr "Missing values are not supported (found in conditioning)."
endif
endif
# construct the matrix input
matrix mY = {target,cause,condvars}
T = rows(mY)
numofvars = cols(mY)
# further checks
if rank(mY) < numofvars
funcerr "Exact linear dependence found in input (duplicated vars?)."
endif
if lagorder >= T
funcerr "Can't have more lags than observations"
endif
if (lagorder*numofvars) >= T
funcerr "No degrees of freedom left (too many lags, too many variables)"
endif
if freqpoints >= $nobs
printf "Warning: too many tested frequencies for graphing in current workfile."
# Because only series can be graphed, and longer series exceed workfile sample.
printf " (But test results will still be stored in return matrix.)\n"
do_graph = 0
else
do_graph = 1
endif
# now do the actual test
matrix BCtestoutput = BCtest_fromto(mY,lagorder,0.01,3.14,freqpoints) # to be returned
# add the critical value as constant series for convenience
numoffreqs = rows(BCtestoutput)
matrix BCtestoutput = BCtestoutput ~ (critical(X,2,siglevel)*ones(numoffreqs,1))
# graph the results
if do_graph = 1
startobs = $t1
# debug: print startobs
endobs = startobs + numoffreqs - 1
smpl startobs endobs # hope it isn't necessary to undo this later,
# because it only affects the inner scope?
series BCfreqs = BCtestoutput[,1]
series BCstats = BCtestoutput[,2]
series BCcvals = BCtestoutput[,3]
graph BCcvals BCstats BCfreqs
endif
return BCtestoutput
```

```
# this function calculates the actual test statistic
# mY: first col target, second col cause, further cols conditioning vars
T = rows(mY)
K = cols(mY)
matrix mXstar = mY[3:T,2] - 2*cos(freqpoint)*mY[2:T-1,2] + mY[1:T-2,2]
# (the Gegenbauer polynomial for non-causality at frequency freqpoint)
matrix mX = mY[lagorder:T-1,] ~ mY[lagorder-1:T-2,] # two lags
# (at least lagorder 2 is required because otherwise no cyclical behavior)
matrix mind = {2, K+2} # index vector refers to first and second lag of causal var
upto = lagorder-2
loop for i=1..upto --quiet # add more lags
matrix mtemp = mY[lagorder-1-i:T-2-i,] # next lag of all variables
matrix mX = (K>2) ? (mX ~ mtemp[,1] ~ mtemp[,3:K]) : (mX ~ mtemp[,1])
# (remove causing var, handled below; K>2 means we have conditioning vars)
matrix mX = mX ~ mXstar[lagorder-1-i:T-2-i]
# (these Gegenbauer lags of the causal variable are non-causal by construction)
# (it seems un-important that the Xstarlags come last, so they don't, here)
endloop
matrix mX = mX ~ ones(T-lagorder,1)
matrix mDepvar = mY[lagorder+1:T,1]
# here the original Gauss code was b=depvar/x
# AFAIK this means regress depvar on x and get the estimates
matrix mU = {} # to hold the residuals
matrix mb = mols(mDepvar,mX,&mU)
scalar sig2 = mU'mU
scalar sig2 = sig2 / (T-lagorder-cols(mX))
matrix mVarb = sig2 * invpd(mX'mX)
# now the vector (2;k+2) was used for indexing, I hope this is what was meant
scalar wald = qform(transp(mb[mind]), invpd(mVarb[mind,mind]))
return wald
```

```
# We expect the return matrix to have numoffreqs rows, but could be slightly different
# so one shouldn't rely on that.
# The idea of having fromfreq and tofreq is that in a batch processing context
# we can save computing time by looking only at the freqs we want.
# For a general graph instead we would look at the entire frequency band.
matrix mcontainer = zeros(0,2) # new w/o the critical value; 1st col freqs, 2nd col stats
T = rows(mY)
step = (tofreq-fromfreq)/(numoffreqs-1)
loop for (f=fromfreq; f<=tofreq; f+=step) --quiet
teststat = caus_freq(mY,lagorder,f)
matrix mcontainer = mcontainer | {f, teststat}
endloop
return mcontainer
```