Statistical biases and very-long-term time stability analysis

The prediction of very-long-term time stability is a key issue in various fields, such as time keeping, obviously, but also navigation and spatial applications. This is usually performed by extrapolating the measurement data obtained by estimators such as the Allan variance, modified Allan variance, Hadamard variance, etc. This extrapolation may be assessed from a fit over the variance estimates. However, this fit should be performed on the log-log graph of the estimates, which corresponds to a least-squares minimization of the relative difference between the variance estimates and the fitting curve. However, a bias exists between the average of the log of the estimates and the log of the true value of the estimated variance. This paper presents the theoretical calculation of this log-log bias based on the number of equivalent degrees of freedom of the estimates, shows simulations over a large number of realizations, and provides a reliable method of unbiased logarithmic fit. Extrapolating this fit yields a more confident assessment of the very-long-term time stability.

T  stability of clocks and oscillators is generally as- sessed by estimators such as the Allan deviation, modified Allan deviation, Hadamard deviation, etc.However, unlike the variance estimates, the deviation estimates are severely biased and should not be used for further calculations, such as fit.
On the other hand, the prediction of very-long-term time stability is usually performed by fitting and extrapolating an Allan variance (AVAR) curve.Such a curve may be modeled as σ ττττ An extrapolation may be achieved by extending the asymptotes of such a curve beyond the largest τ value (see Fig. 1).
The C i coefficient are related to the noise levels h α of the power law noise model by Table I, [1]: and to the linear frequency drift D 1 : where η(t) represents the pure random component of the frequency deviation y(t).
However such an extrapolation is very sensitive to estimation errors.It is then of importance to perform a reliable estimation of the C i coefficients.The weighted relative least squares are an effective method [2] which is weighted, to take into account the uncertainty ∆ j over each estimate ˆ, σ τ yj 2 () and relative, to be equivalent to a classical least square fit on the log-log plot by minimizing the relative distance between the curve and the estimate [(estimate − curve)/estimate], It is important to emphasize that for the relative fit, i.e., the fit of ∆ ˆσσ 22 / and not of ∆ σ 2 , that 1) it ensures that the fit will minimize the relative differences between the curve and the estimates (e.g., 10% for an estimate at the 10 −14 level as well as 10% for an estimate at the 10 9 level), whereas a classical fit would minimize the absolute differences, overwhelming the small estimates; and 2) because dx/x = d ln x, it is exactly equivalent to a classical fit on the log of the estimates.
As a consequence, this extrapolation method is based on the logarithm of the estimates and not the estimates themselves.We know that the estimate σ τ

A. Statistics of the Allan Variance and the Allan Deviation
The definition of the Allan variance is [3], [4] The estimate of the Allan variance may be calculated as F. Vernotte is with UTINAM, Observatory Terre-Homme-Environnement-Temps-Astronomie (THETA) of Franche-Comté, University of Franche-Comté/Centre National de la Recherche Scientifique (CNRS), Besançon, France (e-mail: francois.vernotte@obs-besancon.fr).
E. Lantz is with Department of Optics, FEMTO-ST, University of Franche-Comté/CNRS, Unité de Formation et de Recherche (UFR) Sciences, Besançon, France.∑ − Because y 2 and y 1 are Gaussian values ).As shown by Fig. 2, the number of EDF reflects the AVAR estimates dispersion.A method for assessing the EDF may be found in [6].

B. Model World and Measurement World
Let us define [see Fig. 3(a)] θ as the model parameter and ξ as the measured data of the parameter.For example, we could consider the true value σ τ Two problems are generally addressed: the direct problem, which attempts to forecast the measurement data knowing the parameter and the inverse problem, which aims to estimate the parameter knowing the measurement data.In the same vein, Tarantola distinguishes the model space, i.e., the space in which the parameter is given, from the data space, i.e., the space in which the measurement data are given [7].In the following, we will use the terms model world and measurement world.

1) Model World (Direct Problem):
In the direct problem, the task is to determine how the measurement data ξ are distributed when the parameter θ 0 is known [see Fig. 3(b)].However, the parameter θ 0 is precisely the unknown quantity that we want to estimate.This parameter is only known in theory and simulations.

2) Measurement World (Inverse Problem):
In the inverse problem, the task is to estimate a confidence interval over θ from the known measurement data ξ 0 [see Fig. 3(c)].This is the usual task of the metrologist.

