Modeling and estimation of multivariate discrete and continuous time stationary processes

In this paper, we give a AR$(1)$ type of characterization covering all multivariate strictly stationary processes indexed by the set of integers. Consequently, we derive continuous time algebraic Riccati equations for the parameter matrix of the characterization providing us with a natural way to define the corresponding estimator under the assumption of square integrability. In addition, we show that the estimator inherits consistency from autocovariances of the stationary process and furthermore, the limiting distribution is given by a linear function of the limiting distribution of the autocovariances. We also present the corresponding existing results of the continuous time setting paralleling them to the discrete case.


Introduction
Stationary stochastic processes provide a significant instrument for modeling numerous temporal phenomena related to different fields of science. In particular, due to the evidence of long dependence structures in the real financial data, stationary processes possessing long-memory have been widely applied in mathematical finance.
When discrete time is considered, stationary data is typically modeled by applying ARMA processes or their extensions. One focal reason for popularity of ARMA processes is that for every stationary process with a vanishing autocovariance γ(·) and for every n ∈ N there exists an ARMA process X such that γ X (k) = γ(k) for |k| ≤ n. For a comprehensive overview of ARMA processes we mention [9], [16] and [31]. The immense ARMA family include for example SARIMA models, where a seasonal ARMA process is obtained by differencing the original data, and different GARCH models originating from [14] and [7] that are commonly used in financial modeling taking account of the time-dependent volatility. ARMA processes, their extensions and estimation have been concerned e.g. in [18], [36], [30], [29], [27], [15], [8], [17], [3], [26] and [43], to mention but a few. Moreover, in [38] we showed that all univariate strictly stationary processes indexed by the integers are characterized by the AR(1) equation where the noise Z belonging to a certain class of stationary process is not necessarily white. Established on the characterization, we proposed an estimation method for φ in the case of a square integrable stationary process that has several advantages over the conventional methods such as maximum likelihood and least squares fitting of ARMA models. Furthermore, in [40] we applied our method in estimation of a generalized ARCH model involving a covariate process that can be interpreted as the liquidity of an asset. In the case of continuous time, the Ornstein-Uhlenbeck process X given by the Langevin equation where θ > 0 and B is a two-sided Brownian motion, can be seen as the analogue of the discrete time AR(1) process. By posing a suitable initial condition, (1) yields a stationary solution. The foregoing can be generalized, for example, by replacing Brownian motion with other stationary increment processes satisfying certain integrability conditions, such as fractional Brownian motion recovering the fractional Ornstein-Uhlenbeck process introduced in [12]. This kind of generalized Ornstein-Uhlenbeck processes are applied e.g. in mathematical finance to describe mean-reverting systems under the influence of shocks, and they are a highly active topic of research. Equations of type (1) with varying driving forces, and estimation in such models have been concerned e.g. in [19], [21], [2], [4], [5], [10], [13], [20], [32], [33], [35], [1] and [28], to mention but a few. Furthermore, in [41] we showed that a generalized multidimensional version of (1) characterizes all multivariate strictly stationary processes with continuous paths. Consequently, we proposed an estimation method for the parameter matrix of (1) under the assumption of square integrability. The method is based on continuous time algebraic Riccati equations (CAREs) written in terms of the autocovariance function of the stationary solution. Algebraic Riccati equations have been studied intensively in the literature and they occur naturally e.g. in optimal control and filtering theory. Real-valued CAREs often take the symmetric form where C and D are symmetric, and symmetric solutions A are to be found. For a general approach to algebraic Riccati equations the reader may consult for example [24]. The existence and uniqueness of a solution to (2) is a well-studied topic, especially when C and D are positive semidefinite (see e.g. [22], [42] or [34]).
The rest of the paper is organized as follows. In Subsection 2.1 we complete our previous investigations of stationary processes by treating the multivariate discrete time case. First, we show that the characterization is now given by a multidimensional AR(1) type of equation. Then, by taking a similar approach as in [41] we obtain a set of symmetric CAREs that serve as a basis for estimation of the model parameter matrix. Finally, we state theorems for consistency and asymptotic distribution of the estimator. In Subsection 2.2 we present the main results of [41] while at the same time comparing them to the results obtained in discrete time. The proofs are postponed to Section 3.

