```
/* 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.
Or for a general graph instead just look at everything from 0 to pi.
*/
matrix mcontainer = zeros(0,2) # 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)
mcontainer = mcontainer | {f, teststat}
endloop
# make sure that we really don't miss pi if tofreq==pi
if (tofreq > $pi-1e-3) && (mcontainer[rows(mcontainer), 1] < $pi-1e-3)
# is it not yet done? then do it now explicitly
teststat = caus_freq(mY, lagorder, $pi)
mcontainer = mcontainer | {$pi, teststat}
endif
return mcontainer
```

```
# This is a kind of "middleware" for the GUI-level function
# "BreitungCandelonTest", linking the user input to the backends.
#
# We assume the following minimal existing bundle members:
# 1) b.mtarg (matrix/vector) -- we work with matrices here
# 2) b.mcause (matrix/vector)
#
# Then there are the optional members:
# 3) b.lagorder (int)
# 4) b.freqpoints (int)
# 5) b.siglevel (scalar)
# 6) b.fromfreq (scalar)
# 7) b.tofreq (scalar)
# 8) b.mcvars (matrix); (coming from "condvars")
# this matrix type is a workaround for the fact that lists
# are not allowed in bundles.
# 9) b.verbosity (int)
# 10) b.doplot (bool)
#
# On exit the constructed matrix is put into the bundle as "b.BCoutput"
# The return value is 0 on success.
# apply default values
if !inbundle(b, "lagorder")
b.lagorder = 3
endif
if !inbundle(b, "freqpoints")
b.freqpoints = 100
endif
if !inbundle(b, "siglevel")
b.siglevel = 0.05
endif
if !inbundle(b, "fromfreq")
b.fromfreq = 0.01 # almost zero
endif
if !inbundle(b, "tofreq")
b.tofreq = 3.13 # almost pi
endif
if !inbundle(b, "mcvars")
matrix b.mcvars = {}
endif
if !inbundle(b, "verbosity")
b.verbosity = 1
endif
if !inbundle(b, "doplot")
b.doplot = 1
endif
if b.verbosity
printf "Chosen lag order: %d\n", b.lagorder
if b.lagorder == 2
print "Warning: The test with only 2 lags is not suitable for I(1) variables."
endif
endif
# construct the matrix input
matrix mY = vec(b.mtarg) ~ vec(b.mcause) ~ b.mcvars
T = rows(mY)
numofvars = cols(mY)
# further checks
if rank(mY) < numofvars
funcerr "Exact linear dependence found in input (duplicated vars?)."
endif
if b.lagorder >= T
funcerr "Can't have more lags than observations"
endif
if (b.lagorder * numofvars) >= T
funcerr "No degrees of freedom left (too many lags, too many variables)"
endif
if (b.fromfreq != 0.01) || (b.tofreq != 3.13)
findmin = 1
if b.verbosity
print "-------------------------------------------------------"
printf "Performing frequency-band Granger test in the band "
printf "from %f to %f.\n", b.fromfreq, b.tofreq
print "The relevant test statistic will be the minimal one."
printf "Calculating/showing all pointwise results in the chosen "
printf "band for convenience.\n"
print "-------------------------------------------------------"
endif
# The adjustment factor for the test stat:
adjust = critical(X, 2, b.siglevel) / critical(X, 1, b.siglevel)
else
findmin = 0
if b.verbosity
print "-------------------------------------------------------"
print "Performing Breitung-Candelon tests over the whole range (0, pi)."
print "The tests are only valid pointwise for a pre-specified frequency."
print "Calculating/showing all pointwise results for convenience."
print "-------------------------------------------------------"
endif
endif
# now do the actual test
matrix BCoutput = BCtest_fromto(mY, b.lagorder, b.fromfreq, b.tofreq, b.freqpoints)
# add the critical value as constant series for convenience
numoffreqs = rows(BCoutput)
BCoutput = BCoutput ~ (critical(X,2, b.siglevel) * ones(numoffreqs,1))
# treat special cases 0 and pi
if BCoutput[1, 1] == 0
if b.verbosity
print "Frequency 0 is included, needs special treatment;"
printf "The raw test statistic with d.o.f.=1 is %f\n", BCoutput[1,2]
endif
# adjust this test statistic
BCoutput[1, 2] *= adjust
endif
if BCoutput[numoffreqs, 1] > 3.13
if b.verbosity
print "Frequency pi is included, needs special treatment;"
printf "The raw test statistic with d.o.f.=1 is %f\n", BCoutput[numoffreqs, 2]
endif
# adjust this test statistic
BCoutput[numoffreqs, 2] *= adjust
endif
# graph the results
colnames(BCoutput, "frequency teststat crit.val.")
if b.doplot
plotBC(BCoutput)
endif
if b.verbosity
if findmin # a range was specified
scalar row = iminc(BCoutput[, 2])
printf "The min test stat (adjusted if needed) is %f,\n", BCoutput[row, 2]
printf " associated with frequency %f.\n", BCoutput[row, 1]
printf "The critical value at significance level %f is %f,\n", b.siglevel, BCoutput[row, 3]
string decision = BCoutput[row, 2] > BCoutput[row, 3] ? "" : " do not"
printf " and the test decision is therefore:%s reject non-causality.\n", decision
endif
printf "The columns in the return matrix (if any) hold\n"
printf "(1) frequencies, (2) test statistics,\n"
printf "(3) critical value (repeated) at %4.3f level.\n", b.siglevel
endif
# copy the result matrix
matrix b.BCoutput = BCoutput
return 0
```

