###### Abstract

Inference and hypothesis testing are typically constructed on the basis that a
specific model holds for the data. To determine the veracity of conclusions
drawn from such data analyses, one must be able to identify the presence of
the assumed structure within the data. In this paper, a model verification
test is developed for the presence of a random walk-like structure in the
variations in the frequency of complex-valued sinusoidal signals measured in
additive Gaussian noise. This test evaluates the joint inference of the random
walk hypothesis tests found in economics literature that seek random walk
behaviours in time series data, with an additional test to account for how the
random walk behaves in frequency space.

###### keywords:

model verification; structure identification; random walk frequency;
sinusoidal signal

1 1 0 13 May 2022 8 August 2022 https://doi.org/ Testing for a Random Walk
Structure in the Frequency Evolution of a Tone in Noise Testing for a Random
Walk Structure in the Frequency Evolution of a Tone in Noise Scarlett
Abramson1,*, William Moran1, Robin Evans 1 and Andrew Melatos 2 Scarlett
Abramson, William Moran, Robin Evans, and Andrew Melatos Abramson, S.; Moran,
W.; Evans, R.; Melatos, A. Correspondence<EMAIL_ADDRESS>
##  1 Introduction

Many analyses done today are model-based and inference results drawn from
these works rely heavily on strong assumptions to hold. One may ask if the
models used are the correct choice; and if true, or if not, what impact that
would have on the veracity of the results. The questions surrounding verifying
decisions in model selection have been popular fields of inquiry. (See Knuth
et al. (2015); Ding et al. (2018) for thorough reviews of model selection
methods.)

This paper considers a specific structure identification (SID) problem: to
determine if a given time series arose from noisy samples of a sinusoid whose
frequency varies according to a random walk (RW). This problem has important
practical applications in several areas of signal processing as discussed
briefly in Section 1.1 and is part of the long-established theoretical
tradition of optimal detection and estimation for tones with randomly varying
frequency. Our current motivation relates to the, as yet unsuccessful, search
for continuous gravitational waves (CGWs). The structure of the frequency
variation of such signals has astrophysical significance.

Since the seminal works of Kolmogorov and Wiener in the 1940s, the majority of
sensor-type signal processing approaches are model-based. Sensor-type signal
processing techniques based on carefully defined models of sensors and signals
have two huge advantages—firstly, optimal processing methods can be developed
by exploiting model structure, and secondly, performance can often be
characterized analytically.

There are, however, certain disadvantages. If the model is not an accurate
description of the sensor, the signals, or both, then the above advantages can
be diminished.

To address this issue, over the past 50 years many techniques for identifying
parameterized models of given structure from measured data have been developed
including adaptive methods which continually adjust the model parameters in
real-time based on the most recent data. Robust techniques have been developed
to maintain near-optimal performance under fairly broad modelling errors
including classes of structural errors.

The next important step in model-based sensor and signal processing problems
is the validation of models. In some situations, model validation can be based
on ground truth, though this is rarely possible.

Recently, there has been growing interest and activity in quantifying the
evidence for a model which is being used in a model-based optimal-sensor
signal processing scenario. That is, given a model and an optimal processing
solution—there is interest in determining the evidence for the model being
used based on observed data and optimally-processed outputs. A comprehensive
review covering methods for evaluating evidence for a model can be found in
the recent paper Llorente et al. (2020).

The theoretical foundations of general SID (as opposed to parameter
estimation) are still not well developed. Existing approaches such as Bayes
factor-based methods, model residual testing, and multiple hypotheses testing
all propose a structure, and then test in different ways how strongly the data
are consistent with the model. While this is an advance on parametric
estimation for a fixed given model, it is still a long way from true
structural discovery.

For the problem considered in this paper, an initial question is how to even
determine whether an RW is present in the data. This is a well-researched
field in the econometrics literature known as random walk hypothesis (RWH)
testing Hui (2012). Such methods try to determine from a recording of time-
series data whether the data behaves as an RW.

The second question is then how to distinguish between an RW and the
trigonometric function of an RW. One key quality that comes into play here is
that the resulting process is bounded so it is trivial to do this given
sufficiently a long recording time. However, many processes are random and
bounded and the RWH techniques would not distinguish the raw trigonometrically
transformed data by design. To distinguish, one must pick out the component
believed to act as an RW (the frequencies of the signal). Conventionally this
is done by the Fourier transform.

Once in Fourier space, one may then naively think to simply apply RWH testing
techniques and yield a viable result. A challenge comes from the fact that
frequency space is cyclic. Given that, it then becomes a challenge to
distinguish between something that would resemble an RW under RWH testing and
the Fourier transform of white Gaussian noise (WGN).

We argue here that should the trigonometrically transformed-RW in question
vary slowly, it could be distinguished from the Fourier transform of WGN and
so be susceptible to RWH testing methodologies.

We will present a non-nested multiple hypothesis procedure to test if a
structure (RW frequency signal) is a good model. The choice of hypotheses is
important in order to avoid structures that are RW-like but not actually RWs.
We create four hypotheses which enable sharp discrimination between RW
structure and structures which generate data close to an RW in frequency
space.

For completeness, we briefly outline a number of approaches to SID. We then
explain the logic underlying the multiple hypothesis structure we analyse in
this paper.

Bellman & Åström Bellman and Åström (1970) set out a criterion for structural
identifiability, and a generalised approach to SID problems. However, a key
challenge is the arbitrary nature of setting criteria for accepting a given
model over others while minimising complexity. One of the earliest attempts to
deal with this complex problem is the Akaike Information Criterion (AIC)
proposed in Akaike (1973), based on a form of dimension-penalization. However,
as Rissanen argued Rissanen (1982), Akaike’s criterion only minimizes total
predictive error in the degenerate case where one model size yields an
estimate with probability one. Rissanen proposed the Minimum Descriptor Length
(MDL) criterion Rissanen (1982), which is invariant under linear coordinate
transformations and consistent with respect to dimension estimation.
Additionally, Willems formalized SID problems through a behavioural approach
Willems (2007) that asks ‘given inputs and outputs to a system, does any
structure exist at all which may model such behaviours?’ ‘Tearing’ [framing
the system as a network of interconnected subsystems], ‘zooming’ [modelling
the subsystems], and ‘linking’ [modelling the connections between subsystems]
offer systematic approaches to the task. More recently, kernel-based methods
implemented in the machine learning methodology have been a popular tool (see
Hofmann et al. (2008) for a thorough review). Few of them have addressed the
issue of generating hypothesis tests on the model structure when applied to an
observed phenomenon. Kernel methods test whether a structure of some sort
works, without identifying it explicitly Hofmann et al. (2008); this is not
what we seek to do. Rather, we ask whether a specific structure accurately
models an observed phenomenon in a statistically significant sense.

We explore an SID problem, which we phrase through the lens of model
detection. This is intrinsically different from a signal detection problem
with a known (or more appropriately assumed) signal model. We wish to decide
whether any sinusoidal signal with random time-based variations in its
instantaneous frequency can provide a fit, with confidence, to the given data.
This, rather, seeks the presence of a structure in a record, not the presence
of a signal with a specified structure. However, what we observe is the signal
itself, and not its underlying parameters. Furthermore, here we seek whether
the data fits into a class of models, and so do not impose any specific model
(target signal) within the said class. Thus, given that within the class the
model parameters may be arbitrary, we assume that they are unknown. In the
case of detection, with unknown parameters, we then must implement estimation
procedures to verify whether the data resembles a sinusoidal signal with time
variation in its tone.

The parameter estimation problem of optimally estimating the frequency of a
constant frequency tone based on noisy measurements of the tone is among the
well-studied problem in signal processing with many papers exploring its finer
details. The general problem itself was formulated in discrete time by Rife
and Boorstyn Rife and Boorstyn (1974), following the work of Slepian Slepian
(1954) in continuous time. In approaching our problem, we assume that the
instantaneous frequency varies linearly over short time intervals (snapshots,
blocks) and follows an RW structure from snapshot to snapshot. A potential
extension would be block-wise (in time) frequency modulation of arbitrary
polynomial order (with respect to time); this would, in effect, extend the
problem to any continuous frequency modulation.

The statistical theory surrounding RWH testing Charles and Darné (2009) has
been shown to be effective at deciding in certain problems whether the stock
market can be said not to be following an RW Lo and MacKinlay (1988); Chow and
Denning (1993); Dickey and Fuller (1981), with implications for forecasting of
market trends.

However, it is argued that such tests can lead to ambiguous conclusions Hoque
et al. (2007) (in the sense of what conclusions drawn mean with regards to the
question) or inconsistent results Deo and Richardson (2003) (in the sense of
robustness). This is thought to arise from the sensitivity of variance-based
inference under the small-sample regimes Daniel (2001). Furthermore, because
of phase-wrapping behaviours, as discussed in the proof of Remark 2.1, such
methods do not discriminate between a signal measured in noise with an RW in
frequency and WGN.

The problem that we consider involves the detection of a model of a
trigonometric transformation of an RW process, rather than the detection of an
RW process. This is inherently different to standard RWH testing, as
trigonometric transformations are nonlinear in nature. The plausibility of
testing the RWH under such a transformation has been considered by other
researchers in the economics literature Granger and Hallman (1991); Wang and
De Jong (2013).

An auto-regressive (AR) time series $\\{\xi_{n}\\}_{n}$ of order 1, known as
an AR(1) process, has a recursive representation of the form
$\xi_{n}=\rho\xi_{n-1}+\varsigma_{n},$ with independent and identically
distributed (IID) driving noise $\\{\varsigma_{n}\\}_{n}$. The Dickey-Fuller
(DF) unit root test Dickey and Fuller (1981), assuming the presence of an
AR(1) process $\\{\xi_{n}\\}_{n}$ (not corrupted by noise) in the data, tests
a null hypothesis of a stationary (in the statistical sense) first-difference
process ($\xi_{n}-\xi_{n-1}$) which would occur if $|\rho|=1$ against an
alternative hypothesis $|\rho|\neq 1$.

It had been shown Wang and De Jong (2013) that the DF test, if applied to the
trigonometric transformation of an AR(1) process—asymptotically may still be
used to test for stationarity of the first difference process—following
analytical Ermini and Granger (1993), and empirical Granger and Hallman (1991)
arguments on the asymptotics of nonlinear transformations on 1st-order
integrated time series I(1), and subsequently on AR(1) processes.

However, this is different to the model we test for: the signal model of a
scaled trigonometric transformation of an AR(1) process with $|\rho|=1$, given
noise-corrupted measurements. As we do not assume, inherently, that a
(transformed) AR(1) time series is present within the data, but rather are
testing the validity of such a model—given the data—the DF test regime does
not apply here, and so we follow a different scheme that will be discussed
below.

In this paper, we advance a scheme for discriminating between WGN and a
sinusoidal signal whose instantaneous frequency executes a discrete-time RW,
which, when combined with RWH methods provides a more general inference on the
existence of an RW structure in frequency. In Section 2.1 the problem is
formulated mathematically. We design an estimator-detector in Sections
2.2–2.3, and discuss its analytical performance in Sections 3.1–4.

In continuous time the RW is replaced by a Wiener process and in this context,
prior to sampling, we are faced with signals of infinite bandwidth. This is,
as is often the case with mathematical models, a useful idealization. Whereas
a typical function of a Wiener process has infinite bandwidth, ultimately the
signal has to be measured and this forces a bandwidth constraint on the
measured signal imposed by the sensor. This has to be taken into account in
the sampling and subsequent processing of the discretized signal. We work
within the discrete time domain here and the bandwidth issues are overcome by
the determination of the sampling rate of the continuous time signal.