Main results
The considered processes are n-dimensional, real-valued and indexed by I ∈ {Z, R}.
t . Equality of the distributions of two random vectors Y t and Z t is denoted by Y t law = Z t . Similarly, equality of two processes Y and Z in the sense of finite dimensional distributions is denoted by Y = (Y t ) t∈I law = (Z t ) t∈I = Z. Throughout the paper, we investigate strictly stationary processes meaning that (X t+s ) t∈I law = (X t ) t∈I for every s ∈ I. Consequently, we omit the word 'strictly' and simply say that X is stationary. By writing A ≥ 0 or A > 0 we mean that the matrix A is positive semidefinite or positive definite, respectively. We denote an eigendecomposition of a symmetric matrix by A = QΛQ ⊤ , where Λ = diag(λ i ). Furthermore, the L 2 vector norm and the corresponding induced matrix norm is denoted by · .
By applying the models of stationary processes, which we introduce in this paper, we consequently obtain symmetric CAREs of the form where C, D ≥ 0, and we are solving the equation for a positive definite A. There exists a vast amount of literature on existence and uniqueness of a solution (see e.g. [22] or [42]) in the described setting. In particular, if C, D > 0, then there exists a unique positive semidefinite solution to (3). Furthermore, there exists several numerical methods for finding the positive semidefinite solution of (3) (see e.g. [11], [25] or monograph [6]).

