1. Introduction

Cesium beam atomic clocks have an accuracy of a few parts in 101 over a period of a day; however, they are not deterministic and undergo stochastic variations in both time and frequency. It was shown by Tryon and Jones [1] that the actual frequency of a clock behaves as a random walk which, over a time period of a day, has a standard deviation of less than four nanoseconds per day. The effect of this random walk in frequency on time measurements is that the individual frequency steps are summed, producing an integrated random walk.

In addition, frequency has a white noise component which integrates into an independent random walk in time with a standard deviation of up to 15 nanoseconds over a time period of a day. Therefore, the deviation of clock time from true time behaves as a random walk plus the sum of a random walk. The possibility of a frequency drift also exists.

When estimating time from an ensemble of clocks, the only observations possible are measurements of time differences between clocks. In addition, observations may occur at unequally spaced time points, and various types of errors are possible. The most common errors are read errors, where a single measurement is incorrect but subsequent measurements are the same as if the read error did

*Division of Biometrics, Box B-119, School of Medicine, University of Colorado, Denver, CO 80262.

**Deceased. Dr. Tryon served with the Center for Applied Mathematics, National Engineering Laboratory.

not occur; time steps, where a jump in time occurs and the clock remains at the new position; and frequency steps, where the frequency takes a step and remains at the new position. This last error appears as a change in the rate of gain of the clock.

This paper describes a state space algorithm for estimating time from an ensemble of cesium beam atomic clocks with unequally spaced observations subject to various errors. Clocks can be added to or deleted from the ensemble. Using iterative calculations over a recent history of the ensemble, maximum likelihood estimates of the unknown variances can be obtained by nonlinear optimization. Confidence intervals on the estimated parameters and test of hypotheses, such as whether parameters are significantly different from zero, can also be obtained from the likelihood function.

is expressed in units of nanoseconds per day. w(t) represents a possible frequency drift in units of nanoseconds/day2, and d(t) is the time interval between t-1 and t in days. ε(t), n(t) and a(t) are random variables with zero mean, uncorrelated with each other, and uncorrelated in time. These are the input to the random walks, and for small d(t), their variances are proportional to the length of the time intervals between observations, ó(t),

If oa2

20, the drift, w(t), does not change with time. The state equation for an ensemble of m clocks is obtained by concatenating the three state elements of each clock into a column vector of length 3m. The state transition matrix, (t), is a 3m by 3m matrix consisting of 3 by 3 blocks on the main diagonal corresponding to the state transition matrix in eq (1), and zero blocks off the main diagonal. This equation of state for m clocks will be written.

X(t) = (t) X(t−1) + U(t).

2

2

(2.3)

2

The random input vector, U(t), is of length 3m and consists of the ɛ(t), n(t), and a(t) for each clock. The covariance matrix of U(t), d(t) Q, is a 3m by 3m diagonal matrix with diagonal elements consisting of the on and o2 for each clock. While reasonable guesses are available for the σ and on, it is necessary to obtain better estimates of these variances from data. The values are characteristics of each clock and will be different for each clock. It is also necessary to determine whether the o are significantly different from zero. If the o2 are not significantly different from zero it would indicate that any drifts that exist are deterministic rather than random walks. In this case the w(t) will be constant with respect to time and these can be tested to determine if the drifts are significantly different from zero.

2

2

and has dimension (m-1) by 3m. H(t) can be time dependent since if observations are missing or deleted because of error detection, rows of H(t) are eliminated. V(t) is a vector of observational errors. These errors are very small when reading clock differences. If the data are truncated to the nearest nanosecond, the variances of the elements V(t) would be

In a companion paper [1] Tryon and Jones used 333 daily observations from seven clocks starting February 16. 1979 to estimate clock parameters. These data were chosen since they consisted of a fairly long record with few errors. Some preprocessing was necessary to remove several outliers giving a set of data that could be considered error-free. These calculations were carried out on a CDC Cyber 170/750 computer.

This paper reports the new statistical procedures and techniques that were incorporated when the algorithm was rewritten to run on the National Bureau of Standards Time and Frequency Division's PDP-11/70. The algorithm was initialized on March 31, 1981 with the first set of data collected on April 1, 1981. At the beginning of each month the daily observations from the previous month are passed to the algorithm without preprocessing. While the data are nominally collected at the same time each day, for various reasons, the time of the observations are sometimes off by several hours giving unequally spaced data. Individual clocks may have read errors, time steps, and frequency steps. In addition, clocks may be added or deleted from the ensemble.