```
/* 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.
mAsis: Contains variables to be included "as is", ie. without lagging, such as deterministics (except constant, which is always added later).
*/
matrix mcontainer = zeros(0,2) # 1st col freqs, 2nd col stats
T = rows(mY)
step = (tofreq - fromfreq) / (numoffreqs-1)
# handle the possibility of as-is-variables
matrix mtrans = {}
if exists(mAsis)
mtrans = mAsis
endif
loop for (f=fromfreq; f<=tofreq; f+=step) --quiet
teststat = caus_freq(mY, lagorder, f, mtrans)
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)
# 7b) b.determ (int)
#
# 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)
# 11) b.pdfpath (string, optional; can also be special value "display")
# 12) b.pltsource (bool, produce gnuplot source always; irrelevant if doplot==0)
#
# On exit the constructed matrix is put into the bundle as "b.BCoutput"
# The return value is 0 on success.
# check the input and sample
if !inbundle(b, mtarg) || !inbundle(b, mcause)
funcerr "Need target mtarg and potential cause mcause"
elif rows(b.mtarg) != rows(b.mcause)
print "Sample lengths:"
printf " Target: %d\n", rows(b.mtarg)
printf " Cause: %d\n", rows(b.mcause)
print "Please adjust (or try the higher-level function"
print " BreitungCandelonTest). Aborting."
funcerr "Unmatched sample lengths"
endif
# 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, "determ")
b.determ = 1 # code for constant only
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 !inbundle(b, "pltsource")
b.pltsource = 0
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
# put together
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 "%g ... %g.\n", b.fromfreq, b.tofreq
print "The relevant test statistic will be the minimal one."
printf "Calculating all pointwise results in 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 whole range (0, pi)."
print "Tests only valid pointwise for a pre-specified frequency."
print "Calculating all pointwise results for convenience."
print "-----------------------------------------------------------------"
endif
endif
# now do the actual test
# with deterministics
if b.determ == 1
print "Deterministics: constant only." # (added later)
matrix BCoutput = BCtest_fromto(mY, b.lagorder, b.fromfreq, b.tofreq, b.freqpoints)
elif b.determ == 2
print "Deterministics: constant and (centered) seasonals."
matrix BCoutput = BCtest_fromto(mY, b.lagorder, b.fromfreq, b.tofreq, b.freqpoints, { seasonals(1,1) } )
else
funcerr "other determs not implemented yet"
endif
# add the critical value as constant series for convenience
numoffreqs = rows(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 %g\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 %g\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
if inbundle(b, "pdfpath")
plotBC(BCoutput, b.pdfpath)
if b.pltsource # special request
plotBC(BCoutput)
endif
else
plotBC(BCoutput) # default behavior within this function
endif
endif
if b.verbosity
if findmin # a range was specified
scalar row = iminc(BCoutput[, 2])
printf "The min test stat (adjusted if needed) is %g,\n", BCoutput[row, 2]
printf " associated with frequency %g.\n", BCoutput[row, 1]
printf "The critical value at significance level %g is %g,\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 stats, "
printf "(3) crit val (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.
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))
print "Warning: missing values found in conditioning, will try to remove obs."
smpl condvars --no-missing
endif
endif
if determ == 2 && $pd == 1
print "Warning: non-seasonal dataset, ignoring chosen seasonals"
determ = 1
endif
print "=========== Frequency-specific Granger causality test ==========="
printf "Sample (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
string b.pdfpath = "display" # new behavior in 2.5
b.pltsource = 1 # produce gnuplot source anyway, as before
loop foreach s lagorder freqpoints siglevel fromfreq tofreq determ -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
if exists(pdfORdisp)
if pdfORdisp != "display"
string tail = pdfORdisp + (strlen(pdfORdisp) - 4)
if tail != ".pdf" && tail != ".PDF"
print "Warning, appending '.pdf' to plot path input"
pdfORdisp ~= ".pdf"
endif
endif
endif
if !exists(pdfORdisp) # default: output gnuplot source file
gnuplot 2 3 1 --matrix=BCoutput --with-lines=2 --with-lp=1 --single-yaxis
else # produce pdf or send to display
gnuplot 2 3 1 --matrix=BCoutput --with-lines=2 --with-lp=1 --single-yaxis --output=@pdfORdisp
endif
```

```
# this function calculates the actual test statistic
# Data: first col target, second col cause, further cols conditioning vars
# Asis: may be null or also {}, insert variables as-is
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 ~= mlag(Data[, ix_noX], seq(1, lagord))
# add the constant in all cases
RHS ~= ones(T, 1)
# add further as-is-Variables if present
if exists(Asis)
if cols(Asis) # could also be empty
RHS ~= Asis
endif
endif
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 ~= 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 ~= 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 polynomial for non-causality at frequency freqpoint)
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: resid, just placeholder)
return qform(transp(b[ix_test]), invpd(Varb[ix_test,ix_test]))
```