Discrete time
In this subsection we extend the characterization of stationary processes of [38] to multivariate settings. Consequently, we derive quadratic equations for the corresponding model parameter matrix providing us with a natural way to define an estimator for the parameter. Finally, we pose theorems for consistency and asymptotic distribution. A strong analogue with the continuous time case I = R covered in [41] is obtained. We start by providing some definitions.
Definition 2.1. Let G = (G t ) t∈Z be a n-dimensional stationary increment process. We define a stationary process ∆G = (∆ t G) t∈Z by As in the univariate case, we define a class of stationary increment process having sub-exponentially deviating sample paths.
Definition 2.2. Let H > 0 and let G = (G t ) t∈Z be a n-dimensional stochastic process with stationary increments and G 0 = 0. If lim l→−∞ 0 k=l e kH ∆ k G exists in probability and defines an almost surely finite random variable, we denote G ∈ G H . Remark 2.3. Lemma 3.1 shows that existence of a logarithmic moment is sufficient for G ∈ G H for all H > 0. Particularly, this is the case if G is square integrable. On the other hand, an example of a stationary increment process G with G 0 = 0, but G / ∈ G H for any H > 0 was provided in [37].
The next theorem characterizes all multivariate stationary processes, including processes possessing long-memory.
Theorem 2.4. Let H > 0 be fixed and let X = (X t ) t∈Z be a n-dimensional stochastic process. Then X is stationary if and only if lim t→−∞ e tH X t P = 0 and for G ∈ G H and t ∈ Z. Moreover, the process G ∈ G H is unique.
Corollary 2.5. Let H > 0 be fixed and let X be stationary. Then X admits an AR(1) type of representation where Φ = e −H and G ∈ G H .
By using (5) and the expression (14) from the proof of Theorem 2.4, it is straightforward to show that ∆G is centred and square integrable if and only if X is centred and square integrable, respectively. In what follows, we assume these two attributes and write in this case also G is centred and square integrable, and we denote v(t) = cov(G t ) = EG t G ⊤ t . We would like to point out that centredness can be assumed without loss of generality (see Remark 2.11).
Under the assumptions, we obtain an expression for γ(t) in terms of the noise process.
Remark 2.6. The autocovariance function γ(t) is given by Furthermore, if G has independent components, we obtain The following lemma writes the quadratic equations for the model parameter Φ = e −H presented in [38] in our multivariate setting.
Lemma 2.7. Let H > 0 be fixed and let X be stationary of the form (5). Then for every t ∈ Z.
Remark 2.8. In the proof of the lemma, we have utilized the increment process ∆G of the noise similarly as in [38] yielding quadratic equations for Φ. However, for a general stationary X, (6) is a symmetric CARE only if t = 0, and even in this case, existence of a unique positive semidefinite solution is not guaranteed.
By applying a similar approach as in [41] by considering the noise G directly, we obtain a set of symmetric CAREs on which we construct an estimator for the model parameter Φ = e −H . For this, we define the following matrix coefficients.
Theorem 2.10. Let H > 0 be fixed and set Θ = I − e −H . Let X = (X t ) t∈Z be stationary of the form (4). Then the CARE is satisfied for every t ∈ N.
Remark 2.11. Equations (6) and (7) are covariance based. Consequently, they hold also when X and G in Theorem 2.4 are not centred.
Remark 2.12. By Lemma 3.6, the matrix Θ is positive definite. Since grows enough in time, D t becomes positive definite (see [41]). This is the case e.g. when the noise has independent components with growing variances.
We give a couple of examples on how some basic multivariate processes of ARMA type can be presented in the form (5), and how to derive the corresponding noise G together with its covariance function v.
Example 2.13. Let X be a n-dimensional stationary AR(1) type of process given by Example 2.14. Let X be a n-dimensional stationary ARMA(1, q) type of process given by Similarly as above, we may set Φ = φ and now ∆G equals to the MA(q) process on the right. Consequently, for t ≥ 1, where θ 0 = I.
In [38] we proposed an estimation method of one-dimensional stationary processes based on equations (6). In particular, we showed that the method is applicable except in some special class of stationary processes. In [39] we provided a comprehensive analysis of the class, and proved that it consists of highly degenerate processes. On the other hand, due to the strong dependence structure, the failure of different estimation methods is expected. Fundamentally, a stationary process X belongs to the class if there exists two values H andH such that the corresponding processes ∆G and ∆G in (5) have identical autocovariance functions. Next, we state a lemma showing that these degenerate processes have a special characteristic also under the new set of equations (7). Lemma 2.15. Let X be a one-dimensional stationary process and let H > 0 be fixed.
If the equation yields the same two solutions Φ,Φ > 0 for every t ∈ Z, then also the equation For estimation, it is desirable that (7) admits a unique positive semidefinite solution guaranteeing convergence to the correct parameter matrix. Hence, we simply assume that t is chosen in such a way that C t , D t > 0, and we omit the subindex t from Equation (7). We have justified the assumption of positive definiteness in detail in continuous time (see Subsection 2.1 and Remark 2.11 in [41]). Furthermore, we assume that v(t) is known and the stationary process X is observed up to the time T > t, and the coefficient matrices B, C and D are estimated from these observations by replacing the autocovariances γ(·) with some estimatorsγ T (·). The coefficient estimators are denoted byB T ,Ĉ T andD T , and we set Next, we define an estimatorΘ T for the matrix Θ = I −Φ = I −e −H . The proofs of the related asymptotic results allow a certain amount of flexibility in the definition. Thus, we give a definition that probably is the most convenient from the practical point of view. Consistency and the rate of convergence ofΘ T are inherited from autocovariance estimatorsγ T (·) of the observed stationary process. In addition, the limiting distribution is obtained as a linear function of the limiting distribution of the autocovariance estimators.
where Z is a (t + 1)n 2 -dimensional random vector, then: (1) LetZ be the permutation of elements of Z corresponding to the order of elements of  Define a linear mapping L 1 : R (t+1)n 2 → R 3n 2 by . . .
where 0 1 is an empty sum. Then (2) If D, C > 0 andΘ T is given by Definition 2.16, then where L 2 : R 3n 2 → R n 2 is a linear mapping expressible in terms of Θ, t and r.