###  1.1 Applications

The concept of a time-varying auto-regressive (TVAR) tone appears across
various fields. Below are some applications of TVAR tones in the literature,
including the motivating application for this paper in gravitational wave (GW)
astronomy.

####  1.1.1 Gravitational Wave Astronomy

Predicted theoretically by Einstein in 1915, GWs are disturbances in space and
time, generated by the acceleration of asymmetric bodies, which propagate at
the speed of light Misner et al. (1973) as perturbations of the metric tensor
$g_{\mu\nu}$, which describes distances in spacetime by the invariant interval
$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}$ (with respect to displacement four-vectors
$dx^{\mu}$).

GWs cause oscillations in the proper displacements of freely falling test
masses and can be detected by long-baseline laser interferometers, such as the
Laser Interferometer Gravitational-Wave Observatory (LIGO) Aasi et al. (2015).

A GW incident normally on the plane of an interferometer stretches and
squeezes distances along the arms Misner et al. (1973); along two
characteristic polarisations ‘plus’ (whose principal axes of action align with
those of the spatial dimensions) and ‘cross’ (whose principal axes of action
are rotated $\frac{\pi}{4}$ with respect to those of the spatial dimensions).

At the time of writing, dozens of GW signals have been detected, including the
first discovery of the chirped tones from a binary black hole merger in 2015
Abbott et al. (2016) and a binary neutron star merger in 2017 Abbott et al.
(2017), inaugurating a new era in astronomy.

Instruments such as LIGO search for a variety of waveforms, not just chirps
from mergers. The theory also predicts the existence of persistent,
sinusoidal, quasimonochromatic CGW signals, believed to originate (for
example) from isolated or binary neutron stars Riles (2013). Under the biaxial
rotor model (a rigid body with two equal principal moments of inertia), a
pulsar emits GWs continuously at $f_{\star}$ and 2$f_{\star}$ Jaranowski et
al. (1998).

For an isolated neutron star, $f_{\star}$ decays monotonically over timescales
$\gtrsim 10^{3}\,$years, but wanders stochastically on smaller timescales
Shapiro and Teukolsky (1983); suspected mechanisms include seismic activity in
the crust or far-from-equilibrium avalanche processes in the superfluid
interior.

Whereas, in pulsar binaries, the stochastic wander of spin frequency is driven
by fluctuations in the accretion torque of gas from the binary companion de
Kool and Anzer (1993); Romanova et al. (2004); the underlying process is
suspected to be a result of transient disk formation Baykal et al. (1991), or
instabilities in the disk-magnetosphere Romanova et al. (2004).

The modelling and detection of such CGWs by interferometric data (such as from
aLIGO) follows the work of Jaranowski Jaranowski et al. (1998) (and the
following parts). The stochastic wander in $f_{\star}$ is incorporated into
detection schemes based on hidden Markov models Melatos et al. (2020);
Suvorova et al. (2016).

####  1.1.2 Structural Analysis

Materials used in construction exist in meta-stable states, and their
properties evolve stochastically over long timescales because of a variety of
factors (such as cracking or the influence of humidity). To ensure structural
stability one must be aware of the material’s natural frequency and harmonics,
to damp against excitation such as from applied wind pressure (for example,
the case of the Tacoma Narrows Bridge, “Galloping Gertie”) or from synchronous
loading (as seen with the London Millennium Bridge). When considering concrete
elements, for example, such a natural frequency wanders randomly in time.
Efforts are made to track this randomly wandering natural frequency in papers
such as Wu et al. (2018).

####  1.1.3 Sound Processing

The tracking of TVAR harmonics in signals is an approach considered in
communications and sound processing problems such as in the tracking of a
time-varying chirp (such as the Doppler shifts in mobile communications). Such
applications have been discussed in papers such as Dubois and Davy (2007) and
others.

##  2 Materials and Methods

###  2.1 Problem Formulation

To begin, we denote the discretisation (for indices $n=0,\ldots,N-1$,
$k=0,\ldots,K-1$) of a general temporal function $[\cdot]$ by:

$\displaystyle[\cdot]_{n}^{(k)}$ $\displaystyle=[\cdot](t_{n}^{(k)}),$ (1)
$\displaystyle t_{n}^{(k)}$ $\displaystyle=[n+k(N-1)]T,$ (2)

so that $t_{N-1}^{(k)}=t_{0}^{(k+1)}.$

This discretisation provides a partition for a data record of data sampled
uniformly at rate $f_{s}=\frac{1}{T}$ into $K$ blocks
$\\{x_{n}^{(k)}\\}^{k=0,\ldots,K-1}_{n=0,\ldots,N-1}$ of $N$ observations,
that is,

$\\{x_{n}^{(k)}\\}^{k}_{n}=\\{\\{x_{n}^{(0)}\\}_{n},\ldots,\\{x_{n}^{(K-1)}\\}_{n}\\}=\\{\\{x_{0}^{(0)},\ldots,x_{N-1}^{(0)}\\},\ldots,\\{x_{0}^{(K-1)},\ldots,x_{N-1}^{(K-1)}\\}\\}$

(where in the temporal discretisation structure noted above, we overlap the
last observation $x_{N-1}^{(k^{\prime})}$ and first observation
$x_{0}^{(k^{\prime}+1)}$ of any two adjoining blocks
$k=k^{\prime},k^{\prime}+1$) to a data record of $K(N-1)+1$ observations, or
$KN$ points including repeated points.

For a block $\\{x_{n}^{(k)}\\}_{n}$ (fixed $k$), we consider $N$ to be the
sample size of that subset

$\\{x_{n^{\prime}=(k-1)(N-1)+0},x_{n^{\prime}=(k-1)(N-1)+1},\ldots,x_{n^{\prime}=(k-1)(N-1)+(N-1)}\\},$

of the full record

$\\{x_{n^{\prime}=0},x_{n^{\prime}=1},\ldots,x_{n^{\prime}=K(N-1)}\\}.$

The parameter $N$ is the _blockwise sample size_. Graphically, the
discretisation (with overlap at endpoints) is constructed as in Figure 1.

Now, for any block (fixed $k$), let $\\{w_{n}^{(k)}\\}_{n}$ be samples of
complex WGN ($\mathbb{V}[\\{w^{(k)}_{n}\\}_{n}]=\sigma^{2}_{w}<\infty$ known),
and let $\\{s_{n}^{(k)}\\}_{n}$: $s(t_{n}^{(k)})=s^{(k)}_{n}$ be a sampled
complex sinusoid defined by,

$\displaystyle s_{n}^{(k)}$ $\displaystyle=Ae^{i\phi_{n}^{(k)}},$ (3)

where, for fixed $k$, $\\{\phi_{n}^{(k)}\\}_{n}$:
$\phi(t_{n}^{(k)})=\phi^{(k)}_{n}$ is the instantaneous phase of the sampled
sinusoid (assuming a quadratic structure in time). That is, for fixed $k$,
$\\{f_{n}^{(k)}\\}_{n}$: $f(t_{n}^{(k)})=f^{(k)}_{n}$ is linearly ramping
between $f_{0}^{(k)}$ and $f_{0}^{(k+1)}$ with respect to pivoting frequencies
$\\{f_{0}^{(k)}\\}_{k}$:

$\displaystyle
f^{(k)}_{n}=f_{0}^{(k)}+\frac{n}{(N-1)}[f^{(k+1)}_{0}-f^{(k)}_{0}],$ (4)

which corresponds to a phase:

$\displaystyle\phi^{(k)}_{n}=$
$\displaystyle\varphi_{{}_{0}}+\frac{1}{2}\frac{n^{2}}{(N-1)}T[f^{(k+1)}_{0}-f^{(k)}_{0}]$
$\displaystyle+\frac{1}{2}(N-1)T[f_{0}^{(k)}-f_{0}^{(0)}],$ (5)

where $\varphi_{{}_{0}}:=\phi^{(0)}_{0}$.

Figure 1: Visualisation of construction of block and blockwise sample indices
$k,n$ with $k=0,\ldots,K-1,$ $n=0,\ldots,N-1$ for partitioning $K$ blocks
$\\{x^{(k)}_{n}\\}_{n}$ (with end-point overlaps
$x^{(k)}_{N-1}=x^{(k+1)}_{0}$, $k$ fixed) of $N$ samples from data record
$\\{x_{n^{\prime}}\\}_{n^{\prime}}$ of $K(N-1)+1$ observations with
$n^{\prime}=0,\ldots,K(N-1)$.

We will say that the pivot frequencies $\\{f_{0}^{(k)}\\}_{k}$ are _stable_ if
they vary ‘slowly’ in the sense that,

$\displaystyle|m_{0}^{(k+1)}-m_{0}^{(k)}|$ $\displaystyle\ll N,$ (6)

where, as defined in (8), $\\{m_{0}^{(k)}\\}_{k}$ are the centre frequencies
of the Fourier bins corresponding to $\\{f_{0}^{(k)}\\}_{k}$ with
$-\frac{f_{s}}{2}$$\leq$$f_{0}^{(k)}$$\leq\frac{f_{s}}{2}.$

We note that this is effectively a band-limiting statement on the signal.
Centre frequencies of Fourier bins of the signal are limited in their
variation between consecutive bins to be significantly less than the length
$N$ of a block. Effectively this prevents wrapping of the frequencies within a
block and aliasing.

By means of RWH testing, one cannot distinguish between an RW in frequency
space that violates condition (6) and the sequence of peak frequencies
$\\{f_{w}^{(k)}\\}_{k}$ of blocks (fixed $k$) of WGN $\\{w_{n}^{(k)}\\}_{n}$
when processing via Fourier transforms.

Existing RWH tests, such as the LOMAC Test Lo and MacKinlay (1988); or the
Chow-Denning Test Chow and Denning (1993), are structured around variance
ratio (VR) testing (for an overview of such tests see Charles and Darné
(2009)). If a process $\\{\varsigma_{n}\\}_{n}$ does not satisfy condition
(6), then for any realistic $N$, the existing RWH tests would view the process
$\\{[f(\varsigma)_{n}]_{f_{s}}\\}_{n}$ as it would view the peak frequencies
$\\{[f(\omega)_{n}]_{f_{s}}\\}_{n}$ of blocks of WGN $\\{\omega_{n}\\}_{n}$.

To clarify, a variance ratio test would not be able to determine whether
$\\{[f(\varsigma)_{n}]_{f_{s}}\\}_{n}$ was drawn from
$\mathcal{N}(0,[\sigma^{2}_{\omega}]_{f_{s}}),$ or from
$\mathcal{N}([\mu_{\xi}(t)]_{f_{s}},[\sigma^{2}_{\xi}(t)]_{f_{s}})$ (where
$\\{\xi(t)\\}$ is some RW on $F$).

Another issue that is faced, is the wrapping behaviour of frequency spaces at
the edges. If the frequency variations do not occur sufficiently far from the
boundaries—then there would also be issues distinguishing an RW frequency
process, that is $f_{s}$-modular, from the peak frequencies of a WGN process,
that are $f_{s}$-modular.

