Mathematics

M. Arato, Academician A. N. Kolmogorov, Ya. G. Sinai

On the Estimation of Parameters of a Complex Stationary Gaussian Markov Process

§ 1. We shall consider a two-dimensional stationary random process whose components \(\xi(t)\) and \(\eta(t)\) satisfy the stochastic differential equations

\[ \begin{aligned} d\xi &= -\lambda \xi\,dt-\omega \eta\,dt+d\varphi,\\ d\eta &= \omega \xi\,dt-\lambda \eta\,dt+d\psi, \end{aligned} \tag{1} \]

where \(\varphi(t)\) and \(\psi(t)\) are two independent Wiener processes with

\[ M\,d\varphi=M\,d\psi=0,\qquad M(d\varphi)^2=M(d\psi)^2=a\,dt . \]

Putting

\[ \zeta=\xi+i\eta,\qquad \chi=\varphi+i\psi,\qquad \gamma=\lambda-i\omega, \]

one can write system (1) in the form of a single equation

\[ d\zeta=-\gamma \zeta\,dt+d\chi . \tag{1a} \]

The complex correlation function of our process has the form

\[ C(\tau)=A(\tau)+iB(\tau)=M[\zeta(t)\overline{\zeta(t+\tau)}] =\sigma^2\exp(-\lambda|\tau|-i\omega\tau), \tag{2} \]

where \(\sigma^2=a/\lambda\).

If the process is observed on the interval \([0,T]\), then one can determine the empirical correlation function

\[ c(\tau)=a(\tau)+ib(\tau)=\frac{1}{T-\tau}\int_0^{T-\tau}\zeta(t)\overline{\zeta(t+\tau)}\,dt . \tag{3} \]

The empirical correlation function, with probability one, has at zero the right derivative

