```
# r: see function getS_C_J
# p: ditto
# output: the phase shifts (0;2pi) in 2nd col (1st the fpoints)
#
# If the rCov-matrix was passed, get the asy standard error of the
# respective phase in the 3rd col. (otherwise the col is zero);
# and (new2016) the three vectors:
# (1-c_a*)^2 + s_a*^2 = |F_a|^2
# s_b*^2 + c_b*^2 = |F_b*|^2
# c_b* (1-c_a*) - s_a* s_b* = D
# to check whether they are bounded away from zero as required.
if isnull(rCov)
matrix SaC = getS_C_J(r, fpoints, p)
else # need Jacobian to produce the variance of the phase
matrix Jac
SaC = getS_C_J(r, fpoints, p, &Jac)
if cols(Jac) != rows(rCov)
funcerr "mismatch between dims of Jac and the input rCov"
endif
endif
matrix Down = SaC[,3] .* (1-SaC[,5]) - SaC[,2] .* SaC[,4]
matrix Up = SaC[,2] .* (1-SaC[,5]) + SaC[,3].* SaC[,4]
matrix phase = atan2m(Up, Down)
matrix out = fpoints ~ phase
if !isnull(rCov) # add the delta-method standard errors
matrix phaseSE
loop i=1..rows(phase) --quiet
phaseSE |= sqrt(Jac[i,] * rCov * Jac[i,]') # just one col
endloop
else
phaseSE = zeros(rows(phase), 1) # dummy content
endif
out ~= phaseSE
# numerical checks
# part 1: (1-c_a*)^2 + s_a*^2 = |F_a|^2
matrix check1 = (1 - SaC[,5]).^2 + SaC[,4].^2
# part 2: |F_b*|^2
matrix check2 = SaC[,2].^2 + SaC[,3].^2
# part 3: D was calculated above
out ~= (check1 ~ check2 ~ Down)
return out
```

```
/*
The return matrix contains:
1st col: frequency points
2nd col: 'raw' phase shifts (0;2pi)
3rd col: standard error of the phase (asympt.)
4th col: the 'unwrapped' phase
5th col: the 'unwrapped' time delay
6th col: lower CI bound of that time delay
7th col: upper CI bound of that time delay
8th col: |F_a|^2 = (1-c_a*)^2 + s_a*^2
9th col: |F_b*|^2 = s_b*^2 + c_b*^2
10th col: D = c_rho = c_b* (1-c_a*) - s_a* s_b*
*/
wc = inlist(cond, const) ? 1 : wc # this is for backward-compatibility
# (explicit const in list "cond" overrules wc setting)
# (but position of const in cond is irrelevant, as gretl orders
# the constant always first)
factor = critical(N, (1-cover)/2)
# construct the regressors of the ARDL equ., and estimate
list ylags = (py == 0) ? null : lags(py,y)
list reg = lags(px,x) ylags
# further vars
reg = reg || cond
# constant term
reg = wc ? (const || reg) : reg
###
print "Using the following ARDL specification:"
string yname = argname(y)
string xname = argname(x)
printf " Dependent variable: %s\n", yname
printf " Driving variable: %s\n", xname
printf " Lagged endogenous terms: %d\n", py
printf " Included lags of driving variable: from 1 to %d\n", px
if !wc
print " Estimating w/o constant as requested."
endif
print "The confidence interval bounds in the returned matrix"
print " (6th & 7th cols 'lower'/'upper') are based on the"
printf " nominal coverage %4.3f.\n", cover
print "Underlying estimates:"
ols y reg
# the interesting coefficients
matrix ix = seq(1 + wc, wc + px + py) # without constant coeff.
matrix r = $coeff[ix]
matrix rCov = $vcv[ix, ix]
# get the phase calculation
matrix phwSE = getPhaseARDL(r, $pi*seq(1,99)'/100, px, &rCov)
matrix f = phwSE[,1]
# unwrap the phase
matrix phuw = phaseunwrap(phwSE[,1:2])
# divide by frequency to get the time shift
matrix tdluw = phuw[,2] ./ f # single col. here
# for the time delay also need to divide SE by the freq.
matrix tdlSE = phwSE[,3] ./ f # 3rd col is the SE
# the numerical checks
matrix checks = f ~ phwSE[,4:5]
# plot the delay
matrix tdlUWwCI = f ~ tdluw ~ (tdluw - factor*tdlSE) ~ (tdluw + factor*tdlSE)
if plotting
print "Producing plots..."
setopt gnuplot persist --with-lines --output=display
colnames(checks, "frequencies |F_a|^2 |F_b*|^2")
gnuplot 2 3 1 --matrix=checks
#
strings plotstrings = defarray("frequencies", sprintf("delay of %s, given %s", yname, xname), sprintf("lower (%4.3f)", (1-cover)/2 ), sprintf("upper (%4.3f)", 1 - (1-cover)/2) )
colnames(tdlUWwCI, plotstrings)
gnuplot 2 3 4 1 --matrix=tdlUWwCI --single-yaxis { set ylabel 'periods'; }
endif
# return everything, but w/o duplicating the freq column
matrix out = phwSE[,1:3] ~ phuw[,2] ~ tdlUWwCI[,2:4] ~ phwSE[,4:6]
colnames(out, "freq phase SE unwr_phase delay lower upper |F_a|^2 |F_b*|^2 D")
return out
```