To clarify, should a process $\\{\xi_{n}\\}_{n}$ in the frequency space $F$
step outside the domain by some small amount $\epsilon>0$: $\epsilon\ll
f_{s}$, it would appear to jump to the other boundary. Because of wrapping
issues, an RWH testing algorithm, utilising Fourier transforms, would read
this as having jumped a length $f_{s}-\epsilon$ either by
$\frac{f_{s}}{2}\downarrow-\frac{f_{s}}{2}+\epsilon$ or
$-\frac{f_{s}}{2}\uparrow\frac{f_{s}}{2}-\epsilon$, rather than the
potentially infinitesimal jump of length $\epsilon$ of
$\frac{f_{s}}{2}\uparrow\frac{f_{s}}{2}+\epsilon$ or
$-\frac{f_{s}}{2}\downarrow-\frac{f_{s}}{2}-\epsilon$. Such a process would
then appear to have variance comparable to the size of the domain itself and
would be indistinguishable from the frequencies of a WGN process.

Recalling that on an unbounded interval, an RW process is also unbounded, and
thus the sinusoid with an RW frequency would have infinite bandwidth in the
asymptotic case. However, the mere act of measurement is a finite process in
finite time—which then imposes a finite bandwidth upon the process.

By standard practice, the band-limiting of the signal would be determined by
an anti-aliasing filter setting the lower bound on the sampling frequency to
the Nyquist limit. Then, one would be able to choose the block-wise sample
size determined by the sampling frequency as per condition (6), and the other
parameters can be determined to optimise the efficiency of the verification
process.

For a more detailed discussion of how to optimally set out the parameters,
refer to Section 4.

We discretise the frequency range $[-\frac{f_{s}}{2},\frac{f_{s}}{2}]$ into
Fourier bins (indexed by $m=0,\ldots,N-1$) of width $L=\frac{1}{NT}$, and
represent the effective frequency $f(m)$ of a bin $m$ by:

$\displaystyle f(m)=\begin{cases}\frac{m}{NT},\quad\quad
m=0,\ldots,\lfloor\frac{N}{2}\rfloor-1,\\\ \frac{m-N}{NT},\quad
m=\lfloor\frac{N}{2}\rfloor,\ldots,N-1.\end{cases}$ (7)

Then the discretisation in time (and subsequently frequency) is equivalent to,

$\displaystyle f_{n}^{(k)}$ $\displaystyle=f(m_{n}^{(k)})+\delta_{n}^{(k)},$
(8)

with $n=0,\ldots,N-1$, and where $\delta_{n}^{(k)}$ is a displacement term
representing uncertainty in phase due to resolution of frequency; if nothing
is known about its distribution we assume
$\delta_{n}^{(k)}\sim\mathcal{U}(-\frac{L}{2},\frac{L}{2})$. This is a non-
physical prior which will impose little to no bias on inferences, used
similarly in discretised methods such as in HMMs Suvorova et al. (2016);
Melatos et al. (2020).

We will say that a process $\\{\xi_{n}\\}_{n}$ is _RW-like_ if, when testing
it against the RWH, its behaviour would not be considered significantly
different from that of an RW, in the statistical sense.

In terms of AR(1) models, of the form $\xi_{n}=\rho\xi_{n-1}+\varsigma_{n},$
$\\{\varsigma_{n}\\}_{n}$ IID, we would say that a process $\\{\xi_{n}\\}_{n}$
is _RW-like_ if, when testing it against the RWH, one cannot say in a
statistically significant sense that $|\rho|\neq 1$.

The importance of appropriately primed SID approaches on data series assumed
to be containing RWs prior to significance tests based on that assumption is
made clear by Durlauf and Phillips Durlauf and Phillips (1988); who
analytically characterised the behaviour of regression coefficients for time
trends in data, assuming trend stationary data, on I(1) processes (a class
under which RW-like series fall).

They found, analytically, that regressions of an RW against a time trend will
yield incorrect inferences of a greater significance of a trend than the
present work, resulting in a stronger bias for hypotheses assuming time trends
than appropriate; supporting previous empirical works done via Monte Carlo
(MC) simulations Nelson and Kang (1981).

For time-series data tested against the RWH, alternative hypotheses of RW-like
series may be distinguished from an RW by spectral methods Durlauf (1991) by
utilising the weak convergence (in the Hilbert sense, that is $\langle
x_{n},y\rangle\to\langle x,y\rangle\forall y$) under the sup
metric$d(f,g)=\sup_{x}\\{d^{\prime}(f(x),g(x))\\}$ of the periodogram
deviation process (the sequence of deviations of the time-series in question
from the periodogram of a true RW process) to a Brownian bridge on $C[0,1]$.

This is not a convergence that is preserved under trigonometric transforms,
and thus we would not be able to apply equivalent tests in our case.

This, taken with Remark 2.1, yields that the following hypothesis testing
structure is enough to detect an RW in frequency space, as per our model:

$\displaystyle\begin{cases}\mathcal{H}_{0}:\ x_{n}^{(k)}=w_{n}^{(k)},\\\
\mathcal{H}_{1}:\ x_{n}^{(k)}=s_{n}^{(k)}+w_{n}^{(k)},\\\
\quad\quad\quad\quad\\{f_{0}^{(k)}\\}_{k}\text{ `stable', RW-like},\\\
\mathcal{H}_{2}:\ x_{n}^{(k)}=s_{n}^{(k)}+w_{n}^{(k)},\\\
\quad\quad\quad\quad\\{f_{0}^{(k)}\\}_{k}\text{ `stable', not RW-like},\\\
\mathcal{H}_{3}:\ x_{n}^{(k)}=s_{n}^{(k)}+w_{n}^{(k)},\\\
\quad\quad\quad\quad\\{f_{0}^{(k)}\\}_{k}\text{ unstable, not RW-
like}.\end{cases}$ (9)

In terms of AR(1) models, of the form $\xi_{n}=\rho\xi_{n-1}+\varsigma_{n},$
$\\{\varsigma_{n}\\}_{n}$ IID, instability can be considered to mean
$|\rho|>1$ and non-RW-likeness can be considered to mean $|\rho|\neq 1$ in a
statistically significant sense.

The reasoning behind these hypothesis tests will be expanded on, once the
necessary tools are extracted from the literature and developed for our
problem’s lens (see the discussion at the end of Section 2.3).

###  2.2 Carrier Frequency Estimation

Given that the parameters ($A$, $f$, $\phi$) of the signal model are assumed
to be unknown, frequency estimation is necessary to apply the hypothesis test
introduced in (9).

Following Niebauer (2013) the carrier frequency of a linear chirp is taken to
be the average value of the lower and upper frequencies of the peak region of
the Fourier-transformed signal.

The ‘slow’ nature of the RW in frequency (6) asserts that the endpoints of the
peak region are separated by a negligible distance with respect to the size of
the frequency space.

Under the assumption of the ‘slowness’ of the RW, the signal’s peak power is
distributed across the Fourier bins in a manner that would resemble a narrow
peak in Fourier space.

The coarse search process given by Rife and Boorstyn Rife and Boorstyn (1974)
estimates the constant-valued frequency of a single-toned sinusoid. Given the
sharpness of the peak in Fourier space, the methods in Rife and Boorstyn
(1974) are an appropriate means of estimating the carrier frequency of the
signal.

We note that the estimation process here is not necessarily optimal, but
rather is implemented for its well-understood statistical properties.

Taking the $\frac{1}{N}$-normalised Discrete Fourier Transform (DFT)

$\displaystyle X_{m}^{(k)}$
$\displaystyle=\mathcal{F}[\\{x_{n}^{(k)}\\}_{n}](m)=\frac{1}{N}\sum_{n=0}^{N-1}x^{(k)}_{n}\exp{[\frac{-2\pi
inm}{N}]},$ (10)

where $m=0,\ldots,N-1$ is the Fourier bin number, we estimate the location of
the carrier frequency

$\displaystyle f^{(k)}_{c}$
$\displaystyle=\frac{f^{(k)}_{0}+f^{(k+1)}_{0}}{2},$ (11)

by the method in Rife and Boorstyn (1974), as:

$\displaystyle\hat{m}_{c}^{(k)}$ $\displaystyle=\underset{0\leq m\leq
N-1}{\text{argmax}}|X_{m}^{(k)}|^{2}.$ (12)

We take $\hat{f}_{c}^{(k)}:=f(\hat{m}_{c}^{(k)})$, where carets indicate an
estimated quantity. We represent the carrier frequency estimator
$\hat{f}_{c}^{(k)}$, at fixed $k$, by the ‘true’ carrier frequency of the
signal displaced by an estimation error:

$\displaystyle\hat{f}_{c}^{(k)}$
$\displaystyle=f_{c}^{(k)}+\varepsilon^{(k)}.$ (13)

The variance of the peak location of a sinusoidal signal in noise, as
estimated according to (12) and (13) is shown to be $\sigma_{\varepsilon}^{2}$
= $\frac{1}{N}\frac{6}{(2\pi NT)^{2}\text{SNR}}$ by Stoica et al. (2005);
Stoica and Nehorai (1989) with SNR = $\frac{A^{2}}{\sigma_{w}^{2}}.$

Given limiting distributions of periodogram frequency estimators Nabeya and
Tanaka (1986) are normal to leading order, and the asymptotic normality of the
Maximum Likelihood Estimator (MLE), we will assume that $\varepsilon^{(k)}$
follows a distribution of the form $\varepsilon^{(k)}$ $\overset{a}{\sim}$
$\mathcal{N}$(0, $\sigma_{\varepsilon}^{2}$) for sufficiently large $K$. Note,
however, that the rate of convergence here is slow—with Theorem 8 of Nabeya
and Tanaka (1986) stating effectively that for the standardised version
$\overline{\varepsilon}$ of $2\pi\varepsilon$:

$\displaystyle\mathbb{P}(\overline{\varepsilon}<\xi)$
$\displaystyle=\Phi(\xi)\left[1+\frac{1}{K^{2/5}}\frac{H_{e_{3}}(\xi)}{3}\right.$
(14)
$\displaystyle\quad+\frac{1}{K^{4/5}}\left(\frac{H_{e_{4}}(\xi)}{4}+\frac{H_{e_{6}}(\xi)}{18}\right)$
$\displaystyle\quad+\frac{1}{K^{6/5}}\left(\frac{H_{e_{5}}(\xi)}{5}+\frac{H_{e_{7}}(\xi)}{12}+\frac{H_{e_{9}}(\xi)}{162}\right)$
$\displaystyle\quad+\frac{1}{K^{8/5}}\left(\frac{H_{e_{6}}(x)}{6}+\frac{47H_{e_{8}}(x)}{480}+\frac{H_{e_{10}}(x)}{72}\right.$
$\displaystyle\quad\left.\left.+\frac{H_{e_{12}}(x)}{1944}\right)+\mathcal{O}(K^{-2})\right],$

where
$H_{e_{r}}(x)=(-1)^{r}e^{\frac{x^{2}}{2}}\frac{\rm{d}^{r}}{\rm{d}x^{r}}e^{-\frac{x^{2}}{2}}$
is the probabilistic Hermite polynomial of order $r$.

As will be noted again later in further discussion, the simulations failed to
provide results for SNRs less than $-20$ dB. We argue that the procedure
provides a reasonable degree of performance for Signal-to-Noise Ratios greater
than $-20$ dB, however, in the lower SNR regime, we would argue, that this is
likely due to the estimation procedures used as a proof-of-concept for the
broader issue of model verification in the problem being considered. More
sophisticated estimation approaches could be applied, however, that is not the
focus of the paper as this is not an estimation-detection problem we do not
need to know whether the model verification lens in this problem can be used,
and the SNR-regime of SNRs greater than $-20$ dB is sufficient to that end.

We do, however, note that Figure 2a,b shows that the parameters ($T$, $K$) do
not affect the performance of the procedure as one would anticipate.

###  2.3 Detection of Random Walk Structure