\[ c'(0)=-a-\frac{1}{T}s_1^2+\frac{1}{T}s_2^2-ir, \]

where \(a\) is the parameter introduced above, characterizing the intensity of the “white noises” \(\varphi'(t)\) and \(\psi'(t)\), and

\[ s_1^2=\frac{1}{2}\left[|\zeta(0)|^2+|\zeta(T)|^2\right],\qquad s_2^2=\frac{1}{T}\int_0^T|\zeta(t)|^2\,dt,\qquad r=\frac{1}{T}\int |\zeta(t)|^2\,d\theta . \]

The integration in the expression for \(r\) is carried out with respect to the angular argument \(\theta\), determined from

\[ \zeta(t)=|\zeta(t)|e^{i\theta(t)} . \]

In Fig. 1 the empirical correlation function \(c(\tau)\) is shown for Chandler variations of the coordinates of the Earth’s pole*.

* The instantaneous axis of rotation of the Earth moves relative to the minor axis of the terrestrial ellipsoid (the so-called “free nutation”). In these displacements there is a periodic component with a one-year period. After their elimination there remain Chandler displacements, which have a tendency toward oscillations with a period of the order of 14 months, but are not strictly periodic, and with large, mainly smooth (waves of the order of 10–20 years), changes in amplitude. Fig. 1 shows that the Chandler component of the displacement of the pole satisfies well the hypotheses set out at the beginning of our note.

§ 2. The parameter \(a\) is determined exactly from the realization. It remains to consider the problem of estimating the parameters \(\lambda\) and \(\omega\). Let \(P\) denote the probability measure in the space of realizations of our process on the interval \([0,T]\). In the same space introduce the standard measure

\[ V=L\times W, \]

where \(L\) is the ordinary Lebesgue measure in the plane of \(\xi(0)\), and \(W\) is the two-dimensional Wiener measure in the space of increments \(\xi(t)-\xi(0)\) with those characteristics

Fig. 1. \(\tau=0.1\cdot n;\ n=0,1,\ldots,156\)

that were adopted for the random process \(x(t)\). It can be shown that (cf. \((^{2,3})\))

\[ \frac{dP}{dV} = C\lambda \exp\left[ -\frac{\lambda^2+\omega^2}{2a}Ts_2^2 -\frac{\lambda}{a}s_1^2 +\lambda T +\frac{\omega}{a}Tr \right], \tag{4} \]

where \(C\) is some constant. Formula (4) shows that the system of three statistics \(s_1^2, s_2^2, r\) is a sufficient system of statistics for the problem. Differentiating

\[ L=\log\frac{dP}{dV} = c' + \log\lambda -\frac{\lambda^2+\omega^2}{2a}Ts_2^2 -\frac{\lambda}{a}s_1^2 +\lambda T +\frac{\omega}{a}Tr \]

with respect to \(\omega\) and \(\lambda\), we obtain the equations

\[ \frac{\partial L}{\partial \omega} = -\frac{\omega}{a}Ts_2^2 +\frac{T}{a}r =0; \tag{5} \]

\[ \frac{\partial L}{\partial \lambda} = \frac{1}{\lambda} -\frac{\lambda}{a}Ts_2^2 -\frac{s_1^2}{a} +T =0 \tag{6} \]

for determining the maximum-likelihood estimates \(\hat{\omega}\) and \(\hat{\lambda}\). From (5) we obtain

\[ \hat{\omega}=\frac{r}{s_2^2}. \tag{*} \]

Fig. 1 was obtained by processing the data of Table 6 in the work of A. Ya. Orlov \((^1)\). From the coordinates \(x(t), y(t)\) of Table 6, the component with an annual period was separated out, and the remainder was taken as \(\xi(t)\) and \(\eta(t)\). In the figure, crosses indicate the points corresponding to increments \(\tau\) of \(0.1\) year. From the figure one can immediately see that the period \(2\pi/\omega\) is approximately equal to 14 months. The correct character of the spiral obtained may lead to the supposition that the parameter \(\lambda\) also admits a very accurate estimate. This, however, is probably not so, as will be explained at the end of the note. A more detailed account of the computational technique and discussion of the results will be published elsewhere.

It can be shown that

\[ \frac{\hat{\omega}-\omega}{\sigma(\hat{\omega})},\qquad \sigma^2(\hat{\omega})=\frac{a}{T s_2^2} \]

is distributed normally as \((0,1)\) (this is an exact, not an asymptotic result). Equation (6) always has a unique positive solution. Denoting \(\lambda T=\varkappa\), \(\hat{\lambda}T=\hat{\varkappa}\), we obtain for \(\hat{\varkappa}\) the equation

\[ h_2\hat{\varkappa}^{\,2}+(h_1-1)\hat{\varkappa}-1=0, \]

where \(h_1=s_1^2/aT,\ h_2=s_2^2/aT\).

The distribution of the statistics \(h_1\) and \(h_2\), and consequently of \(\hat{\varkappa}\), depends only on the parameter \(\varkappa\). Since the distribution of \(\hat{\varkappa}\) is continuous, for any \(\alpha\), \(0<\alpha<1\), and \(\varkappa\), \(0<\varkappa<\infty\), one can find such a \(k\) that

\[ \mathbf{P}\{\hat{\varkappa}>k\mid \varkappa\}=\alpha . \tag{7} \]

Inverting the dependence

\[ k=k_\alpha(\varkappa), \]

we obtain the dependence

\[ \varkappa=\varkappa_\alpha(k) \]

(it has been established that, as \(\varkappa\) varies from \(0\) to \(\infty\), the function \(k_\alpha(\varkappa)\) varies monotonically from \(0\) to \(\infty\), so that inversion is possible and unique). Obviously,

\[ \mathbf{P}\{\varkappa<\varkappa_\alpha(\hat{\varkappa})\mid \varkappa\}\equiv \alpha . \tag{8} \]

We have arranged the computation of the function \(\varkappa_\alpha(\hat{\varkappa})\) for \(\alpha=0.1;\ 0.05;\ 0.025;\ 0.01;\ 0.005;\ 0.001;\ 0.9;\ 0.95;\ 0.975;\ 0.99;\ 0.995;\ 0.999\). The results will be published after completion of the computations.

For small \(\hat{\varkappa}\), (8) is equivalent to the relation

\[ \mathbf{P}(\hat{\varkappa}<u\varkappa)=\exp\left(-\frac{1}{u}\right), \tag{9} \]

i.e. the ratio \(\hat{\varkappa}/\varkappa\) has a \(\chi^2\) distribution with two degrees of freedom. For large \(\hat{\varkappa}\), (8) is equivalent to

\[ \mathbf{P}\left(\varkappa<\hat{\varkappa}+u\sqrt{\hat{\varkappa}}\right)\simeq \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{u} e^{-t^2/2}\,dt, \tag{10} \]

i.e. the estimate \(\hat{\varkappa}\) is asymptotically normal with variance

\[ \sigma^2(\hat{\varkappa})\sim \hat{\varkappa}. \tag{11} \]

§ 3. For the case of the motion of the Earth’s pole mentioned at the beginning of the note, on the basis of observations over \(T=60\) years, the following were obtained\(^*\)

\[ \hat{\omega}=5.274,\qquad \hat{\varkappa}=3.6,\qquad 2\pi:\hat{\omega}=1.191,\qquad \sigma(2\pi:\hat{\omega})=0.006. \]

\(^*\) Introducing Wiener \(\varphi\) and \(\psi\) into equations (1), i.e. perturbations of the “white-noise” type, is, of course, a rough idealization in the case of the motion of the Earth’s poles. It would be more correct to write

\[ \xi'=-\lambda\xi-\omega\eta+f,\qquad \eta'=\omega\xi-\lambda\eta+g. \]

However, the data from (1) show that the values of the functions \(f(t)\) and \(g(t)\) at time instants \(t\) separated from one another by several years are practically mutually independent, so that replacing the functions \(f\) and \(g\) by “equivalent white noise” is legitimate. The error in determining the intensity of this equivalent white noise is, apparently, sufficiently small not to affect substantially the results of estimating the parameter \(\lambda\). The value of \(\hat{\omega}\) was calculated by the discrete analogue of formula (*), obtained by the method of maximum likelihood for a “scheme with discrete time.”

For estimation of the parameters \(\lambda\) and \(\omega\) for the motion of the Earth’s pole, see also \((^6)\). In \((^6)\), results close to ours are given: \(\lambda=1/15,\ 2\pi:\omega=1.193\). Close values were indicated by Jeffreys \((^7)\), but in papers \((^5,^8)\) sharply different values \(\lambda=0.3\) and \(\lambda=0.01\) were indicated. The reasons for these discrepancies will be explained elsewhere.

The asymptotic formula (11) gives

\[ \sigma^2(\hat{\chi})=3.6. \]

Since \(\chi\) is known to be positive, while formula (10), for \(\alpha<0.03\), yields a negative estimate of \(k_\alpha\), it is clear that the asymptotic formula (10) is not yet applicable.

The calculation carried out by us yields the estimates

\[ \begin{gathered} \chi_{0.90}=5.5,\qquad \chi_{0.95}=6.2,\qquad \chi_{0.975}=7.8,\\ \chi_{0.10}=1.27,\qquad \chi_{0.05}=0.82,\qquad \chi_{0.025}=0.46, \end{gathered} \]

which correspond, for \(\lambda\), to the estimates

\[ \begin{gathered} \lambda_{0.90}=0.09,\qquad \lambda_{0.95}=0.10,\qquad \lambda_{0.975}=0.13,\\ \lambda_{0.10}=0.02,\qquad \lambda_{0.05}=0.01,\qquad \lambda_{0.025}=0.008. \end{gathered} \]

Moscow State University
named after M. V. Lomonosov

Received
20 II 1962

REFERENCES

¹ A. Ya. Orlov, Latitude Service, Publishing House of the Academy of Sciences of the USSR, 1958.
² Ch. T. Striebel, Ann. of Math. Stat., 30, No. 2, 559 (1959).
³ M. Arato, DAN, 145, No. 1 (1962).
⁴ E. P. Fedorov, Polar Motion, Gravimetry, Publ. 2, 3 (1948).
⁵ N. I. Panchenko, Proc. of the 14th Astronomical Conference of the USSR, Publishing House of the Academy of Sciences of the USSR, 1960, p. 232.
⁶ W. H. Munk, G. J. F. Macdonald, The Rotation of the Earth, Cambridge Monographs on Mechanics and Applied Mathematics, 1960.
⁷ H. Jeffreys, Monthly Not. Roy. Astr. Soc., Geophys. Suppl., 100, 139 (1942).
⁸ A. M. Walker, A. Young, Monthly Not. Roy. Astr. Soc., Geophys. Suppl., 115, 443 (1955); 117, 119 (1957).