```
/* An atan2 function for many inputs in one go (m= matrix).
U: Numerator (vector z x 1, "Up")
D: Denominator (z x 1, "Down")
Does "unwrapping" of arctan according to the signs of N and D:
1) D>0, U>0: standard, angle in (0, \pi/2) according to arctan
2) D>0, U<0: angle in (3*\pi/2, 2*\pi); arctan(U/D) + 2*\pi
3) D>0, U=0: depends on 'zerono' setting, 0 or 2*\pi
4) D=0, U>0: \pi/2
5) D=0, U<0: 3* \pi/2
6) D<0 : angle in (\pi/2, 3*\pi/2), namely arctan(U/D)+\pi
'zero2pi': Maps to angles in (0; 2*\pi) by default;
if zero2pi==0, then maps to (-\pi; \pi) instead.
(not yet implemented...)
'zerono': By default interprets a zero angle as 2*\pi, not zero.
*/
if !zero2pi
funcerr "Not yet implemented (zero2pi=0)"
endif
U = vec(U)
D = vec(D)
if rows(U) != rows(D)
funcerr "U and D must match"
endif
matrix angles = {}
loop i=1..rows(U) --quiet
if D[i] == 0 # can this actually hold exactly in reality?
# cases 4 and 5
angles = U[i] > 0 ? (angles | $pi/2) : (angles | 3*$pi/2)
else
UoD = U[i] / D[i]
if D[i] < 0 # case 6
angles |= atan(UoD) + $pi
elif U[i] > 0 # D > 0, case 1
angles |= atan(UoD)
elif U[i] < 0 # case 2
angles |= atan(UoD) + 2*$pi
elif U[i] == 0 # case 3
angles = zerono ? (angles | 2*$pi) : (angles | 0)
else
funcerr "shouldn't happen"
endif
endif
endloop
return angles
```

```
# r: input vector of stacked coefficients a and b
# p: the first p coeffs are from b(L), rest of a*(L)
# rCov: the V-Cov-matrix of r; if it is provided,
# then
#
# ARDL: (1-a*(L))y_t = b(L)x_{t-1}
# <=> y_t = a*(L)y_t + b*(L)x_t
#
# Output will have 5 columns:
# 1) f (frequencies)
# 2) s_\beta* = \sum^p b_j sin(j*f)
# -> sa (because of notation change in 2016)
# 3) c_\beta* = \sum^p b_j cos(j*f)
# -> ca
# 4) s_\alpha* = \sum^n a_j sin(j*f)
# -> sb
# 5) c_\alpha* = \sum^n a_j cos(j*f)
# -> cb
if (cols(r) > 1) || (cols(fpoints) > 1)
funcerr "need column vector inputs"
endif
matrix a = r[1:p,]
matrix b = r[p+1:,]
n = rows(b)
matrix argp = fpoints ** seq(1,p)
matrix vspT = sin(argp) ## each vsp is a row here!
matrix vcpT = cos(argp)
matrix sa = vspT * a # transp. w.r.t paper because all freqs
matrix ca = vcpT * a
matrix argn = fpoints ** seq(1,n)
matrix vsnT = sin(argn)
matrix vcnT = cos(argn)
matrix sb = vsnT * b
matrix cb = vcnT * b
if !isnull(Jacob) # also calculate the Jacobian
# each row will be J'; so what is vertically stacked
# in the paper is horizontally stacked here
matrix Denom1 = sa.^2 + ca.^2
matrix Denom2 = sb.^2 + (1-cb).^2
matrix Left = (ca ** ones(1,p)) .* vspT
Left -= (sa ** ones(1,p)) .* vcpT
Left = Left ./ (Denom1 ** ones(1,p))
matrix Right = ((1-cb) ** ones(1,n)) .* vsnT
Right += (sb ** ones(1,n)) .* vcnT
Right = Right ./ (Denom2 ** ones(1,n))
Jacob = Left ~ Right
endif
return (fpoints ~ sa ~ ca ~ sb ~ cb)
```

```
## Step 4: unwrap phases,
# i.e. adjust the previous phase when there are sudden jumps
#
# 'fphases' must have phases in 2nd col, fpoints in 1st col.
# Output:
# 2nd col the unwrapped phases,
# also returns in the 3rd col the step function for further use
# find the jump locations
matrix locdown = (diff(fphases[,2]) .< -thresh)
# (attn: spurious first entry in the above!)
locdown[1] = 0
matrix locup = (diff(fphases[,2]) .> thresh)
# construct the associated step function
matrix jstep = cum(locdown) - cum(locup)
# add all the necessary jumps
matrix phiuw = fphases[,2] + 2*$pi * jstep
return fphases[,1] ~ phiuw ~ jstep
```