Continuous time
We have collected the main results (Theorems 2.20, 2.27 and 2.32) of [41] considering continuous time stationary processes into this subsection. In addition, in order to complete the analogue between discrete and continuous time, we derive quadratic equations for the model parameter by using the increments of the noise process (Lemma 2.23). We assume that the processes have continuous paths almost surely and hence, the related stochastic integrals can be interpreted as pathwise Riemann-Stieltjes integrals. Again, we start by defining the class G H of stationary increment processes for H > 0. As in discrete time, it can be shown that existence of some logarithmic moments ensure that G ∈ G H for all H > 0. In particular, square integrability of G suffices, which is the case in our second moment based estimation method.
The next theorem is the continuous time counterpart of Theorem 2.4 showing that all stationary processes are characterized by the Langevin equation, whereas in discrete time, the characterization was given by an AR(1) type of equation.
Theorem 2.20. Let H > 0 be fixed and let X = (X t ) t∈R be a n-dimensional stochastic process. Then X is stationary if and only if for G ∈ G H and t ∈ R. Moreover, the process G ∈ G H is unique.

Corollary 2.21. From Theorem 2.20 it follows that X is the unique stationary solution
to (10).
In order to apply Theorem 2.20 in estimation, we pose the assumption This guarantees that G ∈ G H for all H > 0, and square integrability of X and G. On the other hand, if X is square integrable, then G is also. In addition and without loss of generality, we assume that the processes are centred. Again, we write γ(t) = EX t X ⊤ 0 and v(t) = EG t G ⊤ t . Now, the autocovariance function of the following stationary process is well-defined.
Definition 2.22. Let G = (G t ) t∈R be a centred square integrable stationary increment process and let δ > 0. We define a stationary process ∆ δ G = (∆ δ t G) t∈R by and the corresponding autocovariance function r δ by As in discrete time (Lemma 2.7), we obtain quadratic equations for the model parameter H in terms of r δ . Lemma 2.23. Let H > 0 be fixed and let X be of the form (11). Then  for every t ∈ R.
Remark 2.24. The advantage of the equations above is that we have to consider γ(s) only for s ∈ [t − δ, t + δ], but as in the discrete case, for a general stationary X we obtain a symmetric CARE only when t = 0. In addition, similarly as above, we could set in discrete time ∆ k t G := G t − G t−k , k ∈ N. However, this would lead to more complicated equations in Lemma 2.7.
A significant difference compared to the discrete time equations (6) occurs in the univariate case. Namely, the first order term with respect to H vanishes.

Corollary 2.25. The univariate case yields
for every t ∈ R.
One could potentially base a univariate estimation method on the above equations without the concern of existence of a unique positive solution. However, since we wish to treat also multivariate settings, we present the most central results of [41] that are obtained from Theorem 2.20 by considering the noise G directly. First, we define matrix coefficients corresponding to Definition 2.9. Consequently, we write symmetric CAREs for the parameter H that are similar to the CAREs (7) for the discrete time parameter Θ.
Theorem 2.27. Let H > 0 be fixed and let X = (X t ) t∈R be stationary of the form (11). Then the CARE is satisfied for every t ≥ 0.
Remark 2.28. As in discrete time, the equations (12) and (13) are covariance based and hence, they hold also for a non-centred stationary X.
Remark 2.29. Contrary to the discrete time equations (7), the first order term vanishes in the univariate setting as in (12).
Again, we assume that t is chosen in such a way that C t , D t > 0 ensuring the existence of a unique positive semidefinite solution. We have discussed this assumption in detail in [41]. We define an estimatorĤ T for the model parameter matrix H identically to the discrete time by replacing the autocovariances γ(·) in the matrix coefficients with their estimatorsγ T (·). The below given definition differs slightly from the definition in [41], but the same asymptotic results still apply. As in discrete time, asymptotic properties ofĤ T are inherited from the autocovariance estimators. However, due to the continuous time setting, instead of pointwise convergence, we have to consider functional form of convergence ofγ T (·). In [41] we have provided sufficient conditions in the case of Gaussian noise G with independent components, under which the assumptions of the following theorems are satisfied. In particular, the results are valid for fractional Brownian motion that is widely applied in the field of mathematical finance. (1) LetỸ s be the permutation of elements of Y s that corresponds to the order of elements of vec (γ T (s) − γ(s)) ⊤ . Then (2) If C, D > 0 andĤ T is given by Definition 2.30, then where L 2 : R 3n 2 → R n 2 is a linear mapping expressible in terms of H, t and the covariance function of G.

