A Mathematical Investigation of a Continuous Covariance Function Fitting with Discrete Covariances of an AR Process †

: In this paper, we want to ﬁnd a continuous function ﬁtting through the discrete covariance sequence generated by a stationary AR process. This function can be determined as soon as the Yule–Walker equations are found. The procedure consists of two steps. At ﬁrst the inverse zeros of the characteristic polynomial of the AR process must be ﬁxed. The second step is based on the fact that an AR process can also be seen as a difference equation. By solving this difference equation, it is possible to determine a class of functions from which a candidate for a continuous covariance function can be determined. To analyze if this function is applicable as a positive deﬁnite covariance function, it is analyzed mathematically in view of the power spectral density compared to the characteristics of the power spectral density for the discrete covariances. Then it is shown that this function is positive semi-deﬁnite. At the end, a simulation of a stationary AR(3) process is elaborated to illustrate the derived properties.


Introduction
In geodesy the observations or analyzed signals are often discrete measurements which repeat at regular distances. For example, deformation observations (repeating in time) or data from satellite missions such as GOCE (repeating in time and space). It is a common way to describe regular and equidistant signals by auto regressive moving average (ARMA) processes [1][2][3][4][5].
Within this contribution we focus on the analysis of the AR part. This AR part defines the causal link between an observation and its predecessors. Additionally, Least Squares Collocation (LSC) (see e.g., [6][7][8]) and kriging [9,10] are benefiting from the use of AR processes. For example, the inverse of a covariance matrix, based on AR process, is a band matrix witch bandwidth equals the order of the AR process (see e.g., [11]).
However, to activate the full potential of LSC a continuous covariance function is indispensable. With this function it will be possible to predict a pseudo signal between the observations. Furthermore, it is possible to predict the signal outside the observation field and not only for a multiple of the sampling rate. Moreover, refs. [12,13] used a continuous covariance function to switch from one functional to another (Like using sea level heights to calculate sea level changes which are proportional to the stream velocity).
In this paper, we want to find an analytical description of a continuous function fitting through the discrete sequence of covariances generated by any stationary AR process. This function is derived from the coefficients of the AR process and the discrete covariances using a system of equations with a unique solution. The resulting function should be positive definite, and its spectrum is expected to correspond to the spectrum of the discrete AR process. It will turn out that this function is the continuous solution of the difference equation and correctly interpolates the discrete covariance sequence with appropriate basis functions, indeed following sampling theory/convolution theorem. In the end an example is given by a simulated AR process and the accompanying continuous covariance function as well as the two spectra are estimated.

Continuous Covariance Function
To find a suitable covariance function for any stationary AR process the definition of AR processes is a good starting point. Especially the transfer from the AR process to the difference equation approach will lead us to the continuous representation we looked for (The following definitions could be found in [14][15][16][17][18]. Here the notation of [17] is used).

Construction of a Continuous Covariance Function
The process S t is called one-dimensional AR process of order p (AR(p) process) if it is described by the recursive equation where α 1 , α 2 , ..., α p are the coefficients of the AR process and E t is a i.i.d. sequence with variance σ 2 E [17] (p. 58, Equation (3.4.31)). We assume that α p = 0, as otherwise the AR(p) process is also an AR(p − 1) process so that AR(p) is not well-defined (In addition, if α p = 0 some formulas in this paper cannot be used).
An important quantity is the zeros (ζ l ) of an AR(p) process defined by the zeros of the characteristic polynomial see e.g., [17] (p. 58, Equation (3.4.32)). An alternative definition is given by the auxiliary equation if we interpret the AR process as a difference equation which has the general solution (see [19], p. 134, Equation (3.33)) These zeros p l only occur as real values or in pairs of complex conjugated zeros. Bear in mind that the zeros of the characteristic polynomial from Equation (2) are linked to the zeros of the auxiliary equation (cf. Equation (3)) by ζ l = 1/p l . AR processes are stationary if and only if the ζ l are outside the unit circle, such that |p l | < 1. In the following we will restrict to p l for simplicity.
With this definitions in mind, the discrete covariances Σ j of an AR(p) process, are linked with each other by the Yule-Walker (YW) equations (see e.g., [17], p. 59, Equation (3.4.36)), i.e., The YW equation of higher order than 0 (Equation (5)) are basically homogeneous difference equations of order p, The general solution to the difference equation can be expressed by the powers of the zeros p l of the auxiliary Equation (3). The particular solution is fixed by the boundary conditions using the discrete covariances determined from the YW equations (cf. Equation (5)), Here A l are coefficients which are complex if and only if the corresponding p l is complex. Furthermore, if there is a pair of complex conjugated p l then A l occur also as complex conjugated pairs (see e.g., [18], p. 134, Equation (3.5.44) or [19], p.163, f.).
At this point a new but now continuous function is defined, which can be seen as the continuous covariance function γ(h) : R → R for any AR(p) process, Here A l and p l are the same as in Equation (7), but the domain changed. j ∈ N 0 is replaced by h ∈ R.
Attentive readers will have noticed that γ(h) is complex if any p l ∈ R − . Then γ(h) ∈ C for all h / ∈ N 0 . One important convention that will help with this inconsistency is the use of the real part Re(γ(h)) of the complex function (see Figure 1). This condition will not have any impact if γ(h) is real (what is mostly the case), and furthermore a covariance function for real valued signals is defined to be a function in R not in C.

