This paper introduces a new statistical model for clock phases assuming a multivariate Gaussian distribution for the clock phase deviations from a common time scale. This model allows us to derive a maximum likelihood estimator for the clock phases, which is consistent with the current methods of computing a common time scale for a collection of clocks. Detailing a statistical model of the clock phases, which assumes a Gaussian distribution allows us to find the MLE for each clock’s phase deviation from a common time scale. For verification, the MLE for the clock phases is shown to be consistent with the result of the existing basic time scale equation. The statistical distribution of the frequency states resulting from this statistical model is Gaussian over a window of past time instants. This property can be used to design a new time scale based on the maximum likelihood estimator of frequency and frequency variances that are alternatives to the exponential filters designed for AT1. With the appropriate number of past frequency samples, this MLE has identical performance to the optimal AT1 algorithm in a nominal context. The statistical distribution of the frequency when the clock suffers a phase jump anomaly is then identified as a Student’s t-distribution. The Student’s t-distribution models the statistics of datasets contaminated by outliers, leading to the derivation of a different MLE that is robust to those outliers. The time scale using the robust MLE provides estimates of each clock’s frequency and frequency variance that are unaffected by phase jump anomalies and improves the long-term frequency stability when each clock in the ensemble experiences phase jump anomalies within some window of time.