Stationary processes are an important tool in many practical applications of time series analysis, and the topic is extensively studied in the literature. Traditionally, stationary processes are modelled by using autoregressive moving average processes or linear processes (see monographs [2, 4] for details).

One of the most simple example of an autoregressive moving average process is an autoregressive process of order one. That is, a process (Xt)t∈Z defined by (1) Xt=ϕXt−1+εt,t∈Z, where ϕ∈(−1,1) and (εt)t∈Z is a sequence of independent and identically distributed square integrable random variables. The continuous time analogue of (1) is called the Ornstein–Uhlenbeck process, which can be defined as the stationary solution of the Langevin-type stochastic differential equation (2) dUt=−ϕUtdt+dWt, where ϕ>0 }0$]]> and (Wt)t∈R is a two-sided Brownian motion. Such equations have also applications in mathematical physics.

Statistical inference for AR(1)-process or Ornstein–Uhlenbeck process is well-established in the literature. Furthermore, recently a generalised continuous time Langevin equation, where the Brownian motion W in (2) is replaced with a more general driving force G, have been a subject of an active study. Especially, the so-called fractional Ornstein–Uhlenbeck processes introduced by [3] have been studied extensively. For the parameter estimation in such models, we mention a recent monograph [6] dedicated to the subject, and the references therein.

When the model becomes more complicated, the number of parameters increases and the estimation may become a challenging task. For example, it may happen that standard maximum likelihood estimators cannot be expressed in closed form [2]. Even worse, it may happen that classical estimators such as maximum likelihood or least squares estimators are biased and not consistent (cf. [1] for discussions on the generalised ARCH-model with liquidity given by fractional Brownian motion). One way to tackle such problems is to consider one parameter model, and to replace the white noise in (1) with some other stationary noise. It was proved in [8] that each discrete time strictly stationary process can be characterised by (3) Xt=ϕXt−1+Zt, where ϕ∈(0,1) . This representation can be viewed as a discrete time analogue of the fact that Langevin-type equation characterises strictly stationary processes in continuous time [7].

The authors in [8] applied characterisation (3) to model fitting and parameter estimation. The presented estimation procedure is straightforward to apply with the exception of certain special cases. The purpose of this paper is to provide a comprehensive analysis of these special cases. In particular, we show that such cases do not provide very useful models. This highlights the wide applicability of characterisation (3) and the corresponding estimation procedure.

The rest of the paper is organised as follows. In Section 2 we briefly discuss the motivating estimation procedure of [8]. We also present and discuss our main results together with some illustrative figures. All the proofs and technical lemmas are postponed to Section 3. Section 4 provides a small simulation study comparing an estimator of quadratic type arising out of (3) with the classical Yule–Walker estimator in the case of an AR (1) -process. We end the paper with discussion.

Let X=(Xt)t∈Z be a stationary process. It was shown in [8] that equation (4) Xt=ϕXt−1+Zt, where ϕ∈(0,1) and Zt is another stationary process, characterises all discrete time (strictly) stationary processes. Throughout this paper we suppose that X and Z are square integrable processes with autocovariance functions γ(·) and r(·) , respectively. Using Equation (4), one can derive the quadratic equations of the Yule–Walker type for the parameter ϕ, which can be solved in an explicit form. Namely, for any m∈Z such that γ(m)≠0 we have (5) ϕ=γ(m+1)+γ(m−1)±(γ(m+1)+γ(m−1))2−4γ(m)(γ(m)−r(m))2γ(m). The estimation of the parameter ϕ is obvious from (5) provided that one can determine which sign, plus or minus, should be chosen. In practice, this can be done by choosing different lags m for which the covariance function γ(m) is estimated. Then one can determine the correct value ϕ by comparing different signs in (5) for different lags m (We refer to [8, p. 387] for detailed discussion). However, this approach fails, i.e. one cannot find suitable lags leading to the correct choice of the sign and only one value ϕ, if, for m∈Z such that γ(m)=0 we also have r(m)=0 , and for any m∈Z such that γ(m)≠0 , the ratio (6) am=r(m)γ(m)=a for some constant a∈(0,1) . The latter is equivalent [8, p. 387] to the fact that γ(m+1)+γ(m−1)γ(m)=b for some constant b with ϕ<b<ϕ+ϕ−1 . This leads to (7) γ(m+1)=bγ(m)−γ(m−1). Moreover, if γ(m)=r(m)=0 for some m, it is straightforward to verify that (7) holds in this case as well. Thus (7) holds for all m∈Z . Since covariance functions are necessarily symmetric, we obtain an “initial” condition γ(1)=b2γ(0) . Thus (7) admits a unique symmetric solution.

By the Cauchy–Schwarz inequality and equality γ(1)=b2γ(0) , it is clear that (7) does not define covariance function for b>2 }2$]]>. Furthermore, since ϕ>0 }0$]]>, we conduct a comprehensive analysis of the special cases by studying the functions given by (7) with b∈[0,2] (we include the trivial case b=0 ). For b=2 Equation (7) corresponds to the case Xt=X0 for all t∈Z which is hardly interesting. Similarly, the case b=0 leads to a process (…,X0,X1,−X0,−X1,X0,X1,…) which again does not provide a practical model. On the other hand, it is not clear whether for some other values b∈(0,2) Equation (7) can lead to some non-trivial model in which the estimation procedure explained above cannot be applied. By our first main theorem, for any b∈[0,2] , Equation (7) defines a covariance function. On the other hand, the resulting covariance function, denoted by γb , leads to a model that is not very interesting either.