Properties of the Continuous Covariance Function
Since with γ(h) from Equation (8) a suitable function is found to fit through the discrete covariances from Equation (7), we want to analyze the power spectral densities of the continuous and the discrete functions. On this basis we can demonstrate that the Fourier transform of the continuous covariance function is positive semi-definite.
Initially the problem restricted to AR processes of order 1 and order 2 with two complex conjugated zeros. On the one hand any AR(p) process can be dissected into a product of AR(1) and AR (2) processes. This is a linear function, so the power spectral function is the product of the corresponding AR(1) processes and AR (2) processes. On the other hand, Equation (8) shows that the covariance function is a weighted sum of the real valued zeros, or pairs of complex conjugated zeros. So, the zeros are also in a linear relation, and so is the Fourier transform. So, it is only necessary to examine the spectrum and the Fourier transform for the first order AR process and second order AR process with two complex conjugated zeros.
For these specific types of AR processes there is an analytical solution to switch from AR coefficients α l to the zeros p l (see [20]), and for order p = 2

Power Spectral Density
The power spectral density for an AR(p) process is well known (see e.g., [16], p. 244, Equation (11.20)) and is described by the transfer function In consideration of Equations (9) and (10) the power spectral density for any AR(1) process and AR(2) process with complex conjugated zero can be calculated explicitly. So, the power spectral density AR(1) process is generated via For the AR(2) process with zeros p 1 = p * 2 the power spectral density is a little more complicated and turns out to be Using the Fourier transform of the covariance function γ(h) is an alternative way to derive the power spectral density (For further derivations of the Fourier transform see Appendix A). However, H 2 (ν) = Γ(ν). Especially Equation (11) shows H 2 (ν) as a function whose only parameter ν arise as power of the complex function e −i2πl . Therefore H 2 (ν) is a repetitive function with period 1. In contrast, Equation (14) shows Γ(ν) is aperiodic function with lim ν→∞ Γ(ν) = 0. To understand this circumstance two theorems are of importance: 1.
The discrete covariances of an AR(p) process (Σ j ) are equivalent to the product of the Dirac comb with the continuous covariance function γ(h).

2.
The convolution theorem shows that multiplication in time domain results in convolution in frequency domain.
Combining these two pieces of information shows indeed that H 2 (ν) = Γ(ν) but where ∑ ∞ k=−∞ δ(x − k) is the Dirac comb of distance 1 and Γ(ν) * ∑ ∞ k=−∞ δ(x − k) is the convolution of Γ(ν) and the Dirac comb. The transitions from continuous functions to discrete sequences as well as the resulting Fourier transforms are shown in Figure 2 (For a more detailed method of calculation for AR(1) and AR(2) processes see Appendix C).

Positive Semi-Definite Function
Equation (15) shows that knowing if H 2 (ν) is positive semi-definite is not sufficient to guarantee that the Fourier transform of the continuous function γ(h) is positive semidefinite too. Therefore, the explicit Fourier transforms of the AR(1) and AR(2) process are derived here. For the case of the AR(1) process it is easy to see that for the Fourier transform of the covariance function γ(h) holds (For the derivation of this formula see Appendix B.1). Neither the squared terms could be less than 0 nor 1 − p 2 or − ln(p) due to the fact that p lies within the unit circle. For the AR(2) process things are not that obvious. In Appendix B.2 it is demonstrated that the Fourier transform of γ(h) is To work with complex valued fractions, it is necessary to eliminate the imaginary part in the denominator. This is done by multiplying each term of the sum with the complex conjugated denominator divided by itself. Afterwards it is simple to pick the real part. To simplify the formula, we use the polar coordinates p 1 = re iφ , and p 2 = re −iφ with 0 < r < 1 and 0 ≤ φ ≤ π. So, it can be shown that the numerator is positive (Γ(ν) ≥ 0) if and only if − ln(r) coth(− ln(r)) The function f (r) and g(φ) are displayed in Figures 3 and 4 (Since g(φ) = g(−φ), the negative values of φ are not needed). On the one hand f (r) is a declining function with infimum 1. On the other hand, the function g(φ) is also declining, with a maximum of 1.

