Detecting atomic clock frequency trends using an optimal stopping method

We consider the problem of rapid detection of anomalies in atomic clocks. We assume that the clock error evolution can be modeled by a Wiener process with a drift changing from 0 to μ≠0 at some random unobservable time θ and we approach the problem by applying the quickest detection method. We describe the scope and the limits of the method and we apply it to clock experimental data.


Introduction
Recently, particularly in space applications such as the global satellite navigation systems (GNSS), anomalous clock behavior has been experimentally observed: phase or frequency jumps, increase of instability, general non-stationary trends. Such anomalous behavior has an important impact on the global system and may affect its overall performance very seriously, as it compromises the reliability of the atomic clocks which is a crucial point in navigation systems.
When a clock has an anomaly, this has to be detected very quickly and the users have to be informed promptly to avoid the use of a failing clock and hence a wrong position estimate. Rapid detection of anomalies is therefore of the utmost importance together with the need for very low false alarm probability.
Here we approach the problem by applying the quickest detection method (see [1,2]). We assume that the clock deviation evolution can be modeled by a Wiener process ( ) ⩾ = X X t t 0 and that we observe a trajectory with a drift changing from 0 to µ ≠ 0 at some random time θ (see [3,4]). Our task is to find a stopping time τ of the process X that is as close as possible to the unknown time θ.
We describe the scope and the limits of the method applied to the detection of the atomic clock anomalies and we present results on the application of this method to clock experimental data.

Detection method
Suppose that we observe a trajectory of a Wiener process ( ) ⩾ = X X t t 0 with a drift changing from 0 to µ ≠ 0 at some random time θ. We assume that θ is a random variable that takes the value 0 with probability π and it is exponentially distributed with parameter λ > 0 given that θ > 0. We are looking for a stopping time τ of the observed process X that is as close as possible to the unknown time θ.
The observed process X can be described as which is equivalent to where ( ) ⩾ = W W t t 0 is a standard Wiener process, µ ≠ 0 and σ > 0 2 are given and fixed constants, and ( ) ⋅ I is the indicator function.

Detecting atomic clock frequency trends using an optimal stopping method
The aim is to find a stopping time of the process X that is as close as possible to the unknown time θ by avoiding false alarms and minimizing the detection delay. This means that we are looking for a stopping time τ that minimizes the risk function is the probability of false alarm, i.e. the probability of declaring the alarm when it is not needed; ( ) τ θ − = π + E D : is the expected detection delay, i.e. the mean delay to detect the disorder, and c is a positive constant that balances the two quantities. By (z) + we denote the positive part of z (which equals z when positive and 0 otherwise).
It is known (see [1,2]) that the problem (3) can be rewritten as an optimal stopping problem with the risk function t 0 is the natural filtration of the process X, which means the history of X up to time t for ⩾ t 0. Observe that Π t takes values in [0, 1] for all ⩾ t 0. The following theorem (see [1,2]) presents a solution to the optimal stopping problem.
is given by for [ ] π ∈ 0, 1 , where A is the unique root of the equation It is known (see [1,2]) that the a posteriori probability process ( ) ⩾ Π = Π t t 0 can be computed as and (11) Remark 2.1. From (9)-(11) we note that the random variable Π t requires knowledge of the process X up to time t only. This means that the corresponding detection algorithm is a real time algorithm.

Remark 2.2.
We emphasize that, in order to apply the theorem, the following assumptions need to be satisfied: • the underlying observed process is a Wiener process which changes its drift at a random time θ; • the new drift µ and the diffusion coefficient σ of the Wiener process are known; • the disorder time θ is exponentially distributed and its parameter λ is known; • the constant c (or, equivalently, the false alarm probability −A 1 ) is known; • π, i.e. the probability that the disorder time θ has occurred at time zero is known.

The algorithm
Here we explain, step by step, the detection algorithm.
(i) Hypothesize that the process X and its parameters µ, σ, λ, π are given. (ii) Now there are two possible approaches: we can first fix either the parameter c or the parameter A. In the first case, we fix c and compute A as the unique root of the equation (7). In the second case, we fix A and, recalling that cf A is a linear function of c, we find c using (7). (iii) We can then apply the detection algorithm, i.e. we are waiting for the first time τ A such that the a posteriori probability process Π defined in (9) hits the boundary A (assuming that π < A). Note that we can build Π t in (8) using (9)-(11). Moreover, to compute the integral in (10) it is useful to make use of the numerical rectangle method [5].

Remark 2.3.
The second case of step (ii) is more useful because the parameter A has a probabilistic interpretation. Indeed, the quantity 1 − A is the probability of false alarm and often, in applications, it is very important that this probability is sufficiently close to zero. For this reason, in what follows, we will assume that the value of A is given and fixed.