C. Study of a χ 2 Distribution With ν = 1 Degree of Freedom 1) Properties:
The properties of a χ 1 2 distribution are discussed in this section [8].The probability density function (PDF) is For the high frequency noises ( f +2 FM and f +1 FM), it is necessary to introduce a high cut-off frequency f h [10].
Fig. 1.Extrapolation of a least square fit.The graph shows the Allan deviation (i.e., the square root of the Allan variance) but the fit was performed over the Allan variance estimates according to the criterion (4) using SigmaTheta-1.2[15].The equation of the curve is [(4.7 The PDF is normalized: The mathematical expectation is We can define the logarithmic mean as the exponential of the mathematical expectation of ln (X ) where C is the Euler constant; C ≈ 0.57722.The cumulative distribution function (CDF) is and the inverse CDF is This relationship is useful for getting the bounds of a confidence interval.How is the reduced variable X of the relationships ( 7) to ( 11) related to the measurement data ξ and the model parameter θ ?
2) Parameter, Measurement, and Reduced Variable: Let us consider the standard χ 1 2 variable X = x 2 where x is a Gaussian centered standard random variable: It is well known that σ τ We can then define the reduced variable X as The differential dX is then 3) Conditional Probabilities in the Model World (Direct Problem): The parameter θ is a constant.Let θ = θ 0 ; the differential dX is then: From ( 17) and the relationships ( 7) to ( 9), we find the conditional probability density function: ξ is unbiased: The logarithmic mean is Thus, the mathematical expectation of the log of the measurement data does not converge to the log of the parameter.Therefore, there is a significant bias for log-log representation (and fit); let us call it the log-log bias.

4) Conditional Probabilities in the Measurement World (Inverse Problem):
The measurement data ξ is a constant.Let ξ = ξ 0 ; the differential dX is then From ( 21) and the relationships ( 7)-( 9) comes the conditional probability density function: There is no mean value, In the measurement world also, the mathematical expectation of the log of the parameter does not converge to the log of the measurement data.However, this bias is exactly the inverse of the log-log bias defined previously.(This property will be demonstrated later in a more general context, see Section II-D.)This implies that we just have to renormalize the measurement data by the log-log bias to get both an unbiased estimator in the model world and a mathematical expectation of the parameter equal to the measurement data in the measurement world.

D. Generalization to Any Positive-Valued Random Process
Let ξ be a positive-valued estimator constructed with M measurement data resulting from a random measurement process.ξ is itself a random variable in the model world, which we assume to have a mean.The most obvious example is Allan variance.Another example is the estimation of the standard deviation.It is quite generally possible to define a new reduced variable, which obeys ( 13), unlike X: Let Z = f (Y) be a monotonic function of Y.A fundamental property of density probabilities is where the subscripts emphasize the fact that the different density probabilities are different functions (these subscripts are, as usual, omitted in the following).An immediate consequence is, for Z = ξ, If X is a standard χ M 2 variable with ν = M degrees of freedom, then Y = X/ν ; with ξ a variance measurement data, we have, in the model world, k = 1 and (28) means that ξ is an unbiased estimator of the true variance θ 0 .In the case of standard deviation estimation, the square root of the Allan variance is a biased estimator, but the value of k can be easily calculated as a function of M (the estimator is convergent: k → 1 if M → ∞) and (28) allows also the definition of an unbiased estimator.
The situation is worse in the measurement world.The variable is, in this case, θ : with, once more, Eqs. ( 29) and (30) give, for k = 1, We have seen previously that this latter expression may not be convergent.
Let Z = ln (Y ).We can easily calculate Eq. ( 33) means that ln (ξ ) − B is an unbiased estimator of ln (θ 0 ).Hence, the log-unbiased estimator of θ 0 appears to be, in the model world, ξ ′, defined as In the measurement world, we have, from (29 (37) To conclude this section, a log-estimator can easily be rendered unbiased by multiplying the linear estimator by the right constant, with the important further advantage that, in the measurement world, the unknown true value has, for its mean, the actual measurement data.This equivalence between the model and measurement worlds for a scale parameter was already noticed 50 years ago [9].

E. Generalization to a χ ν 2 Distribution
The results of Section II-C can be extended to any χ ν 2 distribution, regardless of the number of degrees of freedom ν.
1) Model World: The reduced variable is .
2) Measurement World: The reduced variable is and The mean (ν > 2) is The logarithmic mean is Table II and Fig. 4 show the biases and the 95% confidence interval versus EDF.
and, with ν = 1, The mathematical expectation of Lθ is then Lθ is then an unbiased estimator of ln(θ 0 ).
On the other hand, from (45), for measurement data ξ 0 : and from (24), with ν = 1, The mathematical expectation of ln(θ ) is equal to ˆ. L θ Therefore, e Lθ is a log-unbiased estimator of θ.For fitting a log-log plot, this log-unbiased estimator is then optimal.

B. Fitting an Allan Variance Curve
Fig. 5 displays the bias estimates (crosses) as well as the log-unbiased estimates (circles) of a simulated oscillator.Fitting the AVAR estimates rather than the logunbiased estimates leads to a significant underestimation of the long-term asymptotes (τ and τ2 asymptotes) and therefore of the very-long-term stability of the oscillator.