Simulation
To visualize the results from Section 2 of an AR(3) process with two complex conjugated zeros was generated as an example. First to guarantee stationarity the roots are chosen as Let the variance of the white noise be σ 2 = 1. After deriving the coefficients α l , l ∈ {1, 2, 3}, using Equation (3), we estimate the discrete covariances (Σ j ) of the AR(3) process by the reorganized YW equations (see [21], p. 32, Equation (183)). With Equation (8) a continuous function γ(h) is fitted through the discrete covariances (see the left upper corner of Figure 5). Subsequently the power spectral density is set first by Equation (11) using the coefficients α l and secondly by Equation (14). For the second step the coefficients A l are estimated by solving the system of equations Here the zeros p l and covariances Σ l are known. Each row represents Equation (7) for l ∈ {1, 2, 3}. The coefficients A l are used in Equation (8) to estimate the power spectral density of the continuous covariance function. It must be mentioned that this example is an extreme one where the Fourier transform of the continuous covariance function has high values for frequencies higher the Nyquist frequency ν n = 0.5. The Fourier transform of the continuous covariance function is not periodic at all (compare right upper corner of Figure 5). However, periodicity is the characteristic of the spectral density of a discrete AR process. Therefore, the periodicity is a result of the convolution of the Dirac comb and Γ(ν).

Conclusions and Outlook
In this paper, it was shown that the choice of a valid continuous covariance function for AR(p) processes is given by the function Here h ∈ R is the lag, p l are the roots of the characteristic polynomial, and A l follows from the unique solution of Equation (14) for p arbitrarily chosen discrete covariances Σ j 1 , Σ j 2 , ..., Σ j p (with j i ∈ N 0 , and j i = j k ⇔ i = k): Due to the convolution theorem, the power spectral density of γ(h) might be different to the power spectral density of the discrete AR(p) process. Nevertheless, the proof was given that γ(h) is still positive semi-definite (cf. Section 2.2.2), and consequently meets all conditions for a suitable covariance function. The Fourier transforms Γ(ν) and H 2 (ν) may not vary much for ν ∈ [−1, 1] and the simulation (see Section 3) is an extreme example. Anyway γ(h) is an exponential function, so it is easy to use it as functional for covariance function propagation or LSC.
In further works the continuous covariance function γ(h) could be extended to a function for an autoregressive moving average (ARMA) process to examine its properties. It is not yet demonstrated that the oscillation of Re(γ(h)) leads towards the minimal frequency if there is a real negative zero (p l = 0 for any l).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. General Fourier Transform of an AR(p) Process
In this part the Fourier transform of a function γ(h) = ∑ p l=1 A l p |h| l with h ∈ R is computed, Please note that A l and p l might be complex, but this will not have any influence on the equations.

Appendix B. Explicit Fourier Transform of the AR(1) Process and AR(2) Process with Two Complex Conjugated Zeros
In this section, the explicit Fourier transform Γ(ν) of the continuous covariance function γ(h) for the orders p = 1 and p = 2 are given as function of h, and the zeros p l .
In this context, f (x) = γ(h) (see Equation (8)) and F(ν) = Γ(ν) (see Equation (14)). For g(x) = ∑ ∞ k=−∞ δ(x − k) (the Dirac comb of distance dx = 1) the Fourier transform is again a Dirac comb of distance dν = 1/dx = 1. So, in the time domain is the same function as in the frequency domain (for ν = x: G(ν) = g(x)). Using these results leads to In step (i) the sum and the integral are exchanged.
Step (ii) represents the ability of the Dirac impulse that ∞ −∞ f (u)δ(u − x)du = f (x). Finally, (iii) uses the frequency shift of the inverse Fourier transform (see e.g., [16], p. 26, Table 2.2). Using this function to compute the power spectral density for an AR(1) process will result in Equation (12) or in the case of an AR(2) process with two complex conjugated zeros in Equation (13).