Expected delay and probability of false alarm
In assessing the performance of the method it is important to study the delay in detection, i.e.
It is a random variable whose characteristics are not described in the literature and we note that it is possible to write its expected value D in closed form.
We first need to express the function ψ in closed form. Computing explicitly the integral in (7), we get is the incomplete gamma function [6]. Using the previous result the expected delay can be written as Note that the expected delay depends on the process parameters µ and σ just through the signal-to-noise ratio σ µ / . This quantity is very important for the reliability of the algorithm as shown in the following examples.
Here we analyse the behavior of the expected delay as a function of the probability of false alarm as the parameters λ and σ µ / change. Figures 1 and 2 show the expected delay D as a function of PFA = 1 − A and the parameter λ with µ = 5, σ = 1. We can see that, fixing PFA, as λ decreases, the expected delay increases. Indeed, a small λ means that the expected waiting time for the change of drift to occur is large and it implies a higher expected delay in the detection, as the method is more cautious in declaring an anomaly which is known to be rare. Nevertheless, even with very small PFA, the expected delay is just one time unit. Figure 3 shows the expected delay D as a function of the PFA = 1 − A and the signal-to-noise ratio σ µ / with / λ = 1 360. We can see that, fixing PFA, as σ µ / decreases, the expected delay increases. Actually, a small signal-to-noise ratio makes a disorder time more difficult to detect and gives a longer expected delay. If σ µ > / 1 the detection algorithm is reliable, 1 and gets closer to zero, the expected delay grows significantly and the result becomes less applicable unless the PFA is enlarged as shown in the application to data from the New Horizons Project [7].

Simulations
In order to illustrate the behavior, the performance and the dependence on parameters of the detection algorithm, in this section we present some simulated examples.
We have already pointed out that the algorithm requires the knowledge of parameters µ , σ, λ, π but, in applications, their values are not known and should be estimated: • the parameter λ may be estimated from previous observations of the drift change of the process; • the diffusion coefficient σ can be estimated using the Allan variance [8] or the increments of the process X; 0 is fixed in advance by the application problem. It is often set to zero and it represents the probability that the anomaly occurs at time zero.   • the parameter µ is the critical point in the problem. If the data are analyzed a posteriori it is possible to estimate its value using the sample mean of the increments but the reliability of this estimate depends on the signal-to-noise ratio. Moreover, it is well known from the statistics literature (dating back to 1940-50s) that it is impossible to estimate it properly without a huge number of observations. Running the algorithm in real time we may not know the new drift value in the incoming change and we may need to guess this value from previous experience.
Hereafter we study the robustness of the method, varying µ and σ, through simulations. Figure 4 shows a sample path of the observed process X with parameters µ = 3, σ = 1, π = 0 and the corresponding a posteriori probability process Π with PFA = 1 − A = 0.01. The disorder time is θ = 25 and the algorithm recognizes τ = 26.63 A as the optimal stopping time. From the figure we can see that, as the trajectory of X increases, the trajectory of Π increases towards the threshold A.
One of the main problems in the application of the detection algorithm is the knowledge of the new drift µ. Indeed, this parameter is often not known. Figure 5 shows a simulated sample path of the process X with parameters µ = 3, σ = 1, / λ = 1 360, π = 0 and the corresp onding process Π with A = 0.97 and different guessed values of the drift µ. Table 1 lists the values of τ A and the theoretical expected detection delay for different values of µ.
In this case we are in a favorable situation as the signalto-noise ratio is greater or equal than 1. Indeed, we can see that the disorder time θ = 300 is recognized by the algorithm even if we set a wrong µ in the application of the method. We have to note that we only apply the method with a 'slightly' wrong value of µ, as at least the order of its magnitude should be known.
Next we test the detection algorithm's response when a wrong choice of λ has been made. Figure 6 shows a simulated sample path of the process X with parameters µ = 3, σ = 1, / λ = 1 360, π = 0 and the corresp onding process Π with A = 0.97 and different values of λ. Table 2 lists the values of τ A and the theoretical expected detection delay for different values of λ.
We can see that the disorder time θ = 300 is recognized by the algorithm even if we set a wrong λ.

Real data
In this section we test the algorithm on some real clock data [9]. For application reasons in the following we set PFA = 1 − A = 10 −7 .
The first data set was extracted by a GLONASS space caesium clock. The sample size of the data set is 8064 and the sampling time is 300 s. The time deviation data and the increments are represented in figures 7(a) and (b). Analyzing   the time series both in the time and in the frequency domain (we omit these results) we can conclude that the time deviation X can be approximated by a Wiener process. Looking at figure 7(a), it seems to be apparent that the observed process X has no zero drift before the disorder time. So we estimate the initial drift µ = × − 2.63 10 0 13 and we normalize the observed data by compensating this first trend through the transformation ˜= −µ X X t t t 0 . The resulting sample path of the process X is shown in figure 7(c). We now estimate the diffusion coefficient σ = × − 6.71 10 12 using Allan variance and the new possible trend µ = × − 1.14 10 12 from the increments and we apply the detection algorithm with different values of λ.
In figure 8 the a posteriori probability process Π is plotted for different values of λ. The values of τ A are indicated by a diamond and are very similar. Indeed, it seems that the value of λ does not affect the results significantly. The values of τ A and of the theoretical expected detection delay for different values of λ are listed in table 3. As λ increases, [ ] / θ λ = E 1 decreases, which implies that the disorder time is smaller. Indeed, we can see that the a posteriori probability process Figure 6. Sample path of the process X with parameters µ = 3, σ = 1 and the corresponding a posteriori probability process Π with / λ = 1 10 (blue), / λ = 1 360 (magenta), / λ = 1 1000 (cyan). The disorder time θ is indicated by a star, while the optimal stopping times τ A are indicated by a diamond for 1 − A = 0.03.   Π is more sensitive to the fluctuations of the process X. The smaller the /λ 1 is in size, the earlier the anomaly is expected (see figure 8(d)). Therefore the algorithm is more eager to declare a change.
The second data set was extracted by a GPS space rubidium clock [9]. The sample size of the data set is 186 and the sampling time is 300 s. The time deviation data and the increments are represented in figures 9(a) and (b). Looking at the increments we note that the change of the drift is not instantaneous, in fact the drift takes some time to stabilise. Even if the analysis of the time series in the time and in the frequency domain does not fully confirm that this data set can be approximated by a Wiener process, we nonetheless apply the detection algorithm. Looking at figure 9(a), it seems to be apparent that the observed process X initially has a negative trend. We estimate the initial drift using the sample mean of the increments µ = − × − 1.00 10 0 12 and we normalize the observed data by compensating this first trend. The resulting sample path is shown in figure 9(c). We estimate the diffusion coefficient σ = × − 9.93 10 12 from the Allan variance and the new possible drift µ = × − 1.38 10 12 from the increments and we apply the detection algorithm with different values of λ.
In figure 10 the a posteriori probability process Π is plotted for different values of λ. The values of τ A are indicated by a diamond. The values of τ A and of the theoretical expected detection delay for different values of λ are listed in table 4. Also in this case the optimal stopping times are very similar, and different values of λ do not affect the result. We can conclude that a rough estimate of λ can be sufficient for efficient detection.

Conclusions
The application of the optimal stopping method to the atomic clock data is quite exciting. The method tested on simulated data shows its capability to detect the insurgence of a new drift with a minimal delay even when setting a very low probability of false alarm. A new linear drift in the clock time deviation corresponds to a jump in the frequency values and we know that this can be a serious drawback particularly in GNSS space clocks.
The exposed version of the theorem asks for a considerable knowledge of the process: (i) it shall be a Wiener process; (ii) the parameters µ, λ, σ, and π are assumed to be known.    Regarding the clock time deviation, the assumption of a Wiener process means that its derivative, i.e. the frequency deviation, is white and this is often the case. Also other cases in which the frequency can be modeled by a Wiener process can be found and the method can be applied as well [7] by interpreting the process X as the frequency deviation. The para meter σ can be easily estimated using the Allan deviation, the parameter /λ 1 indicating the mean occurrence of the anomaly can be guessed as well by the knowledge of the clock, and we have seen that the estimates can also be quite inaccurate and the answers of the method are not that different. The only parameter which really represents an issue is the value of the upcoming drift µ. In fact we know that a frequency jump may occur and that this would lead to the insurgence of a different linear slope on the time deviation, but if we want to apply the detector in real time, we cannot know the actual value of the upcoming new drift. We have seen by simulation that the method can cope with a postulated value of the upcoming drift µ that is not exact, however at least the order of its magnitude needs to be determined in some way.
In addition, since real measurement data, particularly from space clocks, are known to be affected by outliers, missing information, and spurious behaviour, all of which corrupt the pure nature of the Wiener process, a suitable preprocessing of data may therefore be necessary.
Nevertheless we have demonstrated a remarkable performance of this method, owing to the existence of an optimal stopping time, which means that under stated conditions, no better estimation is possible. We are currently working on the extension of the method to a random upcoming drift value and studying the analytical properties of the delay in detection to offer a more general method with desirable features in the analysis of atomic clock frequency trends.