Proofs
In the following, we denote the smallest eigenvalue of H > 0 by λ min . Consequently e kH = Qdiag(e λ i k )Q ⊤ = e λ min k for a negative k.

Discrete time
The proof of the next lemma follows the lines of the proof of Theorem 2.2. in [37] that concerns the one-dimensional continuous time case. However, in our setting, we obtain a weaker sufficient condition for G ∈ G H for all H > 0.
Lemma 3.1. Let G = (G t ) t∈Z be a stationary increment process with G 0 = 0. Assume that Proof. Let H > 0. We apply the Borel-Cantelli lemma together with Markov's inequality to show that e kH ∆ k G → 0 almost surely as k → −∞. Let ǫ > 0 be fixed below.
since ∆G is stationary and G 0 = 0. Furthermore, k( log ǫ k − λ min ) ≥ −Ck for some C > 0 and k ≤ k ǫ . Thus, for k ≤ k ǫ , giving the wanted result. We conclude the proof by noting that We next extend the concept of self-similarity to discrete time multivariate processes.
Definition 3.2. Let H > 0 and let Y = (Y e t ) t∈Z be a n-dimensional stochastic process.
for every s ∈ Z.
The following transform and the corresponding theorem giving one-to-one correspondence between self-similar and stationary processes were originally introduced by Lamperti in the univariate continuous time setting ( [23]). Definition 3.3. Let H > 0, and let X = (X t ) t∈Z and Y = (Y e t ) t∈Z be n-dimensional stochastic processes. We define Proof. First, let X be stationary and set Z e t = (L H X) e t . Then for every m ∈ N, t ∈ Z m and s ∈ Z. Hence, Z is stationary completing the proof.
Before the proof of Theorem 2.4 we state an auxiliary lemma.
Lemma 3.5. Let H > 0 and let (Y e t ) t∈Z be a n-dimensional H-self-similar process. We define a process G = (G t ) t∈Z by Proof. It is straightforward to verify that ∆ t G = e −tH ∆ t Y e t for every t ∈ Z. In addition where by self-similarity of Y Hence, we set Proof of Theorem 2.4. Assume that lim t→−∞ e tH X t P = 0 and (4) holds for G ∈ G H . Then, by using (4) repeatedly for every n ∈ N. Since, as n → ∞, the limit of the sum above is well-defined, and lim n→∞ e (t−n−1)H X t−n−1 P = 0, we obtain that Let m ∈ N, t ∈ Z m and s ∈ Z.
Defining G as in Lemma 3.5 completes the proof of the other direction.
To prove uniqueness, we use (14). Assume that, for G,G ∈ G H , for every t ∈ Z. Then e tH X t − e (t−1)H X t−1 = e tH ∆ t G = e tH ∆ tG .
Since e tH is invertible and both processes start from zero, we conclude that G = G. Proof of Lemma 2.7. We have that Taking expectations yields Proof of Theorem 2.10. Let Then for t ∈ N we have Proof of Lemma 2.15. Assume thatΦ = e −H satisfies (8) for every t ∈ Z and set Consequentlỹ wherer(t) is the autocovariance function of (∆ tG ) t∈Z . Now, since G 0 =G 0 = 0, we obtain that for all t ∈ N. Hence, both Θ andΘ are solutions to (9).
In order to show thatΘ T is consistent, we simply need to find suitable bounds for ∆ T B, ∆ T C and ∆ T D in terms of the autocovariance estimators. After that, the same strategy as in [41] can be applied. Now, since v(t) is known, Moreover Proof of Theorem 2.17. The result follows by replacing sup s∈[0,t] γ T (s) − γ(s) with M t,T in Corollary 3.14 and in the proof of Theorem 2.9 of [41]. The details are left to the reader.
Proof of Theorem 2.18. For the first part of the theorem, we notice that where −1 0 and 0 1 are interpreted as empty sums. Now we have that

Continuous time
We only provide the proof of Lemma 2.23, while the other proofs can be found from [41].