Estimation of a stationary covariance model

Let X: \Omega \times \cD \rightarrow \Rset^d be a multivariate stationary normal process of dimension d. We only treat here the case where the domain is of dimension 1: \cD \in \Rset (n=1). If the process is continuous, then \cD=\Rset. In the discrete case, \cD is a lattice. X is supposed a second order process with zero mean. It is entirely defined by its covariance function C^{stat}: \cD \rightarrow \mathcal{M}_{d \times d}(\Rset), defined by C^{stat}(\tau)=\Expect{X_sX_{s+\tau}^t} for all s\in \cD. In addition, we suppose that its spectral density function S : \Rset \rightarrow \mathcal{H}^+(d) is defined, where \mathcal{H}^+(d) \in \mathcal{M}_d(\Cset) is the set of d-dimensional positive definite hermitian matrices. The objective is to estimate C^{stat} from a field or a sample of fields from the process X, using first the estimation of the spectral density function and then mapping S into C^{stat} using the inversion relation (9), when it is possible. As the mesh is a time grid (n=1), the fields can be interpreted as time series. The estimation algorithm is outlined hereafter. Let (t_{0},t_{1},\hdots t_{N-1}) be the time grid on which the process is observed and let also (\vect{X}^0, \dots, , \vect{X}^{M-1}) be M independent realizations of X or M segments of one realization of the process. Using (9), the covariance function writes:

(1)C_{i,j}^{stat}(\tau) = \int_{\mathbb{R}}\exp\left\{ 2i\pi f \tau \right\} S_{i,j}(f)\, df

where C_{i,j}^{stat} is the element (i,j) of the matrix C^{stat}(\tau) and S_{i,j}(f) the one of S(f). The integral (1) is approximated by its evaluation on the finite domain \Omega \in \Rset:

(2)C_{i,j}^{stat}(\tau) = \int_{\Omega}\exp\left\{ 2i\pi f \tau \right\} S_{i,j}(f)\, df

Let us consider the partition of the domain as follows:

  • \Omega =[-\Omega_c, \Omega_c] is subdivided into M segments \Omega = \cup_{k=1}^{M} \mathcal{M}_k with \mathcal{M}_k=[f_k - \frac{\Delta f}{2}, f_k + \frac{\Delta f}{2}]

  • \Delta f be the frequency step, \Delta f = \frac{2 \Omega_c}{M}

  • f_k be the frequencies on which the spectral density is computed, f_k = -\Omega_c + \left(k - \frac{1}{2} \right) \Delta f = \left( 2 k - 1 - M \right) \frac{\Delta f}{2} with k=1,2,\hdots,M

The equation (2) can be rewritten as:

C_{i,j}^{stat}(\tau) = \sum_{k=1}^{M}\int_{\mathcal{M}_k}\exp\left\{ 2i\pi f \tau \right\} S_{i,j}(f)\, df

We focus on the integral on each subdomain \mathcal{M}_k. Using numerical approximation, we have:

\int_{\mathcal{M}_k}\exp\left\{ 2i\pi f \tau \right\} S_{i,j}(f)\, df \approx \Delta f S_{i,j}(f_k) \exp\left\{ 2i\pi f_k \tau \right\}

\tau must match with frequency values with respect to the Shannon criteria. Thus the temporal domain of estimation is the following:

  • \Delta t is the time step, \Delta t = \frac{1}{2 \Omega_c} such as \Delta f \Delta t = \frac{1}{M}

  • \tilde{\mathcal{T}} =[-T, T] is subdivided into M segments \tilde{{\mathcal{T}}} = \cup_{m=1}^{M} \mathcal{T}_m with \mathcal{T}_m=[t_m - \frac{\Delta t}{2}, t_m + \frac{\Delta t}{2}]

  • t_m be the time values on which the covariance is estimated, t_m = -\frac{M}{2 \Omega_c} + \left(m - \frac{1}{2} \right) \Delta t = \left(2 m - 1 - M \right) \frac{\Delta t}{2}

The estimate of the covariance value at time value \tau_{m} depends on the quantities of form:

(3)\int_{\mathcal{M}_k}\exp\left\{ 2i\pi f \tau_{m} \right\} S_{i,j}(f)\, df \approx \Delta f S_{i,j}(f_k) \exp\left\{ 2i\pi f_k \tau_{m} \right\}

We develop the expression of f_k and \tau_{m} and we get:

\left\{ \begin{array}{l} 2m - 1 - M = 2 (m-1) - (M-1) \\ 2k - 1 - M = 2 (k-1) - (M-1) \end{array} \right. Thus,

(2m - 1 - M) (2k - 1 - M) = 4 (m-1)(k-1) - (M-1)(2m -1 -M) - 2 (k-1)(M-1)

and:

(2m - 1 - M) (2k - 1 - M)\frac{\Delta t}{2}\frac{\Delta f}{2} = \frac{(m-1)(k-1)}{M} - \frac{(M-1)(2m -1 -M)}{4M} - \frac{(k-1)(M-1)}{2M}

We denote :

\left\{ \begin{array}{l} \delta(m) = \exp\left\{-i \frac{\pi}{2M} (M-1)(2m -1 -M) \right\}\\ \phi_k = \exp\left\{-i \frac{\pi}{M} (k-1)(M-1) \right\} S_{i,j}(f_k) \end{array} \right.

Finally, we get the following expression for integral in (3):

\int_{\mathcal{M}_k}\exp\left\{ 2i\pi f \tau_{m} \right\} S_{i,j}(f)\, df \approx \Delta f \exp\left\{2i \frac{\pi}{M} (m-1)(k-1) \right\} \delta(m) \phi_k

It follows that:

(4)C_{i,j}^{stat}(\tau_m) = \Delta f \delta(m) \sum_{k=1}^{M} \phi_k * \exp\left\{2i \frac{\pi}{M} (m-1)(k-1) \right\}

API:

Examples: