```
# 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.
matrix mY = {target,cause,condvars}
# some checks:
T = rows(mY)
numofvars = cols(mY)
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
smpl 1 numoffreqs # 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 = mX ~ mtemp[,1] ~ mtemp[,3:K]
# (remove causing var, handled below)
/* old code:
matrix mX = mX ~ mY[lagorder-1-i:T-2-i, 1]
if K>2
# why is this thing nested in the lagorder>2 cond?
# -> only because the 2nd col (causing var) is not covered here
matrix mX = mX ~ mY[lagorder-1-i:T-2-i, 3:K]
endif
*/
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
/* old code:
loop for i=1..upto --quiet
# (loop again because the mXstar things must come at the end)
matrix mX = mX ~ mXstar[lagorder-1-i:T-2-i]
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
```