Under $\mathcal{H}_{1},$ by (11) we know that $f_{c}^{(k)}$ is the average of
two RW elements $f_{0}^{(k)},$ $f_{0}^{(k+1)}$ ($k$ fixed) sampled from an RW
$\\{f_{0}^{(k)}\\}_{k}$ of length $K$. Recalling that the variance of an RW of
finite length is its length, then, under $\mathcal{H}_{1},$ we assert:

$\displaystyle\mathbb{V}[f_{c}^{(k)}|\mathcal{H}_{1}]$ $\displaystyle\propto
K,$ (15)

where we use $\mathbb{V}$[$\cdot$] to denote the ensemble variance operator.

Now, (15) tells us that the same simplifications used in the analysis of RWs
as in Charles and Darné (2009) should be used here for ‘regularity’ (in a
statistical sense). In line with this, we define the first difference process
$\\{y_{c}^{(k)}\\}_{k}$ (for $k=1,\ldots,K-1$),

$\displaystyle y_{c}^{(k)}$ $\displaystyle:=f_{c}^{(k)}-f_{c}^{(k-1)}.$ (16)

We also define, at lag $1\leq\tau\leq\frac{K}{2},$ a general difference
process $\\{\Delta_{\tau}^{(k)}\\}_{k}$ (for $k=\tau,\ldots,K-\tau$):

$\displaystyle\Delta_{\tau}^{(k)}$
$\displaystyle:=f_{c}^{(k)}-f_{c}^{(k-\tau)}.$ (17)

The power of WGN (i.e., $\\{|X_{m}^{(k)}|^{2}\\}_{m}$ under $\mathcal{H}_{0}$)
is equally distributed between Fourier bins. Through some algebraic
manipulation, we find that, under $\mathcal{H}_{0}$:

$\displaystyle\hat{f}_{c}^{(k)}=\left.\begin{cases}-\frac{1}{2T}\\\
-\frac{1}{2T}+\frac{1}{NT}\\\ \quad\quad\vdots\\\ -\frac{1}{NT}\\\ \quad\ 0\\\
\quad\frac{1}{NT}\\\ \quad\quad\vdots\\\ \quad\frac{1}{2T}-\frac{2}{NT}\\\
\quad\frac{1}{2T}-\frac{1}{NT}\end{cases}\right\\}\text{with probability
}\frac{1}{N}.$ (18)

Then, under $\mathcal{H}_{0},$ by standard definitions and properties of the
distributions (see Appendix A), we find that
$\mathbb{V}\left[\hat{f}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]=\frac{1}{12}\left[\left(\frac{1}{T}\right)^{2}-\left(\frac{1}{NT}\right)^{2}\right].$
And so, by (16) define,

$\displaystyle\sigma^{2}_{0}$
$\displaystyle:=\mathbb{V}\left[\hat{y}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]=\frac{1}{6}\left[\left(\frac{1}{T}\right)^{2}-\left(\frac{1}{NT}\right)^{2}\right].$
(19)

Now, combining (8), (11), (13) and (16), we may represent
$\hat{y}^{(k)}_{c}=\hat{f}_{c}^{(k)}-\hat{f}_{c}^{(k-1)}$ ($k$ fixed) under
$\mathcal{H}_{1},$ $\mathcal{H}_{2},$ or $\mathcal{H}_{3},$ by:

$\displaystyle\hat{y}^{(k)}_{c}=$
$\displaystyle\frac{1}{2}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\right]$
$\displaystyle+\frac{1}{2}\left[\delta^{(k+1)}_{0}-\delta^{(k-1)}_{0}\right]$
$\displaystyle\ \ +\left[\varepsilon^{(k)}-\varepsilon^{(k-1)}\right],$ (20)

We assume that, under $\mathcal{H}_{1},$ we may use a representation of the
form:

$\displaystyle m_{0}^{(k)}$ $\displaystyle=m_{0}^{(k-1)}+u^{(k)},$ (21)

where the discretised driving term of the RW-like structure is defined to be a
sequence of integer-valued jumps between Fourier bins, represented by
$\\{u^{(k)}\\}_{k}$.

To illustrate the range of behaviours, we consider, at extreme ends, two
different forms of $\\{u^{(k)}\\}_{k}$. The first case is:

$\displaystyle u^{(k)}$ $\displaystyle=\left.\begin{cases}+1\\\ +0\\\
-1\end{cases}\right\\}\text{ with probability }\approx\frac{1}{3}.$ (22)

This discrete uniform jump structure represents the case where minimal
information is assumed for the jumps satisfying (6).

The second case considered is:

$\displaystyle u^{(k)}$ $\displaystyle=\text{nint}\left[z^{(k)}\right],$ (23)

where nint[$\cdot$] is the nearest integer to [$\cdot$], and $z^{(k)}$ $\sim$
$\mathcal{N}$(0,1).

Then,

$\displaystyle\sigma_{1}^{2}$
$\displaystyle:=\mathbb{V}\left[\hat{y}_{c}^{(k)}\middle|\mathcal{H}_{1}\right]$
(24) $\displaystyle=\frac{1}{N}\frac{12}{(2\pi
NT)^{2}\text{SNR}}+\begin{cases}\frac{9}{24}(\frac{1}{NT})^{2},\quad\text{ for
uniform jumps},\\\ \frac{14}{24}(\frac{1}{NT})^{2},\quad\text{ for normal
jumps}.\end{cases}$

Thus, for this problem, the ‘slowly’ varying nature, in the sense of condition
(6), of the RW in frequency space under $\mathcal{H}_{1},$ is manifest by
$\sigma_{0}^{2}>>\sigma_{1}^{2}.$

Given the presence of the discretisation displacement terms in representation
(20), the data is non-gaussian regardless of the jump structure considered.
Motivated by the asymptotic properties of variance estimators for non-Gaussian
data, as shown in O’Neill (2014), we construct the test statistic,

$\displaystyle\hat{\chi}_{0}^{2}:=\hat{\text{DoF}}\frac{\hat{\sigma}^{2}}{\sigma_{0}^{2}}\overset{a}{\sim}\chi^{2}(\hat{\text{DoF}}),$
(25)

where $\hat{\text{DoF}}=\frac{2(K-1)}{\hat{\kappa}-\frac{K-4}{K-2}}$ is the
sample Degrees of Freedom estimator, with $\hat{\kappa},$ $\hat{\sigma}^{2}$
the sample kurtosis and sample variance estimators respectively, for
$\\{y_{c}^{(k)}\\}_{k}$. As explained in Douillet (2009), we require $K\gtrsim
11$ for the asymptotic properties in (25) to hold to at least 1 sigma.

To normalise the test statistic, we perform a Wilson-Hilferty transform,

$\displaystyle\hat{Z}_{0}:=\frac{\sqrt[3]{\frac{\hat{\sigma}^{2}}{\sigma_{0}^{2}}}-(1-\frac{2}{9\hat{\text{DoF}}})}{\sqrt{\frac{2}{9\hat{\text{DoF}}}}}\overset{a}{\sim}\mathcal{N}(0,1).$
(26)

We reject $\mathcal{H}_{0}$ to Probability of False Alarm (PFA)
$\alpha\in$(0,1) if we have that $\hat{Z}_{0}<\Phi^{-1}$($\alpha$), where
$\Phi$($\cdot$) denotes the Cumulative Distribution Function (CDF) of the
Standard Normal distribution.

However, given that the above procedure seeks _controlled_ (in the sense of a
‘slow’ wander as defined in (6)), but RW-like, jumps in the frequency
estimators $\\{\hat{f}_{c}^{(k)}\\}_{k}$; any stable (in the sense of a mean-
reversion type behaviour; $|\rho|\leq 1$) random process bounded on a
sufficiently small domain would also be flagged by this test as statistically
significant. Thus, not only does $\hat{Z}_{0}$ distinguish between
$\mathcal{H}_{0}$ and $\mathcal{H}_{1}$, but also $\mathcal{H}_{2}$ and
$\mathcal{H}_{3}$.

The hypothesis problem (9) addressed here prompts three questions to be
tested:

  1. 1.

Does the measured frequency exhibit a structure resembling an RW?

  2. 2.

Assuming that the measured frequency does exhibit an RW-like structure, does
that structure wander in a ‘slowly’ (in the sense of definition (6)), rather
than randomly jumping in a manner that could be characterised by randomly
selecting Fourier bins in a uniformly-distributed manner? That is to say,
assuming that the answer to Question 1 was affirmative, then, is that a true
RW-like structure, or an artefact of the estimation procedure acting on pure
noise?

  3. 3.

Does the measured frequency revert to the mean, as in an Ornstein-Uhlenbeck
process? That is to say, does $\\{f_{0}\\}^{(k)}$ wind down? (For good
resources on processes like Ornstein-Uhlenbeck processes, we refer the reader
to Klebaner (2005)). Mean reversion is likely in astrophysical and GW
applications, e.g., signals from rotating neutron stars Mukherjee et al.
(2018).

In the economics literature, VR-based RWH testing of time series data (such as
the LOMAC Test Lo and MacKinlay (1988); or the Chow-Denning Test Chow and
Denning (1993)), determines whether the ratio of the variance of the first
difference process (16) to the variance of the difference process (17) at lag
$2\leq\tau\leq\lfloor\frac{K}{2}\rfloor$ ($\tau$ fixed) departs from unity in
a statistically significant sense. This determines whether mean-reversion
behaviours are present, as seen in an Ornstein-Uhlenbeck process.

Whereas alternative RWH testing approaches seek the presence of a unit root
against stationarity (such as the DF test Dickey and Fuller (1981), or
Bhargava’s von Neumann ratio test Bhargava (1986) for the finite sample
regime), there are stationary processes that appear like RWs from the lens of
unit root testing Nicolau (2002).

It has been shown in Diebold and Rudebusch (1991) (and other papers) that the
DF unit root test statistic $\hat{\rho}$ has low power against edge cases
(near the unit root) and against trend-stationary processes; resulting in
incorrect conclusions being drawn on testing hypotheses for discriminating
between an RW and mean-reverting processes (such as Ornstein-Uhlenbeck
processes) via unit root tests.

Furthermore, many traditional detectors (such as F test statistics, Hausman
test statistics, or Durbin-Watson test statistics) which may be framed to test
for (or test against) stationarity tend to yield inconclusive results when
tested in an RWH framework Durlauf and Phillips (1988).

Lastly, it has been shown that the LOMAC and Chow-Denning test statistics are
considerably more powerful than the Dickey-Fuller test statistic and other
unit root tests Lo and MacKinlay (1988); Chow and Denning (1993) against even
borderline cases such as distinguishing between AR(p) and ARIMA(p,d,q) models;
and the Chow-Denning test statistic more powerful than the LOMAC Chow and
Denning (1993).

Thus, for the purposes of this study, we consider VR-based approaches (such as
LOMAC or Chow-Denning) and not unit root-based approaches (such as DF).

The first question is one inherent in the RWH testing literature, and so, such
methodologies try to address this question.

The second question can be covered by the Controlled Variation Test defined
above (26)—which checks that the first difference process in frequency is
varying in a sufficiently controlled manner—that it is more than likely not
just randomly picking tones in a uniform manner across the bins.

The third question, similarly to the first question, can be covered by an RWH
test in determining if the variance of the lagged differences changes over
time in a statistically significant sense.

If the first difference process (16) is serially uncorrelated, the frequencies
do not form an Ornstein-Uhlenbeck process. Then, one could select a sequence
$\\{\tau_{1},\ldots,\tau_{J}\\}$ of $J$ of lags
($2\leq{}\tau_{j}{}\leq{}\lfloor\frac{K}{2}\rfloor$), and construct a Chow-
Denning Chow and Denning (1993) test statistic:

$\displaystyle V_{1}=\max_{1\leq j\leq J}|M_{1}(\tau_{j})|,$ (27)

where $M_{1}(\tau)$ is the LOMAC Lo and MacKinlay (1988) test statistic (2
$\leq$ $\tau$ $\leq$ $\lfloor\frac{K}{2}\rfloor$):

$\displaystyle M_{1}(\tau)=\frac{VR(\tau)-1}{\varphi(\tau)}$ (28)

with asymptotic variance $\varphi(\tau)$ defined, as in Charles and Darné
(2009); Lo and MacKinlay (1988) by:

$\displaystyle\varphi(\tau)=\frac{2(2\tau-1)(\tau-1)}{3\tau},$ (29)

and $VR(\tau)$ is the variance ratio defined, as in _loc. cit._ by:

$\displaystyle
VR(\tau)=\frac{\sum_{k=\tau}^{K-\tau}(y_{c}^{(k)}+\cdots+y_{c}^{(k-\tau+1)}-\tau\hat{\mu})^{2}/\tau}{\sum_{k=1}^{K-1}(y_{c}^{(k)}-\hat{\mu})^{2}},$
(30)

where $\hat{\mu}$ is the sample mean estimate for $\\{y_{c}^{(k)}\\}_{k}$.

One could infer between the two tests, taking similar assumptions of
independence as taken by Chow and Denning Chow and Denning (1993), with
appropriate Šidák corrections. Then, to PFA $\alpha\in(0,1),$ we argue that:

  * •

The tones vary in a sufficiently controlled manner by test (26) if
$\hat{Z}_{0}<\Phi^{-1}(\alpha^{*})$; and,

  * •

We reject the RWH by test (27) if $V_{1}>\Phi^{-1}(1-\frac{\alpha^{*}}{2})$,

where $\alpha^{*}=1-(1-\alpha)^{1/J}$ is the Šidák correction for the joint
inference of a Chow-Denning Test (27) of size $J$. $J$ is taken as the Šidák
correction factor, and not $J+1$ (as would be implied by the additional
independent test from the Controlled Variations test), due to the non-nested
nature of the tests Ulloa and Pesaran (2006).

Set hypothesis tests for the Controlled Variations Test (26) and the (size
$J$) Chow-Denning Test (27), respectively, as:

$\displaystyle\begin{cases}\mathcal{H}_{0}^{1}:\ \\{f_{0}^{(k)}\\}_{k}\
\text{`unstable'},\\\ \mathcal{H}_{1}^{1}:\ \\{f_{0}^{(k)}\\}_{k}\
\text{`stable'},\end{cases}$ (31)
$\displaystyle\begin{cases}\mathcal{H}_{0}^{2}:\ \\{u^{(k)}\\}_{k}\ \text{is
serially uncorrelated},\\\ \mathcal{H}_{1}^{2}:\ \exists\text{ lag }\tau\text{
so that }\\{f_{c}^{(k)}-f_{c}^{(k-\tau)}\\}_{k}\\\ \quad\quad\quad\quad\text{
is serially correlated at lag }\tau.\end{cases}$ (32)

The joint inference of the above hypothesis tests determines which of the four
hypotheses holds as described in (9).

Given that the two tests are independent of each other (as per the assumptions
that the Chow-Denning test statistic is constructed upon Chow and Denning
(1993)), our testing structure falls under the non-nested hypothesis testing
methodologies Ulloa and Pesaran (2006).

The four possible results could arise from this joint inference problem:

  * •

If the test (27) asserts that one cannot reject the RWH and the test (26)
asserts that the tones vary in a sufficiently controlled sense: then one could
argue that there is an RW in the frequency.

  * •

If the test (27) asserts that the data rejects the RWH and the test (26)
asserts that the tones vary in a sufficiently controlled sense: then one could
argue that the frequency follows a stable process that does not resemble an RW
(an Ornstein-Uhlenbeck process, for example).

  * •

If the test (27) asserts that one cannot reject the RWH and the test (26)
asserts that the tones do not vary in a controlled sense: then one could argue
there is just noise.

  * •

If the test (27) asserts that the data rejects the RWH and the test (26)
asserts that the tones do not vary in a controlled sense: then one could argue
that the frequency follows an unstable process that does not resemble an RW.

To elaborate, should the Chow-Denning tests (27) argue that it is not
statistically significant to be able to reject the RWH, then by Definition 2.1
the central instantaneous frequency process $\\{f_{c}^{(k)}\\}^{k}$ as
estimated from the observation process $\\{x_{n}^{(k)}\\}_{n}^{k}$ is RW-like,
that is that it has some component behaving in a manner that is
indistinguishable from that of an RW under RWH testing methods.

As RW-like is defined (Definition 2.1 in a binary sense, as is ‘stability’
(Definition 2.1), we may take a Fischer-esque philosophical interpretation of
the results (i.e., if it is not statistically significant to say that it is A,
then it is statistically significant to say it is not A).

The joint inference by the two above tests, are summarised with respect to the
hypotheses (9), (31) and (32) in Table 1.

Table 1: Joint inference table. |  $\mathcal{H}^{2}_{0}$ Holds |  Reject $\mathcal{H}^{2}_{0}$
---|---|---
$\mathcal{H}^{1}_{0}$ Holds |  $\mathcal{H}_{0}$ holds |  $\mathcal{H}_{3}$ holds
Reject $\mathcal{H}^{1}_{0}$ |  $\mathcal{H}_{1}$ holds |  $\mathcal{H}_{2}$ holds

Each case explicitly covers one of the above four possible cases on the
constructed problem. Thus, the hypothesis tests as described in (9) are fully
described by the possible interpretations proposed above. It then follows that
the above-proposed tests are sufficient to fully specify the considered
problem given the assumptions upon which it was constructed.

##  3 Results

###  3.1 Performance of RW Structure Detection

When assuming a normally-distributed jump structure (23), we propose that,
given sufficiently large blockwise sample size $N$, the controlled variations
test, as defined in (26), asymptotically generates a family of analytical ROC
curves ($\mathbb{P}_{D}=\mathbb{P}_{D}(\alpha)$, $\alpha\in(0,1)$):

$\displaystyle\mathbb{P}_{D}$
$\displaystyle=\Phi\left(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\text{--}1)(1\text{--}\frac{2}{9(K-2)})}{\sqrt{\frac{2}{9(K-2)}}}\right),$
(33)

on the parameter space
($J$,$N$,$K$,$T$,SNR)$\in$$\mathbb{N}^{3}$$\times$$\mathbb{R}^{2}_{>0}$.

###### Proof.

Under $\mathcal{H}_{1},$ the probability distribution
$\mathcal{L}(y^{(k)}_{c})$ of $\\{y^{(k)}_{c}\\}_{k}$ is a perturbation of the
probability distribution $\mathcal{L}(u^{(k)})$ of the process
$\\{u^{(k)}\\}_{k}$ (scaled by a factor as per (7)), perturbed by the
uniformly-distributed variables $\\{\delta_{0}^{(k)}\\}_{k}$. Furthermore,
under $\mathcal{H}_{1}$, $\mathcal{L}(\hat{y}^{(k)}_{c})$ is a perturbation of
$\mathcal{L}(y^{(k)}_{c})$ by normally-distributed variables
$\\{\varepsilon^{(k)}\\}_{k}$ Stoica et al. (2005).

Then, by local asymptotic normality LeCam (1960); O’Neill (2014), the Degrees
of Freedom estimator $\hat{\text{DoF}}$ should converge to the Gaussian case
$\text{DoF}=K-2,$ for sufficiently large blockwise sample sizes $N$, if the
degree of perturbation is not too large. So, we assert that
$\hat{\text{DoF}}\underset{N\to\infty}{\longrightarrow}K-2,$ and so,
$\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}\hat{\chi}_{0}^{2}\underset{N\to\infty}{\longrightarrow}\hat{\chi}_{1}^{2}$
where,

$\displaystyle\hat{\chi}_{1}^{2}=(K-2)\frac{\hat{\sigma}^{2}}{\sigma_{1}^{2}}\sim\chi^{2}(K-2).$
(34)

Two constructions for characterising the jump structure of the RW-like
structure on the pivot frequencies $\\{f_{0}^{(k)}\\}_{k}$ are considered for
the above modelling. We demonstrated the following convergence properties
under simulations.

As shown in Figure 4a,b, the uniform jump model (22) for the RW-like structure
does not result in a convergence of
$\hat{\text{DoF}}\underset{N\to\infty}{\longrightarrow}K-2,$ whereas the
normal jump model (23) for the RW-like structure does (see Figure 4c,d).

Consider a ‘slow’ branching process generated by normally-distributed
variables for the jumps characterising the RW-like structure. Assuming normal
jumps (23), the probability of detection, $\mathbb{P}_{D}$, may then be
determined. By applying the Wilson-Hilferty transform to (34), one can relate
the test statistic (26) with respect to the distribution under
$\mathcal{H}_{1}$ by:

$\displaystyle\hat{Z}_{1}=\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\hat{Z}_{0}+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9(K-2)})}{\sqrt{\frac{2}{9(K-2)}}}\overset{\text{a}}{\sim}\mathcal{N}(0,1).$
(35)

Through algebraic manipulation (see Appendix B) of the test statistic
$\hat{Z}_{0}$, this generates a family of (asymptotic) analytical ROC curves
for the parameters ($\alpha$, $J$, $N$, $K$, $T$, SNR):

$\displaystyle\mathbb{P}_{D}$
$\displaystyle\underset{N\to\infty}{\longrightarrow}\Phi\left(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9(K-2)})}{\sqrt{\frac{2}{9(K-2)}}}\right).$
(36)

∎

Given the analytical performance limits as presented in Proposition 3.1, we
then ask for what sample sizes and signal-to-noise ratios such performance
bounds reflect the performance of the model detection methodology considered.

Figure 7a,b reflect the rate of convergence of the performance with respect to
increasing blockwise sample sizes, given a signal-to-noise ratio.

As in Figure 7a,b specifically, the analytical performance bounds presented in
Proposition 3.1 represent the performance of the detectors, in practice, in
the discrete sample regime ($N\leq 128$) for signal-to-noise ratios above
$-10$ dB.

However, simulations failed to provide results for signal-to-noise ratios of
$-20$ dB or less, the underlying reason likely due to non-optimal estimation
procedures used (as this was just a proof of concept for the verification
method). Further, the simulated ROC curves did not appear to converge to their
theoretical counterparts in the very-low SNR regime. Then, it would follow
that it is likely that the analytical performance bounds presented in
Proposition 3.1 fail to represent the performance of the detectors for signal-
to-noise ratios of $-20$ dB or less.

Now, the Chow-Denning test statistic $V_{1}$ (27) of size $J$ constructed on
the pivot frequencies across a recording of $K$ blocks follows a Studentised
Maximum Modulus $SMM(\mu,\nu)$ distribution with parameters $\mu=J$, $\nu=K-1$
Chow and Denning (1993); Charles and Darné (2009).

Following the same logic that the Chow-Denning test statistic assumes the
independence of the LOMAC test statistics from which it was constructed Chow
and Denning (1993), assume the independence of the Chow-Denning test statistic
$V_{1}$ (27) and the Controlled Variations test statistic $\hat{Z}_{0}$ (26).
Given this assumption, we may define the detection distribution of the joint
tests, under appropriate Šidák corrections to the false alarm rate, as the
product distribution
$\mathbb{P}_{D}^{*}:=\mathbb{P}(\hat{Z}_{0}<\gamma_{1},V_{1}>\gamma_{2}|\mathcal{H}_{1})$:

$\displaystyle\mathbb{P}_{D}^{*}(\alpha)$
$\displaystyle=\mathbb{P}(\hat{Z}_{0}<\gamma_{1}|\mathcal{H}_{1})\mathbb{P}(V_{1}>\gamma_{2}|\mathcal{H}_{1}),$
$\displaystyle=\mathbb{P}_{D}(\alpha)\mathbb{P}\left(V_{1}>\Phi^{-1}\left(1-\frac{\alpha^{*}}{2}\right)\right),$
(37)

recalling that, as mentioned above, $V_{1}\sim SMM(J,K-1)$.

##  4 Discussion

When performing a model verification test as shown above, to a specified
$\alpha\in(0,1),$ the analytical performance bounds (as defined in (33), and
illustrated in Figures 5a and 6) allow one to optimally partition the data
across blocks, to better prime the data for analysis, by the following
procedure:

  * •

Determine the sampling rate $f_{s}$ given the Nyquist frequency of the
hypothesized signal;

  * •

Given an SNR$>-20$ dB, fix $N$ $\in$ [$N_{undersampling}$, $N_{oversampling}$]
for which the estimator converges to within a neighborhood less than the
discretisation;

  * •

Take minimal $K\geq 11$ (while maximising
$J\in\\{1,\ldots,\lfloor\frac{K}{2}\rfloor-1\\}$) such that $\mathbb{P}_{D}$
is optimal.

If the signal-to-noise ratio, however, is not greater than $-20$ dB, then a
different estimation procedure would need to be considered.

However, as mentioned in Section 2, the drawback is that the detection method
only works given a few key assumptions:

  * •

The blockwise phases $\\{\phi_{n}^{(k)}\\}_{n}$ are quadratic;

  * •

The model on which the jumps generate the RW-like structure must be standard
normal;

  * •

The period in which the blockwise frequency modulation effectively occurs is
known. That is to say, the expected value of
$\\{t_{0}^{(k+1)}-t_{0}^{(k)}\\}_{k}$ is known; and,

  * •

The pivot frequencies $\\{f_{0}^{(k)}\\}_{k}$ are uniformly spaced, that
var[$\\{t_{0}^{(k+1)}-t_{0}^{(k)}\\}_{k}$] is sufficiently small.

The temporal location of (at least) one of the pivots $f_{0}^{(k)}$ must be
known with sufficient accuracy to align the record with the structure in the
derivation. That is, the pivot frequencies do indeed occur at
$\\{t_{0}^{(k)}\\}_{k}$.

However, note that, while the degree of freedom estimator did not converge to
the theoretical values under the discrete uniform jump model, the ROC curves
still provide intuition on how to partition the data for analysis.

##  5 Conclusions

A model verification test has been designed, which has been shown in
simulations to operate as anticipated. The test generalises the temporal RWH
testing from the economics literature to test for an RW in frequency, by
introducing an additional test to account for phase-wrapping behaviours, and
how the peak frequencies of sampled WGN behave.

The use of such a model verification test should reduce the need for more
heuristic forms of vetting false-positive detections.

We have yet to explore the impact of biases (such as a mean-reverting
behaviour as exhibited under $\mathcal{H}_{2}$) on the detection procedure
above. Further simulations and modelling are needed to understand how biases
affect the ability to determine the presence of RW-like structures in tonal
variations.

The above modelling is done for the case where pivots in frequency space are
discrete in time. It would be worthwhile to extend the above model
verification procedure for Brownian motion-like (BM) structures, with
continuous randomness in tonal variations.

Additionally, the above derivations are only formulated for linear blockwise
frequency modulations. It would not be complicated to generalise to polynomial
blockwise frequency modulations to refine the model verification problem to a
SID problem with respect to model order.

Lastly, one normally considers the apparent randomness of behaviours in
periodic parameters (such as tones) when interested in forecasting (estimating
and/or detecting) the potential of a trend in the signal. That is, for the
model presented on the underlying object creating the perceived phenomenon,
does the data admit a long-term trend?

A key issue in this problem would be the large number of similar models that
could be determined to occur under the same case from the hypothesis structure
explored in this paper.

For the following paper, we intend to look into how the question of trend-
based behaviours could come into each of the outcomes of the hypothesis test,
how one would distinguish between these different trend-based behaviours, and
test whether the data support such models.

Conceptualization, R.E., W.M. and A.M.; methodology, S.A.; validation, R.E.
and W.M.; formal analysis, S.A.; investigation, S.A.; writing—original draft
preparation, S.A.; writing—review and editing, S.A., R.E., W.M. and A.M.;
supervision, R.E., W.M. and A.M.; project administration, S.A. All authors
have read and agreed to the published version of the manuscript.

The authors would like to acknowledge that this project would not have been
possible without the funding provided by the Australian Commonwealth
Government and the University of Melbourne. S. Abramson is supported by an
Australian Government Research Training Program Scholarship. This work was
supported by the Australian Research Council (ARC) Centre of Excellence for
Gravitational Wave Discovery (OzGrav), grant number CE170100004, and an ARC
Discovery project, grant number DP170103625.

The authors declare no conflict of interest. The funders had no role in the
design of the study; in the collection, analyses, or interpretation of data;
in the writing of the manuscript, or in the decision to publish the results.

yes

##  Appendix A First Difference Variance Derivations

Under $\mathcal{H}_{0}$, the carrier frequency estimators
$\left\\{\hat{f}^{(k)}_{c}\right\\}_{k}$ are distributed according to (18),
and so, taking the first two moments yields:

$\displaystyle\mathbb{E}\left[\hat{f}^{(k)}_{c}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\frac{1}{N}\sum_{n=0}^{\frac{N}{2}-1}\frac{n}{NT}+\frac{1}{N}\sum_{n=\frac{N}{2}}^{N-1}\frac{n-N}{NT}$
$\displaystyle=-\frac{1}{2NT},$ (38)
$\displaystyle\mathbb{E}\left[\left(\hat{f}^{(k)}_{c}\right)^{2}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\frac{1}{N}\sum_{n=0}^{\frac{N}{2}-1}\left(\frac{n}{NT}\right)^{2}+\frac{1}{N}\sum_{n=\frac{N}{2}}^{N-1}\left(\frac{n-N}{NT}\right)^{2}$
$\displaystyle=\frac{1}{12}\left[\left(\frac{1}{T}\right)^{2}+2\left(\frac{1}{NT}\right)^{2}\right].$
(39)

Thus, the variance, under $\mathcal{H}_{0}$, follows (38) and (39) as:

$\displaystyle\mathbb{V}\left[\hat{f}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{E}\left[\left(\hat{f}_{c}^{(k)}-\mathbb{E}\left[\hat{f}_{c}^{(k)}\right]\right)^{2}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{E}\left[\left(\hat{f}_{c}^{(k)}\right)^{2}\middle|\mathcal{H}_{0}\right]-\left(\mathbb{E}\left[\hat{f}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]\right)^{2}$
$\displaystyle=\frac{1}{12}\left[\left(\frac{1}{T}\right)^{2}-\left(\frac{1}{NT}\right)^{2}\right].$
(40)

Given that under $\mathcal{H}_{0}$ the central frequency estimators
$\left\\{\hat{f}_{c}^{(k)}\right\\}_{k}$ are serially independent, then,
following the algebra of variances, the first difference estimator’s variance
under $\mathcal{H}_{0}$ is found to be:

$\displaystyle\sigma_{0}^{2}:$
$\displaystyle=\mathbb{V}\left[\hat{y}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{V}\left[\left(\hat{f}_{c}^{(k)}-\hat{f}_{c}^{(k-1)}\right)\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{V}\left[\hat{f}_{c}^{(k)}\middle|\mathcal{H}_{0}\right]+\mathbb{V}\left[\hat{f}_{c}^{(k-1)}\middle|\mathcal{H}_{0}\right]$
$\displaystyle=\frac{1}{6}\left[\left(\frac{1}{T}\right)^{2}-\left(\frac{1}{NT}\right)^{2}\right].$
(41)

Next, under $\mathcal{H}_{1}$, $\mathcal{H}_{2}$, or $\mathcal{H}_{3}$, by
combining (8), (11), (13) and (16), the first difference estimator process
$\left\\{\hat{y}_{c}^{(k)}\right\\}_{k}$ has representation:

$\displaystyle\hat{y}^{(k)}_{c}=$
$\displaystyle\frac{1}{2}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\right]$
$\displaystyle+\frac{1}{2}\left[\delta^{(k+1)}_{0}-\delta^{(k-1)}_{0}\right]$
$\displaystyle\ \ +\left[\varepsilon^{(k)}-\varepsilon^{(k-1)}\right],$ (42)

where the discretisation uncertainty terms
$\left\\{\delta_{0}^{(k)}\right\\}_{k},$ estimation errors
$\left\\{\varepsilon^{(k)}\right\\}_{k},$ and pivot Fourier cell estimators
$\left\\{\hat{m}_{0}^{(k)}\right\\}_{k},$ are all mutually independent
processes.

Noting that the effective frequency of a Fourier cell $f(\cdot),$ as defined
by (7), is a linear transformation on the indices, then the independence of
the pivot Fourier cell estimators $\left\\{\hat{m}_{0}^{(k)}\right\\}_{k},$
from the two ‘noise’ processes $\left\\{\delta_{0}^{(k)}\right\\}_{k},$
$\left\\{\varepsilon^{(k)}\right\\}_{k},$ commutes to the process
$\left\\{f\left(\hat{m}_{0}^{(k)}\right)\right\\}_{k}$ of effective pivot
frequency estimators.

Following the serial independence of the two ‘noise’ processes
$\left\\{\delta_{0}^{(k)}\right\\}_{k},$
$\left\\{\varepsilon^{(k)}\right\\}_{k},$ the algebra of variances, and
representation (42) we take the variance, not under $\mathcal{H}_{0},$ to be:

$\displaystyle\mathbb{V}\left[\hat{y}^{(k)}_{c}\middle|\text{ not
}\mathcal{H}_{0}\right]=$
$\displaystyle\frac{1}{4}\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\text{
not }\mathcal{H}_{0}\right]$
$\displaystyle+\frac{1}{4}\mathbb{V}\left[\delta^{(k+1)}_{0}-\delta^{(k-1)}_{0}\middle|\text{
not }\mathcal{H}_{0}\right]$ $\displaystyle\ \
+\mathbb{V}\left[\varepsilon^{(k)}-\varepsilon^{(k-1)}\middle|\text{ not
}\mathcal{H}_{0}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{4}\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\text{
not }\mathcal{H}_{0}\right]$
$\displaystyle+\frac{1}{2}\mathbb{V}\left[\delta^{(k)}_{0}\middle|\text{ not
}\mathcal{H}_{0}\right]+2\mathbb{V}\left[\varepsilon^{(k)}\middle|\text{ not
}\mathcal{H}_{0}\right].$

Now, given the distributions of the two ‘noise’ processes
$\left\\{\delta_{0}^{(k)}\right\\}_{k},$
$\left\\{\varepsilon^{(k)}\right\\}_{k},$ we have:

$\displaystyle\mathbb{V}\left[\delta^{(k)}_{0}\middle|\text{ not
}\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{V}\left[\delta^{(k)}_{0}\right]$
$\displaystyle=\frac{1}{12}\left(\frac{1}{NT}\right)^{2},$
$\displaystyle\mathbb{V}\left[\varepsilon^{(k)}\middle|\text{ not
}\mathcal{H}_{0}\right]$
$\displaystyle=\mathbb{V}\left[\varepsilon^{(k)}\right]$
$\displaystyle=\frac{1}{N}\frac{6}{(2\pi NT)^{2}\text{SNR}}.$

Combining the above values to the variance of the first central difference
when not under $\mathcal{H}_{0}$ yields:

$\displaystyle\mathbb{V}\left[\hat{y}^{(k)}_{c}\middle|\text{ not
}\mathcal{H}_{0}\right]$
$\displaystyle=\frac{1}{4}\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\text{
not }\mathcal{H}_{0}\right]$
$\displaystyle+\frac{1}{24}\left(\frac{1}{NT}\right)^{2}+\frac{1}{N}\frac{12}{(2\pi
NT)^{2}\text{SNR}}.$ (43)

Now, specifically under $\mathcal{H}_{1},$ we assume that we may use a
representation of the form:

$\displaystyle m_{0}^{(k)}$ $\displaystyle=m_{0}^{(k-1)}+u^{(k)},$ (44)

where we define the discretised driving term of the RW-like structure, as a
sequence of integer-valued jumps between Fourier bins, represented by
$\\{u^{(k)}\\}_{k}$.

Now, assuming we do not cross a boundary in Fourier space, on a given step
range $m_{0}^{(k-1)}$ to $m_{0}^{(k+1)},$ then,

$f\left(m_{0}^{(k+1)}\right)-f\left(m_{0}^{(k-1)}\right)=f\left(m_{0}^{(k+1)}-m_{0}^{(k-1)}\right).$

For large $N,$ under the assumption of a ‘slow’ wander (6), on average, a
boundary in Fourier space is not crossed across two steps. Thus, it is
reasonable to assume that:

$\mathbb{E}\left[f\left(m_{0}^{(k+1)}\right)-f\left(m_{0}^{(k-1)}\right)\right]=\mathbb{E}\left[f\left(m_{0}^{(k+1)}-m_{0}^{(k-1)}\right)\right].$

Given the linearity of expectations,
$\mathbb{E}\left[f\left(\cdot\right)\right]=f\left(\mathbb{E}\left[\cdot\right]\right).$
Utilising this linearity, and the representation (44),

$\displaystyle\mathbb{E}\left[f\left(\hat{m}_{0}^{(k+1)}\right)-f\left(\hat{m}_{0}^{(k-1)}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\mathbb{E}\left[f\left(\hat{m}_{0}^{(k+1)}-\hat{m}_{0}^{(k-1)}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\mathbb{E}\left[f\left(\hat{m}_{0}^{(k-1)}+u^{(k+1)}+u^{(k)}-\hat{m}_{0}^{(k-1)}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\mathbb{E}\left[f\left(u^{(k+1)}+u^{(k)}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=f\left(\mathbb{E}\left[u^{(k+1)}+u^{(k)}\middle|\mathcal{H}_{1}\right]\right).$

By definitions considered of $\left\\{u^{(k)}\right\\}_{k},$ it has zero
expectation under any structure assumed, which then yields above that:

$\displaystyle\mathbb{E}\left[f\left(\hat{m}_{0}^{(k+1)}\right)-f\left(\hat{m}_{0}^{(k-1)}\right)\middle|\mathcal{H}_{1}\right]=0.$

Next, following the above result, the variance of this difference is merely
its second moment. Given the above assumptions:

$\displaystyle\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\mathbb{E}\left[\left(f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\right)^{2}\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\mathbb{E}\left[\left(f\left(\hat{m}^{(k+1)}_{0}-\hat{m}^{(k-1)}_{0}\right)\right)^{2}\middle|\mathcal{H}_{1}\right].$

Now, for any two arbitrary (but fixed) cells $m_{1},$ $m_{2}$, that satisfy
the assumption (6), $f\left(m_{2}-m_{1}\right)=\frac{1}{NT}(m_{2}-m_{1}).$

Thus, under the above assumptions:

$\displaystyle\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\left(\frac{1}{NT}\right)^{2}\mathbb{E}\left[\left(u^{(k+1)}+u^{(k)}\right)^{2}\middle|\mathcal{H}_{1}\right].$

So, by algebra of variances, and the serial independence of the jump process
$\left\\{u^{(k)}\right\\}_{k}$,

$\displaystyle\mathbb{V}\left[f\left(\hat{m}^{(k+1)}_{0}\right)-f\left(\hat{m}^{(k-1)}_{0}\right)\middle|\mathcal{H}_{1}\right]$
$\displaystyle=2\left(\frac{1}{NT}\right)^{2}\mathbb{E}\left[\left(u^{(k+1)}\right)^{2}\middle|\mathcal{H}_{1}\right].$

For distributions considered for the jump process
$\left\\{u^{(k)}\right\\}_{k}$,

$\displaystyle\mathbb{E}\left[\left(u^{(k+1)}\right)^{2}\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\begin{cases}\frac{2}{3},\quad\text{ for uniform jumps},\\\
\frac{13}{12},\quad\text{ for normal jumps}.\end{cases}$

Thus,

$\displaystyle\sigma_{1}^{2}$
$\displaystyle:=\mathbb{V}\left[\hat{y}_{c}^{(k)}\middle|\mathcal{H}_{1}\right]$
$\displaystyle=\frac{1}{N}\frac{12}{(2\pi
NT)^{2}\text{SNR}}+\begin{cases}\frac{9}{24}(\frac{1}{NT})^{2},\quad\text{ for
uniform jumps},\\\ \frac{14}{24}(\frac{1}{NT})^{2},\quad\text{ for normal
jumps}.\end{cases}$ (45)

##  Appendix B Probability of Detection Derivation

To verify $\mathcal{H}_{0},$ at a false alarm rate $\alpha,$ we check
$\hat{Z}_{0}$ against some threshold $\gamma$ such that if
$\hat{Z}_{0}\geq\gamma$ we argue that $\mathcal{H}_{0}$ is statistically
significant to false alarm rate $\alpha.$

Given that under $\mathcal{H}_{0}$,
$\hat{Z}_{0}\overset{\text{a}}{\sim}\mathcal{N}(0,1),$ we set the threshold to
be $\gamma=\Phi^{-1}(\alpha^{*}),$ where $\Phi(\cdot)$ is the standard normal
CDF, and $\alpha^{*}=1-(1-\alpha)^{1/(J+1)}$ is the Šidák correction for the
joint inference of a Chow-Denning Test (27) of size J and the Controlled
Variations Test (26).

However, given
$\hat{\chi}^{2}_{1}=\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}\hat{\chi}^{2}_{0},$
under $\mathcal{H}_{1}$,

$\displaystyle\hat{Z}_{1}$
$\displaystyle=\frac{\sqrt[3]{\frac{\hat{\chi}^{2}_{1}}{\hat{\text{DOF}}}}-(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}\overset{\text{a}}{\sim}\mathcal{N}(0,1),$

Now,

$\displaystyle\hat{Z}_{0}$
$\displaystyle=\frac{\sqrt[3]{\frac{\hat{\chi}^{2}_{0}}{\hat{\text{DOF}}}}-(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}.$

So, by algebraic manipulation, we may relate $\hat{Z}_{0},$ and $\hat{Z}_{1}$:

$\displaystyle\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\hat{Z}_{0}$
$\displaystyle=\frac{\sqrt[3]{\frac{\hat{\chi}^{2}_{1}}{\hat{\text{DOF}}}}-\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}$
$\displaystyle=\frac{\sqrt[3]{\frac{\hat{\chi}^{2}_{1}}{\hat{\text{DOF}}}}}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}-\frac{\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}$
$\displaystyle+\frac{(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}-\frac{(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}$
$\displaystyle=\frac{\sqrt[3]{\frac{\hat{\chi}^{2}_{1}}{\hat{\text{DOF}}}}-(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}$
$\displaystyle-\frac{\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}(1-\frac{2}{9\hat{\text{DOF}}})-(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}$
$\displaystyle=\hat{Z}_{1}-\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}.$

Now, we consider the probability of detection of $\mathcal{H}_{1}$ against
$\mathcal{H}_{0}$, that is that $\hat{Z}_{0}<\gamma$ given that
$\mathcal{H}_{1}$ holds. Mathematically:

$\displaystyle\mathbb{P}_{D}$
$\displaystyle:=\mathbb{P}(\hat{Z}_{0}<\gamma|\mathcal{H}_{1})$
$\displaystyle=\mathbb{P}(\hat{Z}_{0}<\Phi^{-1}(\alpha^{*})|\mathcal{H}_{1})$
$\displaystyle=\mathbb{P}\left(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\hat{Z}_{0}<\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})\middle|\mathcal{H}_{1}\right)$
$\displaystyle=\mathbb{P}\left(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\hat{Z}_{0}+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}\right.$
$\displaystyle\left.\quad\quad\quad<\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}\middle|\mathcal{H}_{1}\right)$
$\displaystyle=\mathbb{P}\left(\hat{Z}_{1}<\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9\hat{\text{DOF}}})}{\sqrt{\frac{2}{9\hat{\text{DOF}}}}}\middle|\mathcal{H}_{1}\right).$

Then, in such a case we have a family of detection curves that take the form:

$\mathbb{P}_{D}=\mathbb{P}_{D}(\alpha^{*},\sigma_{0}^{2},\sigma_{1}^{2},K).$

However, if assuming a normal jump structure:

$\hat{\text{DOF}}\underset{N\to\infty}{\longrightarrow}K-2.$

Given that under $\mathcal{H}_{1}$,
$\hat{Z}_{1}\overset{\text{a}}{\sim}\mathcal{N}(0,1),$ if assuming a normal
jump structure, we would then have that:

$\displaystyle\mathbb{P}_{D}$
$\displaystyle\underset{N\to\infty}{\longrightarrow}\Phi\left(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}\Phi^{-1}(\alpha^{*})+\frac{(\sqrt[3]{\frac{\sigma_{0}^{2}}{\sigma_{1}^{2}}}-1)(1-\frac{2}{9(K-2)})}{\sqrt{\frac{2}{9(K-2)}}}\right).$
(46)

Lastly, recalling that for the case of a normal jump structure, the first
difference variances take the form:

$\displaystyle\sigma_{0}^{2}$
$\displaystyle=\frac{1}{6}\left[\left(\frac{1}{T}\right)^{2}-\left(\frac{1}{NT}\right)^{2}\right],$
$\displaystyle\sigma_{1}^{2}$ $\displaystyle=\frac{1}{N}\frac{12}{(2\pi
NT)^{2}\text{SNR}}+\frac{14}{24}\left(\frac{1}{NT}\right)^{2}.$

Thus, noting that $\alpha^{*}=1-(1-\alpha)^{1/(J+1)},$ the detection curves
are parameterised by:

$\mathbb{P}_{D}=\mathbb{P}_{D}(\alpha,J,N,T,\text{SNR},K).$

##  Appendix C Figures

a

b

Figure 2: Estimation Error. (a) Estimation Error, as % of total cells $N$,
versus Random Walk Lengths $K$. (b) Estimation Error, as % of total cells $N$,
versus Intersampling Periods $T$.

a

b

Figure 3: Cont.

a

b

Figure 4: Degrees of Freedom Estimator. (a) Degrees of Freedom Estimator
$\hat{\text{DoF}}$, Blockwise Sample Sizes $N$, SNR = $0$ dB, $T$ = $1$ s;
uniform jumps. (b) Degrees of Freedom Estimator $\hat{\text{DoF}}$, SNR, $N$ =
$64$, $T$ = $1$ s; uniform jumps. (c) Degrees of Freedom Estimator
$\hat{\text{DoF}}$, Blockwise Sample Sizes $N$, SNR = $0$ dB, $T$ = $1$ s;
normal jumps. (d) Degrees of Freedom Estimator $\hat{\text{DoF}}$ SNR, $N$ =
$64$, $T$ = $1$ s; normal jumps.

a

b

c

d

Figure 5: ROC Curves. (a) ROC Curves for Blockwise Sample Sizes $N$, $J=7$,
$K=16$, $T=0.1$ s, SNR = $-20$ dB. (b) ROC Curves for Random Walk Lengths $K$,
$J=7$, $N=16$, $T=0.1$ s, SNR = $-20$ dB. (c) ROC Curves for Intersampling
Periods $T$, $J=7$, $N=16$, $K=16$, SNR = $-20$ dB. (d) ROC Curves for Signal-
to-Noise Ratios SNR, $J=7$, $N=16$, $K=16$, $T=0.1$ s. Figure 6: ROC Curves
for Chow-Denning Test Sizes $J$, $N=16$, $K=16$, $T=0.1$ s, SNR = $-20$ dB.

a

b

Figure 7: ROC Convergence Results. (a) ROC Curve Convergence for SNR = $0$ dB,
$J=7$, $K=16$, $T=0.1$ s. (b) ROC Curve Convergence for SNR = $-2.5$ dB,
$J=7$, $K=16$, $T=0.1$ s.

References

## References

  * Knuth et al. (2015) Knuth, K.H.; Habeck, M.; Malakar, N.K.; Mubeen, A.M.; Placek, B.  Bayesian evidence and model selection.  Digit. Signal Process. 2015, 47, 50–67.
  * Ding et al. (2018) Ding, J.; Tarokh, V.; Yang, Y.  Model selection techniques: An overview.  IEEE Signal Process. Mag. 2018, 35, 16–34.
  * Llorente et al. (2020) Llorente, F.; Martino, L.; Delgado, D.; Lopez-Santiago, J.  Marginal likelihood computation for model selection and hypothesis testing: An extensive review.  arXiv 2020, arXiv:2005.08334.
  * Hui (2012) Hui, T.  Testing for Random Walk Hypothesis with or without Measurement Error.  Master’s Thesis, Imperial College London, London, UK, 2012.
  * Bellman and Åström (1970) Bellman, R.; Åström, K.J.  On structural identifiability.  Math. Biosci. 1970, 7, 329–339.
  * Akaike (1973) Akaike, H.  Information theory and an extension of the maximum likelihood principle. In Proceedings of the 2nd International Symposium on Information Theory, Tsahkadsor, Armenia, 2–8 September 1971;  Akademiai Kiado: Budapest, Hungary, 1973.
  * Rissanen (1982) Rissanen, J.  Estimation of structure by minimum description length.  Circuits Syst. Signal Process. 1982, 1, 395–406.
  * Willems (2007) Willems, J.C.  The behavioral approach to open and interconnected systems.  IEEE Control Syst. Mag. 2007, 27, 46–99.
  * Hofmann et al. (2008) Hofmann, T.; Schölkopf, B.; Smola, A.J.  Kernel methods in machine learning.  Ann. Stat. 2008, _36_ , 1171–1220.
  * Rife and Boorstyn (1974) Rife, D.; Boorstyn, R.  Single tone parameter estimation from discrete-time observations.  IEEE Trans. Inf. Theory 1974, 20, 591–598.
  * Slepian (1954) Slepian, D.  Estimation of signal parameters in the presence of noise.  Trans. IRE Prof. Group Inf. Theory 1954, 3, 68–89.
  * Charles and Darné (2009) Charles, A.; Darné, O.  Variance-ratio tests of random walk: An overview.  J. Econ. Surv. 2009, 23, 503–527.
  * Lo and MacKinlay (1988) Lo, A.W.; MacKinlay, A.C.  Stock market prices do not follow random walks: Evidence from a simple specification test.  Rev. Financ. Stud. 1988, 1, 41–66.
  * Chow and Denning (1993) Chow, K.V.; Denning, K.C.  A simple multiple variance ratio test.  J. Econom. 1993, 58, 385–401.
  * Dickey and Fuller (1981) Dickey, D.A.; Fuller, W.A.  Likelihood ratio statistics for autoregressive time series with a unit root.  Econom. J. Econom. Soc. 1981, _49_ , 1057–1072.
  * Hoque et al. (2007) Hoque, H.A.; Kim, J.H.; Pyun, C.S.  A comparison of variance ratio tests of random walk: A case of Asian emerging stock markets.  Int. Rev. Econ. Financ. 2007, 16, 488–502.
  * Deo and Richardson (2003) Deo, R.S.; Richardson, M.  On the asymptotic power of the variance ratio test.  Econom. Theory 2003, _19_ , 231–239.
  * Daniel (2001) Daniel, K.  The power and size of mean reversion tests.  J. Empir. Financ. 2001, 8, 493–535.
  * Granger and Hallman (1991) Granger, C.W.; Hallman, J.  Nonlinear transformations of integrated time series.  J. Time Ser. Anal. 1991, 12, 207–224.
  * Wang and De Jong (2013) Wang, C.H.; De Jong, R.M.  Unit root tests when the data are a trigonometric transformation of an integrated process.  S. Afr. Stat. J. 2013, 47, 83–90.
  * Ermini and Granger (1993) Ermini, L.; Granger, C.W.  Some generalizations on the algebra of I (1) processes.  J. Econom. 1993, 58, 369–384.
  * Misner et al. (1973) Misner, C.W.; Thorne, K.S.; Wheeler, J.A.  Gravitation; Macmillan: New York, NY, USA, 1973.
  * Aasi et al. (2015) Aasi, J.; Abbott, B.; Abbott, R.; Abbott, T.; Abernathy, M.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.; et al.  Advanced ligo.  Class. Quantum Gravity 2015, 32, 074001.
  * Abbott et al. (2016) Abbott, B.P.; Abbott, R.; Abbott, T.; Abernathy, M.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.; et al.  Observation of gravitational waves from a binary black hole merger.  Phys. Rev. Lett. 2016, 116, 061102.
  * Abbott et al. (2017) Abbott, B.P.; Abbott, R.; Abbott, T.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.; Adya, V.; et al.  GW170817: Observation of gravitational waves from a binary neutron star inspiral.  Phys. Rev. Lett. 2017, 119, 161101.
  * Riles (2013) Riles, K.  Gravitational waves: Sources, detectors and searches.  Prog. Part. Nucl. Phys. 2013, 68, 1–54.
  * Jaranowski et al. (1998) Jaranowski, P.; Krolak, A.; Schutz, B.  Data analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection.  Phys. Rev. D 1998, 58, 063001.
  * Shapiro and Teukolsky (1983) Shapiro, S.; Teukolsky, S.  Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects; John Wiley & Sons: Hoboken, NJ, USA, 1983.
  * de Kool and Anzer (1993) de Kool, M.; Anzer, U.  A simple analysis of period noise in binary X-ray pulsars.  Mon. Not. R. Astron. Soc. 1993, 262, 726–734.
  * Romanova et al. (2004) Romanova, M.M.; Ustyugova, G.V.; Koldoba, A.V.; Lovelace, R.V.  The propeller regime of disk accretion to a rapidly rotating magnetized star.  Astrophys. J. Lett. 2004, 616, L151.
  * Baykal et al. (1991) Baykal, A.; Alpar, A.; Kiziloglu, U.  A shot noise model for a two-component neutron star.  Astron. Astrophys. 1991, 252, 664–668.
  * Melatos et al. (2020) Melatos, A.; Dunn, L.; Suvorova, S.; Moran, W.; Evans, R.  Pulsar glitch detection with a hidden Markov model.  arXiv 2020, arXiv:2005.09388.
  * Suvorova et al. (2016) Suvorova, S.; Sun, L.; Melatos, A.; Moran, W.; Evans, R.  Hidden Markov model tracking of continuous gravitational waves from a neutron star with wandering spin.  Phys. Rev. D 2016, 93, 123009.
  * Wu et al. (2018) Wu, B.; Wu, D.; Gao, W.; Song, C.  Time-variant random interval natural frequency analysis of structures.  J. Sound Vib. 2018, 414, 284–298.
  * Dubois and Davy (2007) Dubois, C.; Davy, M.  Joint detection and tracking of time-varying harmonic components: A flexible Bayesian approach.  IEEE Trans. Audio Speech, Lang. Process. 2007, 15, 1283–1295.
  * Durlauf and Phillips (1988) Durlauf, S.N.; Phillips, P.C.  Trends versus random walks in time series analysis.  Econom. J. Econom. Soc. 1988, _56_ , 1333–1354.
  * Nelson and Kang (1981) Nelson, C.R.; Kang, H.  Spurious periodicity in inappropriately detrended time series.  Econom. J. Econom. Soc. 1981, _49_ , 741–751.
  * Durlauf (1991) Durlauf, S.N.  Spectral based testing of the martingale hypothesis.  J. Econom. 1991, 50, 355–376.
  * Niebauer (2013) Niebauer, T.  Analytic signal demodulation of phase-modulated frequency-chirped signals.  Appl. Opt. 2013, 52, 1838–1846.
  * Stoica et al. (2005) Stoica, P.; Moses, R.L.  Spectral Analysis of Signals; Pearson Prentice Hall: Upper Saddle River, NJ, USA, 2005.
  * Stoica and Nehorai (1989) Stoica, P.; Nehorai, A.  Statistical analysis of two nonlinear least-squares estimators of sine-wave parameters in the colored-noise case.  Circuits Syst. Signal Process. 1989, 8, 3–15.
  * Nabeya and Tanaka (1986) Nabeya, S.; Tanaka, K.  Aproximate Distributions of the Periodogram and Related Statistics under Normality.  Econom. Theory 1986, 2, 33–65.
  * O’Neill (2014) O’Neill, B.  Some useful moment results in sampling problems.  Am. Stat. 2014, 68, 282–296.
  * Douillet (2009) Douillet, P.L.  Sampling distribution of the variance.  In Proceedings of the 2009 Winter Simulation Conference (WSC), Austin, TX, USA, 13–16 December 2009; IEEE: Piscataway, NJ, USA, 2009; pp. 403–414.
  * Klebaner (2005) Klebaner, F.C.  Introduction to Stochastic Calculus with Applications; World Scientific Publishing Company: Singapore, 2005.
  * Mukherjee et al. (2018) Mukherjee, A.; Messenger, C.; Riles, K.  Accretion-induced spin-wandering effects on the neutron star in Scorpius X-1: Implications for continuous gravitational wave searches.  Phys. Rev. D 2018, 97, 043016.
  * Bhargava (1986) Bhargava, A.  On the theory of testing for unit roots in observed time series.  Rev. Econ. Stud. 1986, 53, 369–384.
  * Nicolau (2002) Nicolau, J.  Stationary processes that look like random walks—The bounded random walk process in discrete and continuous time.  Econom. Theory 2002, 18, 99–118.
  * Diebold and Rudebusch (1991) Diebold, F.X.; Rudebusch, G.D.  On the power of Dickey-Fuller tests against fractional alternatives.  Econ. Lett. 1991, 35, 155–160.
  * Ulloa and Pesaran (2006) Ulloa, M.R.D.; Pesaran, M.H.  Nonnested hypotheses. In _Econometrics_ ; Springer: Berlin/Heidelberg, Germany, 2006; pp. 167–173.
  * LeCam (1960) LeCam, L.  Locally asymptotically normal families of distributions.  Univ. Calif. Publ. Statist. 1960, 3, 37–98.