In many applications of stationary processes, it is assumed that the covariance function γ(·) vanishes at infinity, or that γ(·) is periodic. Note that the latter case corresponds simply to the analysis of finite-dimensional random vectors with identically distributed components. Indeed, γ(m)=γ(0) implies Xn=X0 almost surely, so periodicity of γ(·) with period N implies that there exists at most N random variables as the source of randomness. By items (2) and (3) of Theorem 2.1, we observe that, for suitable values of b, (7) can be used to construct covariance functions that are neither periodic nor vanishing at infinity. On the other hand, in this case there are arbitrary large lags m such that γb(m) is arbitrary close to γb(0) . Consequently, due to the strong dependency structure, it is expected that different estimation procedures will fail. Indeed, even the standard covariance estimators are not consistent. A consequence of Theorem 2.1 is that only a little structure in the noise Z is needed in order to apply the estimation procedure of the parameter ϕ introduced in [8], provided that one has consistent estimators for the covariances of X. The following is a precise mathematical formulation of this observation.

In most situations, a natural assumption regarding the covariance of the noise is r(m)→0 as m→∞ . In this case, Theorem 2.2 gets obviously satisfied. We end this section by visual illustrations of the covariance functions defined by (7). In these examples we have set γb(0)=1 . In Figures 1a and 1b we have illustrated the case of item (1) of Theorem 2.1. Note that in Figure 1a we have b=2sin(13π2)=1 . Figures 1c and 1d demonstrate how k can affect the shape of the covariance function. Finally, Figures 1e and 1f illustrate the case of item (2) of Theorem 2.1.

Throughout this section, without loss of generality, we assume γb(0)=1 . We also drop the subscript and simply write γ for γb . The following first result gives explicit formula for the solution to (7).

Before proving our main theorems we need several technical lemmas.

We are now ready to prove Theorem 2.1 and Theorem 2.2. Proof the Theorem 2.1.

Throughout the proof we denote a2≡a1(mod2π) if a2=a1+2kπ for some k∈Z . That is, a1 and a2 are identifiable when regarding them as points on the unit circle. By a3∈(a1,a2)(mod2π) we mean that a3≡a(mod2π) for some a∈(a1,a2) . 1.

Since arcsin(b2)=klπ2 , the first claim follows from Proposition 3.1 together with the fact that functions sin(·) and cos(·) are periodic. In particular, we have γ(4l+m)=γ(m) for every m∈Z .

In this section we present a simulation study in order to compare the classical Yule–Walker estimator with our quadratic type estimator in the case of an AR (1) -process.

If (Zt)t∈Z is chosen to be white noise in the characterisation (4), the process is an AR (1) -process with ϕ>0 }0$]]> and equations (5) provide natural estimators for the AR (1) -parameter. In this case, it can be verified that the minus sign in (5) is the correct choice whenever m≠0 (see the discussion about determining the correct sign based on the ratio r(m)γ(m) in [8]). If m=0 , then the discriminant of (5) equals to zero yielding (16) ϕ=γ(1)γ(0), which is the classical Yule–Walker equation for the model parameter. The same equation is also given by Theorem 2 of [8]. We would also like to point out that, when k=0 is chosen, the other Yule–Walker equation γ(0)=ϕγ(1)+r(0) related to AR(1)-processes is given by Equation 7 of [8].

Figure 3 displays histograms comparing efficiencies of the Yule–Walker estimator given by (16) and the alternative estimator given by (5) with m=1 . We simulated data from an AR (1) -process with ϕ=0.5 . The used sample size and number of iterations were 10000 and 1000, respectively. The sample mean and sample variance of the alternative estimates (3a) are 0.5000318 and 0.0001030415. For classical Yule–Walker estimates (3b), the corresponding sample statistics are 0.4998632 and 7.771528e-05. Thereby it seems that, in this setting, the variances of the two estimators are of the same order. Moreover, the slightly better performance of the Yule–Walker estimator is something that could have been expected. Indeed, the autocovariance function of an AR (1) -process is exponentially decreasing and consequently, the denominator γ(m) acts as an variance increasing factor in the estimators as m grows. For asymptotic distributions of estimators given by (5) and a more extensive simulation study, we refer to [8].

We have shown (Theorem 2.2) that the estimation method arising out of the AR (1) -characterisation introduced in [8] is applicable except in a very special class of processes. These processes are highly degenerate, namely, they are either driven by two random variables only (Remark 3.7) or their covariance functions can be approximated with covariance functions of such processes (proof of Theorem 2.1 item 3).

The discussed estimation procedure has recently been applied in practice in [5], where we considered a generalized ARCH model with stationary liquidity. As mentioned in [1], the usual approaches of LS and ML methods fail even in the case of liquidity given by the squared increments of a fractional Brownian motion with H≠12 . However, with our method, we were able to derive estimators for the model parameters in the case of liquidity given by a general class of stationary processes. In a more general context, it could be argued that whenever it is possible to derive the maximization problem related, e.g., to ML and QL methods in an adequate way, then these methods provide more efficient estimators. However, deriving ML or QL estimators may turn out to be a difficult task to accomplish. Moreover, unlike our method, these approaches require that the practitioner knows the underlying distribution up to the parameters of interest.

In addition to the above mentioned generalized ARCH model, it would be interesting to study whether our method can be applied in modeling and estimation of different temporal (stationary) models. For example, one could consider GARCH-X models or even integer valued processes such as INAR (1) . However, caution should be taken when comparing different methods and interpreting the results since, in general, the parameters arising out of the AR (1) -characterisation do not coincide with the parameters of the original model.