```
# This is the top-level user-interface function.
# It calls the middleware FreqGranger, which then 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.
#
/* reverted: 14 feb 17: take into account the lag of x when checking for missings
#
# some checks:
list allwithlags = lags(lagorder, cause)
if nelem(condvars)
allwithlags += condvars
endif
allwithlags += target lags(lagorder, target)
if sum(missing(allwithlags))
print "Warning: missing values found, will try to remove obs."
smpl allwithlags --no-missing
endif
*/
if sum(missing(target)) > 0
print "Warning: missing values found in target, will try to remove obs."
smpl target --no-missing
endif
if sum(missing(cause)) > 0
print "Warning: missing values found in cause, will try to remove obs."
smpl cause --no-missing
endif
if nelem(condvars) > 0
if sum(missing(condvars)) # old: sum(ok(condvars)-1) < 0
# old: missing() doesn't work for lists; the -1 seems to broadcast fine
# old: the lower-than-term is correct because everything is inverted here
print "Warning: missing values found in conditioning, will try to remove obs."
smpl condvars --no-missing
endif
endif
printf "Sample size (including lags, T+p): %d; from obs %s to obs %s.\n", $nobs, obslabel($t1), obslabel($t2)
# populate the transfer bundle
bundle b
matrix b.mtarg = {target}
matrix b.mcause = {cause}
matrix b.mcvars = nelem(condvars) ? {condvars} : {}
b.verbosity = 1
b.doplot = 1
loop foreach s lagorder freqpoints siglevel fromfreq tofreq -q
b.$s = $s
endloop
# the call to the middleware
res = FreqGranger(&b)
if !res
return b.BCoutput
else
print "Uh-oh, something went wrong; returning empty matrix"
return {}
endif
```

```
if cols(BCoutput) < 3 || rows(BCoutput) < 1
funcerr "need a bigger matrix"
endif
gnuplot 2 3 1 --matrix=BCoutput --suppress-fitted --single-yaxis
```

```
# this function calculates the actual test statistic
# Data: first col target, second col cause, further cols conditioning vars
T = rows(Data)
K = cols(Data)
# define an index vector that leaves out the causal var:
matrix ix_noX = (K > 2) ? ( {1} ~ seq(3,K) ) : {1}
# Regressors: start with x(-1) for all cases
matrix RHS = mlag(Data[, 2], 1)
# include y(-1) to y(-p) as well as conditioning lags
RHS = RHS ~ mlag(Data[, ix_noX], seq(1, lagord))
# add the constant in all cases
RHS = RHS ~ ones(T, 1)
if fpoint == 0 || (fpoint >= $pi-1e-3 && fpoint <= $pi+1e-3)
# (special treatment needed for the corner cases)
matrix ix_test = {1} # which regressor to test (x(-1) here)
# and different filters for f==0 or f==pi
if fpoint == 0
# estimate a kind of ECM form
# dx(-1) to dx(-(p-1))
RHS = RHS ~ mlag(diff(Data[, 2]), seq(1, lagord-1))
else # so fpoint==pi
# estimate with (1+L) filter instead of (1-L) or 1-2cos(f)L+L^2
# (1+L)x(-1) to (1+L)x(-(p-1))
RHS = RHS ~ mlag(Data[, 2] + mlag(Data[, 2], 1), seq(1, lagord-1))
endif
else # interior standard case with 2 restrictions
# insert second lag of causal, x(-2), into new 2nd position
RHS = RHS[, 1] ~ mlag(Data[, 2], 2) ~ RHS[, 2:]
matrix ix_test = {1, 2} # index picking 1st and 2nd lag of causal var
if lagord > 2
# add further restricted lags of causal
matrix Xstar = Data[, 2] - 2 * cos(fpoint) * mlag(Data[, 2], 1) + mlag(Data[, 2], 2)
# (the Gegenbauer polynomial for non-causality at frequency freqpoint)
RHS = RHS ~ mlag(Xstar, seq(1, lagord - 2))
# (prev. line had a bug before, seq(3,lagord)
# (the Xstarlags need not come last, so they come before the const here)
endif
endif
matrix Depvar = Data[lagord+1: T, 1]
# here the original Gauss code was b=depvar/x
# AFAIK this means regress depvar on x and get the estimates
matrix temp, Varb
matrix b = mols(Depvar, RHS[lagord+1: T, ], &temp, &Varb) # temp: residuals
# scalar sig2 = (temp'temp) / (T - lagord - cols(RHS))
# matrix Varb = sig2 * invpd(RHS'RHS)
return qform(transp(b[ix_test]), invpd(Varb[ix_test,ix_test]))
```