C. Simulation Results
1) Pure Random-Walk FM: We first simulated pure random-walk FM, because it is often the dominating noise for the very-long-term.a) Principle of the simulation: We performed one million simulations of a sequence of 1024 frequency deviation samples of pure random-walk FM with a noise level (model parameter) within the interval: 10 −3 < h −2 < 10 +9 with a uniform distribution of the exponent 2 between −3 and +9.Thus, this variation range can cover all of the model parameters which may give a measurement data ĥ−2 equal to 1.It is then possible to select either the realizations for which the model parameter h −2 is close to 1 at 5% (i.e., within the interval [0.952, 1.05]), or the realizations for which the measurement data ĥ−2 is close to 1 at 5%. b) Study of AVAR for the higher τ value: It is well known that an AVAR measurement data for the higher τ value is χ 1 2 distributed.Because each sequence contains 1024 samples, the higher τ argument of AVAR is τ max = 1024/2 = 512 arbitrary units (a.u.).From Table I, the theoretical AVAR result (model parameter) for h −2 = 1 a.u. and τ = 512 a.u. is Let us now consider the theoretical σ y   surement data.To more easily compare these values, we normalize them by the coefficient K n given by (51): where Y 1 and Y 2 are, respectively, the average of the first 512 frequency deviation samples and the average of the last 512 frequency deviation samples of a considered simulated sequence.Thus, we can select either the sequences for which θ = θ 0 = 1 ± 5% (model world) or the sequences for which ξ = ξ 0 = 1 ± 5% (measurement world).We can then easily compare the histogram of the ξ measurement data, knowing that θ 0 = 1 ± 5%, to the theoretical conditional probability given by (18) as well as the histogram of the θ parameter, knowing that ξ 0 = 1 ± 5%, to the theoretical conditional probability given by ( 22).Fig. 6 shows the very good agreement between histograms and theoretical conditional probability for both the model world and measurement world.
Again, the results of this simulation are in perfect agreement with the theoretical considerations we have developed previously.The estimation of the random-walk noise level is clearly better with the fit over the log-unbiased estimates, but the estimation of h −2 over the classical estimates is only slightly underestimated because it is based on all estimates, including those obtained for low τ values that are little biased (see Fig. 7).
2) White FM + Random Walk FM: We simulated 10 000 sequences of a combination of white FM and random-walk FM in such a way that the random-walk FM is the dominating noise over the last decade.For the simulations, we used model parameter: h 0 = 1 a.u., h −2 = 1.5 • 10 −5 a.u.; fit over classical estimates: 〈 〉 ĥ0 = 0.969 a.u.(−3.1%), 〈 〉 − ĥ 2 = 1.01 • 10 −5 a.u.(−33%); and fit over log-unbiased esti-  In this case, the estimation of the random-walk FM noise level is drastically improved by the use of the logunbiased estimates because the random walk asymptote is fitted over the very last τ values, which are significantly affected by the log-log bias (see Fig. 8).For comparison, Theo1 estimates [14] were added on Fig. 8.
IV. C AVAR estimates are not biased and are easily fitted by weighted relative least squares.However, ln(AVAR) estimates are biased, and this log-log bias affects the fits performed on log-log plots.
A log-unbiased estimator has been defined in such a way that its log is an unbiased estimator of the log of the model parameter.On a log-log AVAR plot, it is then recommended to fit the log-unbiased estimates rather than the classical AVAR estimates.
For this purpose, we developed the SigmaTheta software package [15] which performs fits according to the criterion (4) including log-unbiased estimator.It is free software developed under a GNU license for UNIX and Windows (Microsoft Corp., Redmond, WA) systems, available at http://theta.obs-besancon.fr/spip.php?article103.R Fig. 8. Example of one realization among 10 000 simulations of an Allan deviation (ADEV) curve for a combination of white FM and randomwalk FM (obtained using SigmaTheta-1.2[15]).As in Fig. 1, the graph shows ADEV but the fit was performed over the Allan variance (AVAR) estimates.The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.For comparison, Theo1 estimates [14] are represented as asterisks.

Fig. 2 .
Fig. 2. Dispersion of Allan variance (AVAR) estimates for a pure white frequency noise and their equivalent degrees of freedom, ν (65 536 samples at 1 s).
back into the model world for k = 1 and using (8): , this estimator ensures that its mean is equal to the actual measurement data:

2 ( 512 )
as the model parameter and the σy 2 (512) computed value as the mea-

Fig. 4 .
Fig. 4. Biases and 95% confidence intervals over the model parameter θ versus equivalent degrees of freedom, ν, in the measurement world where the measurement data are ξ 0 = 1.The circles represent the exponential of the mathematical expectation of the log of the model parameter and the crosses represent the mathematical expectation of the model parameter (which does not converge for ν ≤ 2).

Fig. 5 .
Fig.5.Example of log-log fit (obtained with SigmaTheta-1.2[15]).As in Fig.1, the graph shows the Allan deviation (ADEV) but the fit was performed over the Allan variance (AVAR) estimates.The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.

Fig. 6 .
Fig. 6.Comparison of the histogram of the simulations and the theoretical conditional probability distribution function (PDF) (a) in the model world, for which the model parameter is θ 0 = 1, and (b) in the measurement world, for which the measurement data are ξ 0 = 1.

Fig. 7 .
Fig.7.Example of one realization among one million simulations of an Allan deviation (ADEV) curve for a pure random-walk FM with a unit noise level (obtained using SigmaTheta-1.2[15]).The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.

TABLE I .
T R   A V   D N T   P L M    L F D.

TABLE II .
B  95% C I B V EDF ν.
Let us define a new estimator Lθ as ˆln exp,