# A Flexible Class of Non-separable Cross-Covariance Functions for Multivariate Space-Time Data

Marc Bourotte<sup>a</sup>, Denis Allard<sup>a</sup> & Emilio Porcu<sup>b</sup>

<sup>a</sup>Biostatistique et Processus Spatiaux (BioSP), INRA, Avignon, France.

<sup>b</sup>Departamento de Matemática, Universidad Técnica Federico Santa María, Valparaíso, Chile

Running title: Non-separable Multivariate Space-Time Cross-Covariance Functions.

Corresponding author: D. Allard. INRA, UR546 BioSP, Site Agroparc

84914 Avignon, FRANCE

Tel: +33 432722171; Fax: +33 432722182

allard@avignon.inra.fr

## Abstract

Multivariate space-time data are increasingly available in various scientific disciplines. When analyzing these data, one of the key issues is to describe the multivariate space-time dependencies. Under the Gaussian framework, one needs to propose relevant models for multivariate space-time covariance functions, i.e. matrix-valued mappings with the additional requirement of non-negative definiteness. We propose a flexible parametric class of cross-covariance functions for multivariate space-time Gaussian random fields. Space-time components belong to the (univariate) Gneiting class of space-time covariance functions, with Matérn or Cauchy covariance functions in the spatial margins. The smoothness and scale parameters can be different for each variable. We provide sufficient conditions for positive definiteness. A simulation study shows that the parameters of this model can be efficiently estimated using weighted pairwise likelihood, which belongs to the class of composite likelihood methods. We
then illustrate the model on a French dataset of weather variables.

**Keywords** Composite likelihood; Matérn covariance; Multivariate Gaussian processes; Separability; Spatio-temporal processes; Spatio-temporal geostatistics.

# 1 Introduction

Environmental and climate sciences provide an increasing amount of multivariate data indexed by space-time coordinates. For statisticians analyzing these data, one of the key issues is to model the space-time dependence structure, not only within each variable, but also between the variables. This requires models that allow for different range and smoothness parameters for each variable, and whose parameters can be accurately estimated. Let us recall some expository material that will be needed for the presentation of our results.

Consider a $p$-dimensional multivariate random field $\mathbf{Y}(\mathbf{x}) = \{Y_1(\mathbf{x}), \dots, Y_p(\mathbf{x})\}^\top$, where $Y_i(\mathbf{x})$ represents the $i$-th variable, $i = 1, \dots, p$, and $\mathbf{x} = (\mathbf{s}, t) \in D \times T \subset \mathbb{R}^{d+1}$, $d \ge 1$, where $\mathbf{s} \in D \subset \mathbb{R}^d$ is a vector of spatial coordinates and $t \in T \subset \mathbb{R}$ is time. Let us further assume that $\mathbf{Y}(\mathbf{x})$ can be decomposed into the sum of a deterministic and a random component,

$$ \mathbf{Y}(\mathbf{x}) = \boldsymbol{\mu}(\mathbf{x}) + \mathbf{Z}(\mathbf{x}), \quad \mathbf{x} \in D \times T, $$

where $\boldsymbol{\mu}(\cdot)$ is a trend function and $\mathbf{Z}(\cdot)$ a zero mean multivariate Gaussian stationary process. Under the assumptions of Gaussianity and stationarity, the process $\mathbf{Z}(\cdot)$ is completely characterized by its matrix-valued covariance function $\mathbf{C}(\mathbf{h}, \mathbf{u}) = [C_{ij}(\mathbf{h}, \mathbf{u})]_{i,j=1}^p$, which depends solely on the space-time lag, $\mathbf{k} = (\mathbf{h}, \mathbf{u}) \in \mathbb{R}^d \times \mathbb{R}$:

$$ \mathrm{Cov}\{Z_i(\mathbf{s}, t), Z_j(\mathbf{s} + \mathbf{h}, t + \mathbf{u})\} = C_{ij}(\mathbf{h}, \mathbf{u}), \quad i, j = 1, \dots, p, \quad \mathbf{s}, \mathbf{s} + \mathbf{h} \in D, \quad t, t + \mathbf{u} \in T. $$

Cross-covariance functions are not symmetric, but they are invariant with respect to the joint exchange of the variables and the sign of the lag $\mathbf{k} : C_{ij}(\mathbf{k}) = C_{ji}(-\mathbf{k})$, for all $\mathbf{k} \in \mathbb{R}^{d+1}$. Full symmetry is a more restrictive assumption for which the following relationships are also
verified: $C_{ij}(\mathbf{k}) = C_{ij}(-\mathbf{k}) = C_{ji}(\mathbf{k})$. For a complete review on multivariate random fields modeling, with a particular focus on spatial cross-covariance functions, we refer the reader to Genton and Kleiber (2015) with the associated discussions.

Our goal is to elaborate valid, flexible parametric classes of space-time matrix-valued covariance functions for $Z(\mathbf{x})$, i.e. matrix-valued covariance functions that verify the well-known requirement of non-negative definiteness: for any $n \in \mathbb{N}$, for any finite set of points $(\mathbf{s}_1, t_1), ..., (\mathbf{s}_n, t_n)$ and for any vector $\mathbf{\lambda} \in \mathbb{R}^{np}$, we have $\mathbf{\lambda}^\top \mathbf{\Sigma} \mathbf{\lambda} \ge 0$, where $\mathbf{\Sigma}$ is a $np \times np$ matrix with $n \times n$ block elements of $p \times p$ matrices $C(\mathbf{s}_\alpha - \mathbf{s}_\beta, t_\alpha - t_\beta)$, with $\alpha, \beta = 1, ..., n$.

Following Genton and Kleiber (2015) and Gelfand and Banerjee (2010), a multivariate space-time covariance model is said to be separable when it is obtained through the product of a $p \times p$ covariance matrix $\mathbf{A} = [A_{ij}]_{i,j=1}^p$ and a valid univariate space-time correlation function $\rho_{ST}(\cdot)$ on $D \times T$,

$$ \mathbf{C}(\mathbf{k}) = \mathbf{A}\rho_{ST}(\mathbf{k}), \quad \mathbf{k} \in D \times T. \qquad (1) $$

When $\rho_{ST}(\cdot)$ is also space-time separable, the covariance matrix of $Z$ reduces to

$$ \mathbf{\Sigma} = \mathbf{A} \otimes \mathbf{C}_S \otimes \mathbf{C}_T, \qquad (2) $$

where $\otimes$ is the Kronecker product and $\mathbf{C}_S, \mathbf{C}_T$ are respectively spatial and temporal covariance matrices (Cressie and Wikle, 2011) associated to univariate spatial and temporal covariance functions. Separability induces reduced number of parameters and faster computation of the inverse and of the determinant of the matrix $\mathbf{\Sigma}$. Often, separability is an overly simplified assumption for weather and climate data. Space-time separability is equivalent to conditional independence between $Z(\mathbf{s}, t)$ and $Z(\mathbf{s}', t')$ given $Z(\mathbf{s}', t)$ (or $Z(\mathbf{s}, t')$), since in this case we have

$$ \mathrm{Var}\{Z(\mathbf{s}, t)\} \mathrm{Cov}\{Z(\mathbf{s}, t), Z(\mathbf{s}', t')\} = \mathrm{Cov}\{Z(\mathbf{s}, t), Z(\mathbf{s}', t)\} \mathrm{Cov}\{Z(\mathbf{s}', t), Z(\mathbf{s}', t')\}, $$
(s, t), (s', t') ∈ D × T. As a consequence, there is a proportional relationship between C(h, u) and C(h', u) for two fixed spatial lags h, h' ∈ R^d, u ∈ R, which implies that separable co-variances cannot capture sophisticated interactions between space and time. Moreover, in a multivariate framework, separability between variables and space-time variations implies that all variables are characterized by the same space-time correlation function. An approach based on the asymptotic distribution of the sample cross-covariance estimator to test for separability and for full symmetry is proposed in Li et al. (2008). This allows the practitioner to select among the important dependence structures and to make the appropriate modeling choice. Often, separability must be rejected and models that do not separate space, time and variable index must be defined.

Gneiting (2002) proposed a class of univariate fully symmetric space-time covariances, that has now become the standard class of models for univariate space-time Gaussian random fields in geostatistical applications; see also Gneiting et al. (2007) and references therein. The Gneiting class of space-time covariances is defined as

$$ \mathcal{G}(\mathbf{h}, u) = \frac{\sigma^2}{\psi(u^2)^{d/2}} \varphi\left(\frac{\|\mathbf{h}\|^2}{\psi(u^2)}\right), \quad (\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}, \qquad (3) $$

where $\varphi: [0, \infty) \to \mathbb{R}$ is completely monotone on the positive real line, with $\varphi(0) < \infty$, $\psi$ is a positive function whose derivative is completely monotone on the positive real line and $\sigma^2$ is a variance parameter. See Porcu and Zastavnyi (2011) for more relaxed necessary conditions and the complete characterization of this class. Figure 1 depicts an example of a Gneiting space-time covariance function.

In spatial context, multivariate models for which space and variables are not separable are easily built by combining separable ones. In the linear model of coregionalization (LMC, Goulard and Voltz (1992); Wackernagel (2003)), the matrix-valued covariance function is defined as a linear combination of separable models: $C(\mathbf{h}) = \sum_{k=1}^{K} \mathbf{B}_k \rho_k(\mathbf{h})$, where $\rho_k(\mathbf{h})$ are spatial covariance functions with $\rho_k(0) = 1$ and $\mathbf{B}_k$ are positive definite matrices. Drawbacks
<figure><img src="image_1.png" /><figcaption>Figure 1: An example of a Gneiting space-time covariance function $\mathcal{G}(\mathbf{h}, u) = (0.8|u|^{1.5} + 1)^{-1} \exp \{-\|\mathbf{h}\|/(0.8|u|^{1.5} + 1)^{-b/2}\}$. The case $b=0$ corresponds to separability.</figcaption></figure>

of this construction have been discussed in Gneiting et al. (2010) and in Daley et al. (2014). De Iaco et al. (2013) proposed a space-time extension of this construction. A major drawback of this approach is that the smoothness of any component of the multivariate random field is that of the roughest process associated to $\rho_k(\cdot)$, $k = 1, \dots, K$. Apanasovich et al. (2012) introduced a valid parametric family of cross-covariance functions for multivariate spatial random fields where each component has a covariance function from a Matérn class.

In space-time context, Apanasovich and Genton (2010) proposed valid non-separable matrix-valued space-time covariance functions through additional latent dimensions, one for each of the $p$ variables, that represent the variables to be modeled. These dimensions are then used in a similar way as that of time within the Gneiting class of covariance functions. Climate and weather variables need flexible space-time models able to capture different regularity properties from one variable to another. Under the framework proposed in Apanasovich and Genton (2010), this can be obtained through a LMC-type construction, at the cost of a significant increase of the number of parameters. In order to keep a reasonably low number of parameters, we argue in favor of a different approach.

The Matérn family (Matérn, 1986) of covariance functions has become extremely popular
in geostatistical applications. Its expression is

$$ \mathcal{M}(\mathbf{h}; r, \nu) = \frac{\sigma^2 2^{1-\nu}}{\Gamma(\nu)} (r\|\mathbf{h}\|)^\nu \mathcal{K}_\nu(r\|\mathbf{h}\|), \quad \mathbf{h} \in \mathbb{R}^d, $$

where $r > 0$ is a scale parameter ($1/r$ is often called the range) and $\mathcal{K}_\nu$ denotes the modified Bessel function of the second kind of order $\nu$ (DLMF, 2015). The smoothness parameter $\nu > 0$ is directly related to the regularity of the underlying random field. Specifically, a random field $Y(\mathbf{x})$ with Matérn covariance is $m$ times mean square differentiable if and only if $\nu \ge m$. As $\nu \to \infty$, the random field is infinitely mean square differentiable and its associated covariance function tends to the so-called Gaussian covariance function $C(\mathbf{h}) = \exp(-r^2\|\mathbf{h}\|^2)$, $\mathbf{h} \in \mathbb{R}$. Small values of $\nu$ yield rougher random fields; in particular $\nu = 1/2$ corresponds to the exponential covariance function $C(\mathbf{h}) = \exp(-r\|\mathbf{h}\|)$ that is mean square continuous but not mean square differentiable at the origin. When $\nu = k + 1/2$ and $k$ is a positive integer, the Matérn model reduces to the product of the negative exponential with a polynomial of degree $k$.

A multivariate version of the Matérn family has been proposed in Gneiting et al. (2010) and further extended in Apanasovich et al. (2012). These constructions have the same general structure with

$$ C_{ij}(\mathbf{h}) = \sigma_i \sigma_j \rho_{ij} \mathcal{M}(\mathbf{h}; r_{ij}, \nu_{ij}), \quad \mathbf{h} \in \mathbb{R}^d, \qquad (4) $$

with different restrictions on the parameters $\{r_{ij}, \nu_{ij}, \rho_{ij}\}_{i,j=1,...,p}$ to ensure the validity of the matrix-valued covariance function $C(\mathbf{h}) = [C_{ij}(\mathbf{h})]_{i,j=1}^p$. In this model, each marginal component of $\mathbf{Z}$, i.e. each covariance function $C_{ii}(\mathbf{h})$, has a different smoothness parameter, while allowing some cross correlations between the variables. Sufficient conditions have been obtained using scalar mixtures (Gneiting et al., 2010; Schlather, 2010; Porcu and Zastavnyi, 2011) while necessary and sufficient conditions are provided in Gneiting et al. (2010) for $p = 2$ only.

As an alternative to the Matérn class, it is also possible to build a multivariate version of
the generalized Cauchy family, having expression

$$C(\mathbf{h}; r, \beta, \lambda) = (1 + r \|\mathbf{h}\|^\lambda)^{-\beta/\lambda}, \quad \mathbf{h} \in \mathbb{R}^d,$$

with $0 < \lambda \le 2$ and $\beta, r > 0$. In Gneiting and Schlather (2004) it is shown that $\lambda$ characterizes the roughness of the associated random field. Higher values of $\lambda$ yielding smoother fields, while $\beta$ parametrizes the dependence at large distances. Another parameterization, which will be used in the rest of this work, is obtained by replacing $\beta/\lambda$ by $\nu > 0$. The isotropic Matérn and Cauchy covariance functions are related to the completely monotone functions of Table 1 in Gneiting (2002) by the equation $C(\mathbf{h}) = \varphi(\|\mathbf{h}\|^2)$ (Schoenberg, 1938).

In this article we propose a class of non-separable multivariate space-time covariance models for any number of variables. The general structure of this class is $C(\mathbf{h}, u) = [C_{ij}(\mathbf{h}, u)]_{i,j=1}^p$ with

$$C_{ij}(\mathbf{h}, u) = \frac{\sigma_i \sigma_j}{\psi(u^2)^{d/2}} \rho_{ij} \varphi_{ij} \left( \frac{\|\mathbf{h}\|^2}{\psi(u^2)} \right), \quad 1 \le i, j \le p, \quad (\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}.$$

It is an extension of the Gneiting class $\mathcal{G}$ to the multivariate case, where the completely monotone functions $\varphi_{ij}$ have different parameters for each variable.

In particular, we offer sufficient conditions for two cases: $\varphi_{ij}(t) = \mathcal{M}(t^2; r_{ij}, \nu_{ij})$ and $\varphi_{ij}(t) = C(t^2; r_{ij}, \nu_{ij}, \lambda)$. The parameterization is inspired by the multivariate Matérn structure in Gneiting et al. (2010) and Apanasovich et al. (2012) for spatial Gaussian random fields. For both classes, each variable has its own range and own degree of smoothness, while allowing for cross-correlation. The remainder of the article is organized as follows. Section 2 presents the new multivariate space-time classes. Sufficient conditions that provide valid matrix-valued covariance functions are established and proved in the Appendix. An estimation technique, based on pairwise composite likelihood, is then presented in Section 3. Section 4 is devoted to a simulation study in order to validate the estimation procedure. In Section 5, the Matérn model is applied to a multivariate space-time dataset of weather variables in France. It is shown that the proposed model offers substantial improvements to separable
models and also to more parsimonious models with equal range and smoothness parameters for all variables.

## 2 A class of non-separable space-time cross-covariance functions

The mixture approach described in Gneiting et al. (2010) to build Matérn matrix-valued covariance functions for multivariate random fields is based on the stability properties of positive definite functions (Reisert and Burkhardt, 2007, Proposition 3.1). It uses the scale mixture representation of the Matérn covariance function (Gradshteyn and Ryzhik, 2007, Eq. 3.471.9). We extend it to multivariate space-time random fields to build valid matrix-valued covariance functions in $\mathbb{R}^d \times \mathbb{R}$. It can also be adapted to other covariance functions allowing mixture representations, such as the Cauchy model. Let us first recall the mixture representation in the space-time framework.

**Lemma 1.** Let the matrix-valued function $\mathbf{m}(\xi) : Q \subset \mathbb{R} \to \mathbb{R}^p \times \mathbb{R}^p$ be symmetric and non-negative definite for all $\xi \in Q$. Let $C_\xi : \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}$ be a univariate space-time covariance function for any fixed $\xi \in Q$. Suppose furthermore that for all $i, j = 1, \dots, p$ and fixed $(\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}$ the product $m_{ij}(\cdot)C_\xi(\cdot)$ is integrable with respect to a positive measure $F$ on $Q$. Then, the matrix-valued function $C(\mathbf{h}, u) = [C_{ij}(\mathbf{h}, u)]_{i,j=1}^p$ defined through

$$C_{ij}(\mathbf{h}, u) = \int_Q m_{ij}(\xi) C_\xi(\mathbf{h}, u) F(d\xi), \quad 1 \le i, j \le p, \qquad (5)$$

is a valid p-variate matrix-valued covariance function on $\mathbb{R}^d \times \mathbb{R}$.

This lemma is a direct application of Theorem 1 in Porcu and Zastavnyi (2011). With appropriate choices for $C_\xi(\cdot)$ and $\mathbf{m}(\cdot)$ the following results can be established.

**Theorem 1.** Let $\psi(t), t \ge 0$, be a positive function with a completely monotone derivative.

The multivariate Gneiting-Matérn space-time model denoted $C^M = [C_{ij}^M(\cdot, \cdot)]_{i,j=1}^p$, with

$$C_{ij}^M(\mathbf{h}, u) = \frac{\sigma_i \sigma_j}{\psi(u^2)^{d/2}} \rho_{ij} M \left( \frac{\mathbf{h}}{\psi(u^2)^{1/2}}; r_{ij}, \nu_{ij} \right), \quad (\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}, \qquad (6)$$
defines a valid matrix-valued covariance function if, for all $i, j = 1, \dots, p$,

$$r_{ij} = \{(r_i^2 + r_j^2)/2\}^{1/2},$$

$$\nu_{ij} = (\nu_i + \nu_j)/2,$$

$$\rho_{ij} = \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2}\Gamma(\nu_j)^{1/2}} \frac{r_i^{\nu_i} r_j^{\nu_j}}{r_{ij}^{2\nu_{ij}}},$$

with $r_i, \nu_i > 0$ for all $i = 1, \dots, p$, and where $\beta = [\beta_{ij}]_{i,j=1}^p$ is a correlation matrix.

**Theorem 2.** Let $\psi(t), t \ge 0$, be a positive function with a completely monotone derivative.

The multivariate Gneiting-Cauchy space-time model denoted $\mathcal{C} = [C_{ij}^{\mathcal{C}}(\cdot, \cdot)]_{i,j=1}^p$, with

$$C_{ij}^{\mathcal{C}}(h, u) = \frac{\sigma_i \sigma_j}{\psi(u^2)^{d/2}} \rho_{ij} \mathcal{C} \left( \frac{h}{\psi(u^2)^{1/2}}; r_{ij}, \nu_{ij}, \lambda \right), \quad (h, u) \in \mathbb{R}^d \times \mathbb{R}, \quad (7)$$

is a valid matrix-valued covariance function if, for all $i, j = 1, \dots, p$,

$$r_{ij} = \{(r_i^{-1} + r_j^{-1})/2\}^{-1},$$

$$\nu_{ij} = (\nu_i + \nu_j)/2,$$

$$\rho_{ij} = \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2}\Gamma(\nu_j)^{1/2}} \frac{r_{ij}^{\nu_{ij}}}{(r_i^{\nu_i} r_j^{\nu_j})^{1/2}},$$

with $r_i, \nu_i > 0$ for all $i = 1, \dots, p$ and $0 < \lambda \le 2$, and where $\beta = [\beta_{ij}]_{i,j=1}^p$ is a correlation matrix.

The proofs are deferred to the Appendix. Each marginal space-time covariance function $C_{ii}^{\mathcal{F}}(\cdot, \cdot)$, $i = 1, \dots, p$, $\mathcal{F} \in \{\mathcal{M}, \mathcal{C}\}$, belongs to the Gneiting class with different regularity parameters $\nu_i$, $i = 1, \dots, p$ for its associated Matérn, respectively Cauchy, spatial covariance function. The Gneiting-Matérn model is a space-time generalization of one particular case of the multivariate spatial “flexible model” presented in Apanasovich et al. (2012). Figure 2 shows a realization of a bivariate Gaussian field with Gneiting-Matérn space-time covariance at instants $t \in \{0, 1, 2\}$. In this example, the first component is smooth, with $\nu_1 = 1.5$ and scale parameter $r_1 = 1$, whereas the second component is rough, with $\nu_2 = 0.5$ and $r_2 = 0.5$. The co-located correlation parameter is $\rho_{12} = 0.5$.
<figure><img src="image_2.png" /><figcaption>Figure 2: Bivariate Gneiting-Matérn space-time Gaussian random field at 3 consecutive instants. Top row: smoother component ($\nu_1 = 1.5$ and $r_1 = 1$); bottom row: rougher component ($\nu_2 = 0.5$ and $r_2 = 0.5$). The co-located correlation parameter $\rho_{12}$ is set to 0.5.</figcaption></figure>
Positive functions with completely monotone derivative, sometimes referred to as Bernstein functions (Porcu and Schilling, 2011), are related to isotropic variograms by the following relationship: if a function $\psi(x)$, $x \ge 0$, belongs to the class of Bernstein functions, the function $\gamma(\mathbf{h}) = \psi(\|\mathbf{h}\|^2) - c, \mathbf{h} \in \mathbb{R}^d$ is a valid isotropic variogram for some constant $c$. For the rest of this work, we will consider the parametric form advocated in Gneiting (2002) and Schlather (2010)

$$ \psi(x) = (\alpha x^a + 1)^b, \quad x \ge 0, $$

with $\alpha > 0, 0 < a \le 1, 0 \le b \le 1$. With this choice (and a slight abuse of notation, since the Matérn covariance does not have a parameter $\lambda$), the above models become

$$ C_{ij}(\mathbf{h}, u) = \frac{\sigma_i \sigma_j}{(\alpha |u|^{2a} + 1)^{bd/2}} \rho_{ij} \mathcal{F} \left( \frac{\mathbf{h}}{(\alpha |u|^{2a} + 1)^{b/2}}; r_{ij}, \nu_{ij}, \lambda \right), \quad (\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}, \quad (8) $$

with $\mathcal{F} \in \{\mathcal{M}, \mathcal{C}\}$ and the parameter restrictions stated in Theorems 1 and 2. Note that in the parameterization (8), the space-time non-separability parameter $b$ acts both on the spatial covariance $\mathcal{F}$ and on the temporal one. A major drawback of this parameterization is that in case of space-time separability, i.e. when $b = 0$, the temporal covariance is equal to 1 for all $u \in \mathbb{R}$. Following Gneiting (2002) a reparameterization is useful. Multiplying (8) by the temporal covariance function $(\alpha|u|^{2a} + 1)^{-\delta/2}, u \in \mathbb{R}$, with $\delta > 0$ and replacing the exponent $\delta + bd/2$ by $\tau$ leads to the parametric family, which will be used in the rest of this work

$$ C_{ij}(\mathbf{h}, u) = \frac{\sigma_i \sigma_j}{(\alpha |u|^{2a} + 1)^{\tau}} \rho_{ij} \mathcal{F} \left( \frac{\mathbf{h}}{(\alpha |u|^{2a} + 1)^{b/2}}; r_{ij}, \nu_{ij}, \lambda \right), \quad (\mathbf{h}, u) \in \mathbb{R}^d \times \mathbb{R}, \quad (9) $$

with $\tau \ge bd/2$. This parameterization provides independent and interpretable parameters for the spatial scale and smoothness governed respectively by the parameters $r_i$ and $\nu_i$, the temporal scale $\alpha$, smoothness $a$ and the separability $b$. Note that the above class of models is fully symmetric, i.e. $C_{ij}(\mathbf{h}, u) = C_{ji}(\mathbf{h}, u) = C_{ij}(\mathbf{h}, -u) = C_{ij}(-\mathbf{h}, u)$ for all $i, j = 1, \dots, p$.
## 3 Estimation

For a Gaussian random field, maximum full likelihood (FL) requires the evaluation of the determinant and the inverse of the $np \times np$ covariance matrix, where $n$ is the number of space-time points and $p$ the number of variables. The computational cost being of the order of $O((np)^3)$ for both operations, a maximum full likelihood is unfeasible for large datasets, and even for datasets of moderate size as soon as the number of variables increases.

Composite likelihood (CL) methods (Lindsay, 1988) have been proposed to perform efficient estimation with less time-consuming computational steps, and to guarantee good asymptotic properties. Composite likelihoods are products of smaller likelihoods defined on certain subsets of data such as marginal or conditional events, which are easy to compute. See Varin et al. (2011) for an overview. Following Bevilacqua and Gaetan (2015), we choose the pairwise marginal Gaussian likelihoods computed on all pairs of data $\{Z_i(s_\alpha, t_\alpha), Z_j(s_\beta, t_\beta)\}$, where $i, j = 1, \dots, p$ and $\alpha, \beta = 1, \dots, n$, $n$ being the number of sites. The negative log-likelihood of such pairs is

$$l(i, j, s_\alpha, s_\beta, t_\alpha, t_\beta; \theta) = \frac{1}{2} \left\{ \log \Delta_{ij,\alpha\beta} + \frac{A_{ij,\alpha\beta}}{\Delta_{ij,\alpha\beta}} \right\}$$

where $\Delta_{ij,\alpha\beta} = \sigma_i^2 \sigma_j^2 - C_{ij}(h, u)^2$ and $A_{ij,\alpha\beta} = \sigma_j^2 Z_i(s_\alpha, t_\alpha)^2 - 2C_{ij}(h, u)Z_i(s_\alpha, t_\alpha)Z_j(s_\beta, t_\beta) + \sigma_i^2 Z_j(s_\beta, t_\beta)^2$, with $h = \|s_\alpha - s_\beta\|$, $u = |t_\alpha - t_\beta|$. In the special case of multivariate symmetric models, the Weighted Pairwise log-Likelihood (WPL) is thus

$$\text{wpl}(\theta) = \sum_{(i,j,\alpha,\beta) \in \Lambda} l(i, j, s_\alpha, s_\beta, t_\alpha, t_\beta; \theta) w_{\alpha\beta}, \quad \theta \in \Theta, \qquad (10)$$

where

$$\Lambda = \{i = j = 1, \dots, p, \quad \alpha = 1, \dots, n-1, \quad \beta = \alpha + 1, \dots, n\} \\ \cup \quad \{i = 1, \dots, p-1, \quad j = i+1, \dots, p, \quad \alpha, \beta = 1, \dots, n\},$$

and $\Theta \subset \mathbb{R}^q$ is the space of parameters for which the model in (9) is valid, with $q = (p + 2)(p+3)/2$ for the Gneiting-Matérn class and $q = (p+2)(p+3)/2+1$ for the Gneiting-Cauchy class.
The computational cost of WPL is of the order of $O((np)^2)$ when considering all possible pairs. It can be significantly reduced with an adequate choice of the weights. In this paper, we have chosen cut-off weights, that is $w_{\alpha\beta} = 1$ if the space-time lag $(\|h\|, u) \le (d_S, d_T)$, where the order is defined pointwise. Otherwise we set $w_{\alpha\beta} = 0$. We use the expression *weighted pairwise likelihood* but it is sometimes called *truncated composite likelihood* or *tapered composite likelihood*. In spite of its simplicity, this weighting scheme provides a significant gain in computational time, and preserves a reasonable level of statistical efficiency, since pairs of observations whose distance is beyond the correlation range are uninformative for the smoothness and range parameters of the covariance function. The choice of $d = (d_S, d_T)$ will be discussed in details in the next section. At this stage, we can notice that $d$ is the same for all variables.

In the univariate setting, Bevilacqua et al. (2012) proposed to seek the “optimal” window $d^*$ such that

$$d^* = \arg\min_{d} \operatorname{tr}(G_{np}^{-1}(d; \theta)), \quad \theta \in \Theta, \qquad (11)$$

where $\theta$ is replaced by a consistent estimator of $\theta$, such as for example the weighted least squares estimate based on empirical estimation of the variogram (Cressie, 1993, p. 91). The value $d^*$ minimizes the sum of the diagonal elements of the inverse of $G_{np}^{-1}(d; \theta)$, which is the asymptotic variance-covariance matrix of the WPL estimators. The matrix $G_{np}(d; \theta)$ is the Godambe information matrix, also referred to as the “sandwich” information matrix (Bevilacqua et al., 2012). The Godambe information matrix of a full log-likelihood function is equal to the Fisher information matrix.

Estimating the parameters from model (9) requires to maximize wpl($\theta$) in $\Theta$, when wpl($\theta$) is computed with the appropriate window $d^*$. Closed form expressions for the theoretical expression of the Godambe matrix are not easily feasible since first and second order derivatives of (10) with respect to all parameters are required. We thus favor an empirical approach, based on a simulation study detailed in the next Section. For a given value of $d^*$, maximiz-
ing wpl($\theta$) in $\Theta$ is by no means trivial, since model in (9) has $(p+2)(p+3)/2$ parameters.
Employing blindly an optimization function is deemed to fail, even for the simplest case with
$p=2$.

As a way to alleviate this problem, wpl($\theta$) is maximized sequentially in subspaces of $\Theta$, corresponding to blocks of related parameters, while keeping all other parameters fixed to the previously attained values. Among the schemes that have been tested, the following combination seemed to lead to the best estimates: for a fixed value of $b$, first seek the maximum corresponding to the $p(p-1)/2$ correlation parameters $\beta_{ij}$; then, for each variable $i$, find the maximum for $\sigma_i$, followed by $(\nu_i, r_i)$; finally, maximize with respect to the two temporal parameters $a$ and $\alpha$. One iteration is achieved when all parameters have been estimated once. The likelihood is iteratively maximized until a stopping criterion is reached. We have chosen to stop the maximization whenever the log-likelihood is increased by less than one unit. Other orders have been tested, leading to very close estimates.

The separability parameter $b$, crucial in this study, is difficult to optimize in the above framework. We often observed unstable maximization and/or estimates reaching the boundary of the parameter space $\Theta$. Since $b$ belongs to the interval $[0, 1]$, we decided to estimate it by a simple grid search. The parameter $b$ is successively fixed to the values $k/10$, with $k = 0, \dots, 10$. For each of these values, the above maximization procedure is performed. Finally $\hat{b}$ is that value of $b$ corresponding to the highest maximized likelihood.

To perform the maximization within the subspaces of $\Theta$, we used the function *optim* with quasi Newton method "BFGS" (published simultaneously in 1970 by Broyden, Fletcher, Goldfarb and Shanno) as implemented in R. Relevant initial values were set by estimating parameters of the corresponding marginal spatial and temporal covariance functions, independently for each variable. Initial values of the correlation coefficients were set to their empirical values.
# 4 Simulation study

We propose a simulation study with three main goals: setting the optimal window $d^*$, assessing the performance of the estimation procedure detailed in Section 3, and measuring the gain in prediction implied by model (9), as compared to non-separable or less flexible models. Because of the wide popularity of the Matérn covariance functions, and also for space considerations, we only report results obtained with the multivariate Gneiting-Matérn class (6) in Theorem 1. The setting of the simulation study mimics the conditions found for the analysis of weather data as detailed in Section 5. Three variables are available during one month (30 days) at 13 locations (the precise locations are shown later, in Figure 4). Among these, 11 sites are used for the first two tasks (circles) while two sites are kept for validation (stars). There are 10 years of data, assumed independent. We thus simulate a sample of size 10 from a Gaussian vector with $30 \times 13 \times 3 = 1170$ components. The elements of the $1170 \times 1170$ matrix are given by the multivariate space-time Matérn covariance function with unit variance for all variables, i.e.

$$C_{ij}(\mathbf{h}, u) = \frac{\rho_{ij}}{\alpha|u|^{2a} + 1} \mathcal{M} \left( \frac{\mathbf{h}}{(\alpha|u|^{2a} + 1)^{b/2}}; r_{ij}, \nu_{ij} \right) \quad i, j = 1, \dots, p. \quad (12)$$

Exact simulations with Cholesky decomposition are feasible, and a total of 100 repetitions are simulated. The reasonable dimension of the matrices ($1170 \times 1170$) allows to compare WPL inference with a usual likelihood approach, referred to as the Full Likelihood (FL) approach. Parameters have been chosen to match those estimated by an exploratory analysis of the spatial and temporal margins of the dataset (see Table 3). In order to account for lack of space-time separability, the parameter $b$ was set equal to 0.8. As already pointed out earlier, we only report results obtained with the Matérn model.

## 4.1 Setting the optimal window

The optimal window $d^*$ given in Eq. (11) minimizes the sum of the estimation variance of every parameter. For several choices of the spatial distance $d_S \in \{250, 500, 750\}$ (expressed
<table><caption>Table 1: Average of the sum of the estimated variances of all parameters for several windows $d = (d_S, d_T)$. 100 synthetic datasets simulated according to model (12). The number in brackets corresponds to the percentage of pairs used in the computation of wpl.</caption><thead><tr><th><b>d</b></th><th>d<sub>T</sub> = 2 days</th><th>d<sub>T</sub> = 5 days</th><th>d<sub>T</sub> = 10 days</th></tr></thead><tbody><tr><td>d<sub>S</sub> = 250 km</td><td><b>0.120</b> (7.2%)</td><td>0.148 (15.1%)</td><td>0.185 (26.2%)</td></tr><tr><td>d<sub>S</sub> = 500 km</td><td>0.127 (13.5%)</td><td>0.162 (28.3%)</td><td>0.203 (49.1%)</td></tr><tr><td>d<sub>S</sub> = 750 km</td><td>0.146 (15.4%)</td><td>0.196 (33.3%)</td><td>0.227 (57.7%)</td></tr></tbody></table>

in kms), and temporal distance $d_T \in \{2, 5, 10\}$ (in days), the estimation of the parameters is performed on each simulated repetition. Then, empirical estimates of the estimation variances are computed. Their sum is an estimate of the criterion in Eq. (11), as recommended in Bevilacqua et al. (2012). Results are reported in Table 1. As the spatial distance increases, the quantity $G_{np}^{-1}(d; \hat{\theta})$ increases for all temporal distances. For all spatial distances, the criterion is minimal for the smallest temporal window. Parameters are best estimated if pairs corresponding to moderate distances and to successive days only are used, which corresponds to 7.2% and 13.5% of the number of pairs for the two lowest values of the criterion. Overall, the window $d = (250\text{km}, 2\text{days})$ is optimal according to this criterion. The criterion computed with $d = (500\text{km}, 2\text{days})$ is very close to the optimal. The estimated mean and variance of the separability parameter $b$ for the different space-time windows, is reported in Table 2. The picture is now slightly different. Bias and estimation variance are minimal for $d = (500\text{km}, 2\text{days})$. Considering the importance of the separability parameter, and since the trace criterion for this window is very close to the optimum, the window is set to $d = (500\text{km}, 2\text{days})$ for the rest of this work. These results are consistent with those reported in Bevilacqua et al. (2012), where it is observed that the relative efficiency of WPL as a function of the distance first increases to a maximum and then decreases as more distant pairs are added in the WPL.
<table><caption>Table 2: Estimated mean and standard deviation (in brackets) of the separability parameter b for several windows d = (d<sub>S</sub>, d<sub>T</sub>). Same simulations as in Table 1.</caption><thead><tr><th>d</th><th>d<sub>T</sub> = 2 days</th><th>d<sub>T</sub> = 5 days</th><th>d<sub>T</sub> = 10 days</th></tr></thead><tbody><tr><td>d<sub>S</sub> = 250 km</td><td>0.694 (0.234)</td><td>0.681 (0.268)</td><td>0.630 (0.305)</td></tr><tr><td>d<sub>S</sub> = 500 km</td><td><b>0.768</b> (<b>0.230</b>)</td><td>0.731 (0.251)</td><td>0.717 (0.259)</td></tr><tr><td>d<sub>S</sub> = 750 km</td><td>0.744 (0.245)</td><td>0.723 (0.257)</td><td>0.696 (0.281)</td></tr></tbody></table>

## 4.2 Efficiency of the estimation procedure

We now assess the ability of the estimation procedure to estimate efficiently the different parameters, the window being set to d = (500km, 2days). This simulation study also allows us to compare the efficiency of WPL with respect to FL for the multivariate Gneiting-Matérn class. We maximize WPL and FL, using the same maximization scheme, on the same set of 100 simulations. As can be observed from Table 3, the overall performances are good for both approaches. For most parameters, the difference between the median and the mean is negligible, and differences between the mean and the true value is relatively small for most parameters. The range and smoothness parameters are a notable exception. Their estimators are less biased and more dispersed when using WPL than when using FL. For all variables $i \in \{1, 2, 3\}$, the smoothness parameters, $\nu_i$, tend to be overestimated, while the scale parameters, $1/r_i$, tend to be underestimated. In Zhang (2004), it is shown that these parameters compensate each other and that the simultaneous estimation of both is difficult. Results shown in Table 3 confirm our assertions. In accordance with the parameterization chosen in Theorems 1 and 2, the correlation coefficients $\beta_{ij}$, with $1 \le i < j \le 3$ must belong to the interval $(-1, 1)$ with $\beta_{ii} = 1$, for $i = 1, \dots, 3$, and the matrix $[\beta_{ij}]^3_{i,j=1}$ must be positive definite. The estimators of the correlations coefficients $\beta_{ij}$, with $1 \le i < j \le 3$, are unbiased for both approaches, but their dispersion is much larger when using WPL than when using FL. We return to this point later. Remember that the separability parameter b is not estimated continuously and that the window d has been chosen to minimize the bias and
the variance of its estimator. Interestingly, it is better estimated using WPL than using FL. With WPL, the median is equal to the true value and it is close to the mean. The separability parameter is estimated to be non null, for all simulated datasets but one.

To complete this comparison, we compute the relative efficiency of WPL with respect to FL. For each approach, and for each of the 15 parameters $\theta_i$, with $i \le i \le 15$, the root mean square error, $rmse_i = \sqrt{bs_i^2 + sd_i^2}$, is computed, where the bias is $bs_i = \tilde{\theta}_i - \theta_i$, and the variance is $sd_i^2 = \sum_{j=1}^{100} (\hat{\theta}_{j,i} - \tilde{\theta}_i)^2 / 100$, with $\tilde{\theta}_i = \sum_{j=1}^{100} \hat{\theta}_{j,i} / 100$. For each variable, we denote $rmse_i^{FL}$ and $rmse_i^{WPL}$ the root mean square error corresponding respectively to FL and to WPL. The root relative efficiency, defined as

$$rre_i = rmse_i^{FL} / rmse_i^{WPL},$$

is reported for each variable in Table 4. When $rre < 1$, FL is more efficient than WPL, and conversely when $rre > 1$. For all parameters, it was found that FL leads to higher bias than WPL, but that the lower variance associated to FL usually more than compensates this effect. On these simulations, WPL is found to be rather inefficient for the estimation of the correlation parameters $\beta_{ij}$, for which the variance associated to WPL is about 30 times that associated to FL. For all other parameters and in particular for regularity and range parameters, WPL provides estimates with a mild loss in efficiency as compared to a full likelihood approach, with a significant gain in terms of computation. Similar results were obtained for the inference of the parameters of Max-Stable processes in Castruccio et al. (2015). Rarely, the relative efficiency is even slightly larger than 1 but this could be due to the moderate number of replicates. In these cases, the lower variance for FL does not compensate the larger bias.

## 4.3 Assessing the predictive performances

Finally, the predictive performance of the model is assessed. In addition to the non-separable model according to which data are simulated, we also consider three more specific models,
<table><caption>Table 3: Summary statistics of the estimated parameters for d = (500 km, 2 days). Same simulations as in Table 1.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">True</th><th colspan="6">Weighted Pairwise Likelihood</th><th colspan="6">Full Likelihood</th></tr><tr><th>Min</th><th>Q<sub>1</sub></th><th>Median</th><th>Mean</th><th>Q<sub>3</sub></th><th>Max</th><th>Min</th><th>Q<sub>1</sub></th><th>Median</th><th>Mean</th><th>Q<sub>3</sub></th><th>Max</th></tr></thead><tbody><tr><td>σ<sub>1</sub></td><td>1.00</td><td>0.89</td><td>0.97</td><td>1.00</td><td>1.00</td><td>1.02</td><td>1.16</td><td>0.93</td><td>0.97</td><td>0.99</td><td>0.99</td><td>1.01</td><td>1.06</td></tr><tr><td>σ<sub>2</sub></td><td>1.00</td><td>0.91</td><td>0.97</td><td>1.00</td><td>1.00</td><td>1.02</td><td>1.11</td><td>0.95</td><td>0.98</td><td>0.99</td><td>1.00</td><td>1.01</td><td>1.06</td></tr><tr><td>σ<sub>3</sub></td><td>1.00</td><td>0.92</td><td>0.97</td><td>0.99</td><td>0.99</td><td>1.02</td><td>1.10</td><td>0.94</td><td>0.97</td><td>0.99</td><td>0.99</td><td>1.00</td><td>1.12</td></tr><tr><td>β<sub>12</sub></td><td>-0.40</td><td>-0.61</td><td>-0.46</td><td>-0.40</td><td>-0.40</td><td>-0.34</td><td>-0.15</td><td>-0.45</td><td>-0.42</td><td>-0.41</td><td>-0.41</td><td>-0.40</td><td>-0.37</td></tr><tr><td>β<sub>13</sub></td><td>-0.40</td><td>-0.62</td><td>-0.47</td><td>-0.41</td><td>-0.41</td><td>-0.34</td><td>-0.21</td><td>-0.45</td><td>-0.42</td><td>-0.41</td><td>-0.41</td><td>-0.40</td><td>-0.36</td></tr><tr><td>β<sub>23</sub></td><td>0.25</td><td>-0.05</td><td>0.19</td><td>0.27</td><td>0.26</td><td>0.33</td><td>0.47</td><td>0.21</td><td>0.25</td><td>0.26</td><td>0.26</td><td>0.27</td><td>0.31</td></tr><tr><td>ν<sub>1</sub></td><td>0.70</td><td>0.55</td><td>0.69</td><td>0.74</td><td>0.75</td><td>0.82</td><td>1.01</td><td>0.58</td><td>0.72</td><td>0.75</td><td>0.76</td><td>0.79</td><td>0.99</td></tr><tr><td>ν<sub>2</sub></td><td>0.80</td><td>0.55</td><td>0.76</td><td>0.83</td><td>0.84</td><td>0.91</td><td>1.17</td><td>0.68</td><td>0.82</td><td>0.87</td><td>0.88</td><td>0.95</td><td>1.17</td></tr><tr><td>ν<sub>3</sub></td><td>0.40</td><td>0.28</td><td>0.39</td><td>0.43</td><td>0.42</td><td>0.46</td><td>0.72</td><td>0.30</td><td>0.43</td><td>0.45</td><td>0.45</td><td>0.47</td><td>0.59</td></tr><tr><td>1/r<sub>1</sub></td><td>250</td><td>135</td><td>205</td><td>229</td><td>235</td><td>262</td><td>388</td><td>147</td><td>199</td><td>216</td><td>220</td><td>238</td><td>333</td></tr><tr><td>1/r<sub>2</sub></td><td>200</td><td>128</td><td>170</td><td>190</td><td>192</td><td>215</td><td>299</td><td>131</td><td>158</td><td>177</td><td>178</td><td>193</td><td>261</td></tr><tr><td>1/r<sub>3</sub></td><td>350</td><td>149</td><td>284</td><td>310</td><td>327</td><td>356</td><td>1017</td><td>195</td><td>247</td><td>277</td><td>283</td><td>293</td><td>1017</td></tr><tr><td>α</td><td>0.90</td><td>0.65</td><td>0.84</td><td>0.90</td><td>0.91</td><td>0.99</td><td>1.24</td><td>0.78</td><td>0.91</td><td>0.95</td><td>0.95</td><td>1.00</td><td>1.12</td></tr><tr><td>a</td><td>0.50</td><td>0.35</td><td>0.46</td><td>0.49</td><td>0.50</td><td>0.54</td><td>0.71</td><td>0.46</td><td>0.49</td><td>0.51</td><td>0.51</td><td>0.53</td><td>0.58</td></tr><tr><td>b</td><td>0.80</td><td>0.00</td><td>0.60</td><td>0.80</td><td>0.77</td><td>1.00</td><td>1.00</td><td>0.30</td><td>0.50</td><td>0.60</td><td>0.64</td><td>0.70</td><td>1.00</td></tr></tbody></table>

<table><caption>Table 4: Bias (bs), variance (sd²), root mean square error (rmse) and root relative efficiency (rre) for each parameter. Same simulations as in Table 1.</caption><thead><tr><th rowspan="2"></th><th colspan="2">σ<sub>1</sub></th><th colspan="2">σ<sub>2</sub></th><th colspan="2">σ<sub>3</sub></th><th colspan="2">β<sub>12</sub></th><th colspan="2">β<sub>13</sub></th><th colspan="2">β<sub>23</sub></th></tr><tr><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th></tr></thead><tbody><tr><td>bs<sup>2</sup> (×10<sup>-5</sup>)</td><td>0.18</td><td>5.13</td><td>1.13</td><td>2.54</td><td>3.48</td><td>17.1</td><td>1.02</td><td>5.77</td><td>6.30</td><td>4.50</td><td>7.96</td><td>5.35</td></tr><tr><td>sd<sup>2</sup> (×10<sup>-3</sup>)</td><td>1.73</td><td>0.84</td><td>1.40</td><td>0.59</td><td>1.24</td><td>0.63</td><td>8.47</td><td>0.26</td><td>7.86</td><td>0.29</td><td>10.8</td><td>0.36</td></tr><tr><td>rmse</td><td>0.042</td><td>0.030</td><td>0.038</td><td>0.025</td><td>0.036</td><td>0.028</td><td>0.092</td><td>0.018</td><td>0.089</td><td>0.018</td><td>0.104</td><td>0.020</td></tr><tr><td>rre</td><td colspan="2">0.72</td><td colspan="2">0.66</td><td colspan="2">0.79</td><td colspan="2">0.19</td><td colspan="2">0.21</td><td colspan="2">0.19</td></tr><tr><th rowspan="2"></th><th colspan="2">ν<sub>1</sub></th><th colspan="2">ν<sub>2</sub></th><th colspan="2">ν<sub>3</sub></th><th colspan="2">α</th><th colspan="2">a</th><th colspan="2">b</th></tr><tr><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th></tr><tr><td>bs<sup>2</sup> (×10<sup>-3</sup>)</td><td>2.63</td><td>3.48</td><td>1.75</td><td>6.30</td><td>0.60</td><td>2.95</td><td>0.132</td><td>2.82</td><td>0.001</td><td>0.018</td><td>1.000</td><td>27.0</td></tr><tr><td>sd<sup>2</sup> (×10<sup>-3</sup>)</td><td>9.71</td><td>4.30</td><td>13.1</td><td>8.95</td><td>3.81</td><td>1.82</td><td>11.4</td><td>4.70</td><td>5.08</td><td>0.699</td><td>52.7</td><td>19.4</td></tr><tr><td>rmse</td><td>0.111</td><td>0.088</td><td>0.122</td><td>0.123</td><td>0.066</td><td>0.069</td><td>0.107</td><td>0.087</td><td>0.071</td><td>0.027</td><td>0.232</td><td>0.215</td></tr><tr><td>rre</td><td colspan="2">0.79</td><td colspan="2">1.01</td><td colspan="2">1.04</td><td colspan="2">0.81</td><td colspan="2">0.38</td><td colspan="2">0.93</td></tr><tr><th rowspan="2"></th><th colspan="2">1/r<sub>1</sub></th><th colspan="2">1/r<sub>2</sub></th><th colspan="2">1/r<sub>3</sub></th><th colspan="2"></th><th colspan="2"></th><th colspan="2"></th></tr><tr><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th><th>WPL</th><th>FL</th></tr><tr><td>bs<sup>2</sup></td><td>218.5</td><td>911.6</td><td>61.8</td><td>473.9</td><td>516.6</td><td>4539.3</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>sd<sup>2</sup></td><td>2300.0</td><td>856.1</td><td>1127.1</td><td>682.4</td><td>9629</td><td>6788.9</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>rmse</td><td>50.2</td><td>42.0</td><td>34.5</td><td>34.0</td><td>100.7</td><td>106.4</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>rre</td><td colspan="2">0.84</td><td colspan="2">0.99</td><td colspan="2">1.06</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
which are special cases of the full model. We compare space-time separable (S) and non-separable (NS) models, and models with equal (E) or different (D) smoothness and scale parameters. We will thus consider the following 4 models:

1. S-E: $b=0$; all components have equal smoothness and scale parameters. The covariance matrix follows (2).

2. NS-E: $0 < b \le 1$; all components have equal smoothness and scale parameters. This model is space-time non-separable, but it has the same space-time covariance for all variables. The covariance matrix-valued function follows (1).

3. S-D: $b=0$; each component has its own smoothness and scale parameters, but it is space-time separable.

4. NS-D: $0 < b \le 1$; each component has its own smoothness and scale parameters, and it is not space-time separable. It is the non-separable model defined through Eq. (12).

Data are simulated according to model NS-D, with parameters as in Table 3, at the 13 stations displayed in Figure 4. Two stations, indicated with a star, have been selected for validation. Data at these two stations are not used to estimate the parameters. For each day $t$, the conditional distribution is computed at these two locations, given the sets $\{\text{all data at time } t-1, \dots, t-q\} \cup \{\text{data at other locations at the same time } t\}$, where $q \in \{2, 3, \dots\}$ is the number of days used for the prediction. Under Gaussianity, this amounts to compute the 6-variate conditional expectation and the conditional covariance matrix.

The four models are compared by means of four different scores: Root Mean Square Error (RMSE), Mean Absolute Error (MAE), the Continuous Ranked Probability Score (CRPS) and Logarithmic Score (LogS) (Gneiting and Raftery, 2007). Let use denote $\tilde{z}_1^t, \dots, \tilde{z}_6^t$ the conditional expectation of the 6 predicted variables ($2$ sites $\times$ $3$ variables) for a given day $t$, $\tilde{\sigma}_1, \dots, \tilde{\sigma}_6$ the corresponding conditional standard deviations, $z_1^t, \dots, z_6^t$ being the observed
values. Note that the conditional standard deviations are independent of the day t, since they only depend on covariance values.

The first two scores, MAE and RMSE, compare the conditional expectation to the true value. The MAE is defined as

$$MAE = \frac{1}{6|T|} \sum_{t \in T} \sum_{j=1}^{6} |\tilde{z}_{j}^{t} - z_{j}^{t}|,$$

where the sum is taken over the set of all testing days, T and |T| denotes the cardinality of the set T. The mean square error (MSE) is

$$MSE = \frac{1}{6|T|} \sum_{t \in T} \sum_{j=1}^{6} (\tilde{z}_{j}^{t} - z_{j}^{t})^{2}.$$

The RMSE is the square root of the MSE and has the advantage of being recorded in the same unit as the data.

The other two scores, CRPS and LogS, assess not only the prediction but its variance as well. They are easily computed in the case of a normal predictive distribution. The CRPS measures the discrepancy between the predictive cumulative distribution function (CDF) and the true value. Specifically, if F is the predictive CDF and z the true value, crps is defined as

$$crps(F, z) = \int_{-\infty}^{\infty} [F(u) - H(u - z)]^2 du,$$

where H(u - z) denotes the Heaviside function, which takes the value 0 when u < z and 1 otherwise. In case of normal CDF Φ<sub>j</sub><sup>t</sup> with expectation ˜z<sub>j</sub><sup>t</sup> and variance ˜σ<sub>j</sub><sup>2</sup>, one can show that

$$crps(\Phi_j^t, z) = \tilde{\sigma}_j^t \{ \tilde{z}[2\Phi(\tilde{z}) - 1] + 2\phi(\tilde{z}) - \pi^{-1/2} \},$$

where ˜z = (z - ˜z<sub>j</sub><sup>t</sup>)/ ˜σ<sub>j</sub> is the normalized prediction error and where φ(·) and Φ(·) denote the density and the CDF, respectively, of the normal distribution with mean 0 and variance 1. The CRPS is

$$CRPS = \frac{1}{6|T|} \sum_{t \in T} \sum_{j=1}^{6} crps(\Phi_j^t, y_j^t).$$
It is easy to show that the CRPS tends to the MAE when $\tilde{\sigma}_j \to 0$, for all $j = 1, \dots, 6$. For this reason, the CRPS can be interpreted as a generalized version of the MAE (Gneiting and Raftery, 2007). The marginal logarithmic score is the negative of the logarithm of the marginal predictive density at the true value

$$ \text{LogS}_1 = \frac{1}{6|T|} \sum_{t \in T} \sum_{j=1}^{6} \left\{ \ln \tilde{\sigma}_j + \frac{(\tilde{z}_j^t - z_j^t)^2}{2\tilde{\sigma}_j^2} \right\}. $$

The multivariate logarithmic score considers the multivariate predictive density computed at the true vector $\mathbf{z}^t = (z_1, \dots, z_6)^\top$

$$ \text{LogS}_6 = \frac{1}{6|T|} \sum_{t \in T} \left\{ \ln \det \tilde{\Sigma}^t + \frac{(\tilde{z}^t - \mathbf{z}^t)^\top \tilde{\Sigma}^{-1} (\tilde{z}^t - \mathbf{z}^t)}{2} \right\}, $$

where $\tilde{\Sigma}$ is the conditional covariance matrix and $\mathbf{a}^\top$ is the transpose of $\mathbf{a}$.

The scores defined above are computed for the set of 100 simulations, for each of the four models, and for both estimation approaches, WPL and FL. The average scores are reported in Table 5, and the boxplots for WPL are shown in Figure 3. We only report results obtained using the two previous days as conditioning data (i.e. $q = 2$), but similar results were obtained for $q = 1$ and $q = 3$. Recall that lower scores indicate a better adequacy between the model and the data.

Some comments are in order: in general, differences in score values are low, but they are consistently observed on all simulations. Therefore, they can be considered as significant. Obviously, the highest scores correspond to the fully separable model, SE, whereas the lowest scores are obtained for the non-separable model, NS-D, for which the difference between WPL and FL is negligible. Differences between models are more pronounced when considering scores that also involve the conditional variance (CRPS and LogS). Accounting for non-separability makes a more important difference than having different smoothness and scale parameters for each variable. It is interesting to note that the difference between the models is more pronounced with WPL than with FL.
<table><caption>Table 5: Mean of the prediction scores RMSE, MAE, CRPS, LogS₁ and LogS₆, according to models S-E, NS-E, S-D and NS-D, and for WPL and FL approaches. Same simulations as in Table 1.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Weighted Pairwise Likelihood</th><th colspan="5">Full Likelihood</th></tr><tr><th>RMSE</th><th>MAE</th><th>CRPS</th><th>LogS<sub>1</sub></th><th>LogS<sub>6</sub></th><th>RMSE</th><th>MAE</th><th>CRPS</th><th>LogS<sub>1</sub></th><th>LogS<sub>6</sub></th></tr></thead><tbody><tr><td>S-E</td><td>0.494</td><td>0.391</td><td>0.278</td><td>0.725</td><td>4.070</td><td>0.489</td><td>0.386</td><td>0.274</td><td>0.690</td><td>3.792</td></tr><tr><td>NS-E</td><td>0.488</td><td>0.386</td><td>0.274</td><td>0.701</td><td>3.903</td><td>0.488</td><td>0.386</td><td>0.273</td><td>0.688</td><td>3.782</td></tr><tr><td>S-D</td><td>0.490</td><td>0.389</td><td>0.277</td><td>0.729</td><td>4.091</td><td>0.484</td><td>0.383</td><td>0.271</td><td>0.680</td><td>3.733</td></tr><tr><td>NS-D</td><td><strong>0.485</strong></td><td><strong>0.383</strong></td><td><strong>0.271</strong></td><td><strong>0.683</strong></td><td><strong>3.789</strong></td><td><strong>0.484</strong></td><td><strong>0.383</strong></td><td><strong>0.271</strong></td><td><strong>0.678</strong></td><td><strong>3.723</strong></td></tr></tbody></table>

In conclusion, this simulation study shows that WPL is a valid procedure for estimating the parameters of the Gneiting-Matérn multivariate space-time model when analyzing large data sets for which a Full Likelihood maximization is not feasible. The fact that the difference in scores between the models is more pronounced with WPL than with FL is a strong indication that it is particularly interesting to apply WPL to the non-separable NS-D model. When selecting the correct model, we note that in terms of prediction performances (RMSE, MAE, CRPS), the differences are negligible when using WPL estimates instead of FL estimates despite the lower estimating efficiency shown in Table 4.

# 5 Western France weather dataset

This section illustrates the use of the Gneiting-Matérn multivariate space-time model for the analysis of a weather dataset which consists of three daily variables (solar radiation, R, temperature, T, and humidity, H) recorded at 13 stations in Western France from 2003 to 2012 (see Figure 4). These data are part of the French National Institute for Agricultural Research (INRA) archive of weather data. The considered domain experiences an oceanic climate, characterized by moist and cool (but not cold) winters. Due to the prevailing westerly winds, some amount of space-time interaction is thus expected, from West to East. The validation stations (indicated with stars in Figure 4) have been selected in the Eastern part of the domain, neither too close, nor too far from the other stations. The first validation
<figure><img src="image_3.png" /><figcaption>Figure 3: Prediction scores RMSE, MAE, CRPS and $LogS_1$ for the models S-E, NS-E, S-D and NS-D. Same simulations as in Table 1. Estimates obtained with WPL.</figcaption></figure>
station, Le Rheu, is located near Rennes, in Brittany. The second one, Bourran, is located near Agen, in Aquitaine. The other 11 stations are used to estimate the parameters. In order to restrict ourselves to data being stationary in time, we selected data recorded in January, from 2003 to 2012.

<figure><img src="image_4.png" /><figcaption>Figure 4: Location of the 13 weather stations over western France. Stars: validation stations.</figcaption></figure>

To model the January data $\mathbf{Y}(\mathbf{s}, t) = \{Y_R(\mathbf{s}, t), Y_T(\mathbf{s}, t), Y_H(\mathbf{s}, t)\}^\top$, we consider a model that is stationary in time and non-stationary in space, namely

$$ Y_i(\mathbf{s}, t) = \mu_i(\mathbf{s}) + \sigma_i(\mathbf{s}) Z_i(\mathbf{s}, t) \quad i \in \{\text{R}, \text{T}, \text{H}\} \quad (\mathbf{s}, t) \in \mathbb{R}^2 \times \mathbb{R}. \qquad (13) $$

The data are first standardized at each location by their averages and standard deviations. The space-time modeling will be carried out on the standardized variables,

$$ \mathbf{Z}(\mathbf{s}, t) = \{Z_R(\mathbf{s}, t), Z_T(\mathbf{s}, t), Z_H(\mathbf{s}, t)\}^\top, $$

which are assumed to be space-time stationary and multivariate Gaussian. The simultaneous observation of the spatial and temporal variograms of the daily temperature standardized variables, depicted in Figure 5, reveals a zonal anisotropy between space and time. While
<figure><img src="image_5.png" /><figcaption>Figure 5: Marginal spatial and temporal variograms for daily January temperature standardized variables.</figcaption></figure>

the temporal variogram reaches a sill equal to the overall variance, equal to 1 due to the standardization, the spatial variogram does not exceed 0.4, even at very large distances. This zonal anisotropy is modeled by considering that the process Z arises as the sum of a zero mean multivariate temporal Gaussian process X accounting for temporal effects being constant in space, and an independent zero mean multivariate space-time Gaussian process W

$$Z_i(s, t) = X_i(t) + W_i(s, t) \quad i \in \{\mathbb{R}, T, H\}, \quad (s, t) \in \mathbb{R}^2 \times \mathbb{R}. \qquad (14)$$

Physically speaking, the model (14) assumes that part of the daily changes in the weather pattern, which accounts for 60% of the overall variability, is regional and affects the whole domain under study. The stationary space-time variations account for 40% of the variability.

This model provides the flexibility to accurately model the zonal anisotropy observed in Figure 5, since

$$\begin{align*} \gamma_{ii}^{(\mathbf{Z})}(\mathbf{h}, 0) &= \gamma_{ii}^{(\mathbf{X})}(0) + \gamma_{ii}^{(\mathbf{W})}(\mathbf{h}, 0) \underset{\mathbf{h} \to \infty}{\longrightarrow} [\sigma_i^{(\mathbf{W})}]^2 \\ \gamma_{ii}^{(\mathbf{Z})}(\mathbf{0}, u) &= \gamma_{ii}^{(\mathbf{X})}(u) + \gamma_{ii}^{(\mathbf{W})}(\mathbf{0}, u) \underset{u \to \infty}{\longrightarrow} [\sigma_i^{(\mathbf{X})}]^2 + [\sigma_i^{(\mathbf{W})}]^2 = 1. \end{align*}$$

We further suppose that the process W is Gaussian, with a Gneiting-Matérn covariance function, as described in Eq. (6) and that the temporal process X has covariance function
$C_{ij}^{(\mathbf{X})}(u) = \sigma_i^{(\mathbf{X})} \sigma_j^{(\mathbf{X})} \beta_{ij}^{(\mathbf{X})} / (\alpha^{(\mathbf{X})}|u|^{2a^{(\mathbf{X})}} + 1)$, with $i, j \in \{\mathrm{R}, \mathrm{T}, \mathrm{H}\}$ and $u \in \mathbb{R}$. Note that this temporal covariance function is the temporal marginal of a Gneiting-Matérn covariance (i.e. when $\mathbf{h} = \mathbf{0}$).

Regarding the process $\mathbf{W}$, we consider the same four space-time models as in the previous section, for which there are at most 15 parameters to be estimated (for model NS-D). There are 8 parameters for the process $\mathbf{X}$. Since the data have been standardized, we impose $[\sigma_i^{(\mathbf{X})}]^2 + [\sigma_i^{(\mathbf{W})}]^2 = 1$, $i \in \{\mathrm{R}, \mathrm{T}, \mathrm{H}\}$, thus leading to a total of at most 20 parameters. Estimates are reported in Table 6. We also find that the model NS-E performs worst than S-E. Figures 6 and 7 and represent the empirical spatial direct and cross-covariance functions and the model NS-D with estimated parameters using WPL. We observe that, in general, the model fits quite well the experimental covariance. In particular, the model is able to fit the different regularities, as well as the space-time interactions. It is slightly better fitted for the spatial dimensions than along time, probably because of much less flexibility in the model along time. Conditional expectation and envelopes of 100 conditional simulations are shown in Figure 8 at the two validation locations. The true values are always contained within the envelopes.

The analysis of the colocated covariance coefficients between the three variables reveals an interesting pattern. We compute the estimated covariance coefficient $\sigma_{ij}^{(\mathbf{X})} = \sigma_i^{(\mathbf{X})} \sigma_j^{(\mathbf{X})} \rho_{ij}^{(\mathbf{X})}$, with $i, j \in \{\mathrm{R}, \mathrm{T}, \mathrm{H}\}$, where $\rho_{ij}^{(\mathbf{X})}$ is the correlation coefficient as defined in Theorem 1. Similar coefficients are computed for $\mathbf{W}$ and $\mathbf{Z}$ and compared to the empirical correlation coefficients in Table 8. The correlation coefficients estimated for $\mathbf{Z}$ are very close to empirical ones. Note that Radiation is always negatively correlated to Temperature and Humidity, which is typical for winter weather conditions. More interestingly, we can observe that the (negative) correlation between Radiation and Humidity is mostly due to the space-time process $\mathbf{W}$, whereas the (positive) correlation between Temperature and Humidity is mostly due to temporal process $\mathbf{X}$.
<table><caption>Table 6: Estimates of the parameters using WPL with $d = (500\text{km}, 2\text{days})$. Unit is kms for distances.</caption><thead><tr><th></th><th>σ<sub>R</sub><sup>(X)</sup></th><th>σ<sub>T</sub><sup>(X)</sup></th><th>σ<sub>H</sub><sup>(X)</sup></th><th>β<sub>RT</sub><sup>(X)</sup></th><th>β<sub>RH</sub><sup>(X)</sup></th><th>β<sub>TH</sub><sup>(X)</sup></th><th>σ<sub>R</sub><sup>(W)</sup></th><th>σ<sub>T</sub><sup>(W)</sup></th><th>σ<sub>H</sub><sup>(W)</sup></th><th>β<sub>RT</sub><sup>(W)</sup></th><th>β<sub>RH</sub><sup>(W)</sup></th><th>β<sub>TH</sub><sup>(W)</sup></th></tr></thead><tbody><tr><td>S-E</td><td>0.38</td><td>0.90</td><td>0.49</td><td>-0.72</td><td>-0.52</td><td>0.66</td><td>0.93</td><td>0.43</td><td>0.87</td><td>-0.46</td><td>-0.45</td><td>-0.02</td></tr><tr><td>NS-E</td><td>0.33</td><td>0.90</td><td>0.45</td><td>-0.79</td><td>-0.52</td><td>0.71</td><td>0.94</td><td>0.44</td><td>0.89</td><td>-0.47</td><td>-0.45</td><td>-0.02</td></tr><tr><td>S-D</td><td>0.20</td><td>0.85</td><td>0.43</td><td>-0.98</td><td>-0.56</td><td>0.72</td><td>0.98</td><td>0.53</td><td>0.90</td><td>-0.53</td><td>-0.44</td><td>0.08</td></tr><tr><td>NS-D</td><td>0.19</td><td>0.87</td><td>0.44</td><td>-0.97</td><td>-0.55</td><td>0.73</td><td>0.98</td><td>0.50</td><td>0.90</td><td>-0.59</td><td>-0.47</td><td>0.06</td></tr></tbody></table><table><thead><tr><th></th><th>ν<sub>R</sub></th><th>ν<sub>T</sub></th><th>ν<sub>H</sub></th><th>1/r<sub>R</sub></th><th>1/r<sub>T</sub></th><th>1/r<sub>H</sub></th><th>α<sup>(X)</sup></th><th>a<sup>(X)</sup></th><th>a<sup>(W)</sup></th><th>b</th><th>-wcl</th></tr></thead><tbody><tr><td>S-E</td><td>0.66</td><td>0.66</td><td>0.66</td><td>184</td><td>184</td><td>184</td><td>0.12</td><td>1.00</td><td>2.35</td><td>0.89</td><td>0</td><td>631995</td></tr><tr><td>NS-E</td><td>0.67</td><td>0.67</td><td>0.67</td><td>188</td><td>188</td><td>188</td><td>0.13</td><td>0.93</td><td>2.08</td><td>1.00</td><td>0.6</td><td>631991</td></tr><tr><td>S-D</td><td>0.75</td><td>0.84</td><td>0.41</td><td>211</td><td>249</td><td>248</td><td>0.09</td><td>1.00</td><td>0.91</td><td>0.69</td><td>0</td><td>631922</td></tr><tr><td>NS-D</td><td>0.64</td><td>0.84</td><td>0.30</td><td>239</td><td>224</td><td>322</td><td>0.12</td><td>1.00</td><td>1.06</td><td>0.75</td><td>0.6</td><td><b>631845</b></td></tr></tbody></table>

<table><caption>Table 7: RMSE, MAE, CRPS and LogS₁ of predicted values at the validation stations, using WPL estimates with $d = (500\text{km}, 2\text{days})$.</caption><thead><tr><th></th><th>RMSE</th><th>MAE</th><th>CRPS</th><th>LogS<sub>1</sub></th></tr></thead><tbody><tr><td>S-E</td><td>0.433</td><td>0.306</td><td>0.222</td><td>0.440</td></tr><tr><td>NS-E</td><td>0.446</td><td>0.313</td><td>0.228</td><td>0.467</td></tr><tr><td>S-D</td><td>0.417</td><td>0.295</td><td><b>0.215</b></td><td>0.405</td></tr><tr><td>NS-D</td><td><b>0.417</b></td><td><b>0.295</b></td><td>0.216</td><td><b>0.397</b></td></tr></tbody></table>

<table><caption>Table 8: Empirical and estimated covariance coefficients $\sigma_{ij} = \sigma_i\sigma_j\rho_{ij}$, with $i, j \in \{\text{R, T, H}\}$ for the model NS-D.</caption><thead><tr><th></th><th>σ<sub>RT</sub></th><th>σ<sub>RH</sub></th><th>σ<sub>TH</sub></th></tr></thead><tbody><tr><td>Process <b>X</b></td><td>-0.16</td><td>-0.05</td><td>0.28</td></tr><tr><td>Process <b>W</b></td><td>-0.29</td><td>-0.39</td><td>0.03</td></tr><tr><td>Process <b>Z</b></td><td>-0.45</td><td>-0.44</td><td>0.31</td></tr><tr><td>Empirical</td><td>-0.38</td><td>-0.42</td><td>0.25</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 6: Spatial direct and cross-covariance functions for R, T and H. Points: empirical values. Solid line: model NS-D with parameters estimated using WPL (see Table 6). Left panel: $u = 0$. Right panel: $u = 1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Temporal direct and cross-covariance functions for R, T and H, as in Figure 6. Left panel: $\mathbf{h} = 0$. Right panel: $||\mathbf{h}|| = 199$ km.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 8: Prediction of Radiation, Temperature and Humidity at the two validation stations, Le Rheu (Brittany) and Bourran (Aquitaine). Black points: true values. Dotted lines: conditional expectation. Shaded area: 90% envelope interval, computed on 100 simulations.</figcaption></figure>
# 6 Discussion

Flexible multivariate space-time models can be built using a mixture approach. We proposed two classes based on the univariate Gneiting class of space-time covariances: the Gneiting-Matérn and the Gneiting-Cauchy multivariate space-time models. A simulation study showed that the relatively large number of parameters could be accurately estimated by maximizing a weighted pairwise likelihood function. Whenever data are simulated from the most general model, the analysis of the estimation bias and variance, as well as that of validation scores revealed that ignoring space-time interaction or the multivariate flexibility consistently leads to non optimal scores. The relative efficiency of estimates obtained when maximizing weighted pairwise likelihood with respect to those obtained when maximizing a full likelihood is, in general, above 70%. It thereby confirms that weighted pairwise likelihood is a reliable alternative to full likelihood when analyzing large datasets.

Our approach could be generalized in several ways. First, with similar arguments to those used in Theorems 1 and 2, we could propose a wider range of flexible multivariate space-time models, including compactly supported models, and models on manifolds such as spheres. One limit of our models is that there is a unique temporal covariance for all variables. By exchanging the role of time and space, one can easily build models that are flexible in time and unique in space. Building valid multivariate space-time models that are flexible both in time and space is still an open challenge. For the time being, the analysis of models requiring different temporal covariances for different variables can be carried out using a sum of models, as in Section 5.

WPL has been shown to provide accurate estimates for univariate spatial and space-time Gaussian random fields (Bevilacqua et al., 2012; Bevilacqua and Gaetan, 2015). We showed that it also provides accurate estimates for multivariate space-time Gaussian random fields. We found, however, that the space-time separability parameter, *b*, was particularly difficult
to estimate. The effect of the separability, which could be formally defined as the ratio between the non-separable space-covariance and the separable one, decreases as $||h|| \to 0$ or $u \to 0$. Therefore, the relative efficiency of the estimator of $b$ decreases when the window $d = (d_S, d_T)$ defining the pairs in WPL is too small. The optimal size of this window is thus a balance between two opposite requirements: accurately estimating $b$, while offering speed gain. Our opinion is that one could probably increase the efficiency of the estimation of the parameters without sacrificing too much computation speed by proposing new forms of composite likelihoods.

# Appendix

Our proofs are based on the scale mixture representation (5), for which, according to Lemma 1, we need the following ingredients: a relevant univariate covariance function $C_\xi(\mathbf{h}, u)$ in $\mathbb{R}^k \times \mathbb{R}$ and a non-negative definite multivariate mixture $\mathbf{m}(\xi) = [m_{ij}(\xi)]_{i,j=1}^p, \xi > 0$. The measure $F$ in Lemma 1 will be set to be the Lebesgue measure on $(0, \infty)$.

The first part of the proof is common to both Theorems. We first verify that the function

$$C_{\xi}(\mathbf{h}, u) = \frac{1}{\psi(u^2)^{k/2}} \exp \left\{ -\xi \left( \frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}} \right)^{\lambda} \right\}, \quad (15)$$

defined on $\mathbb{R}^k \times \mathbb{R}$ with $0 < \lambda \le 2$, is a valid covariance function for any $\xi > 0$ and for $\psi(t), t \ge 0$, being a positive function with a completely monotone derivative. Indeed, the mapping $t \mapsto \exp(-ct^\gamma)$, $t \ge 0$, $c > 0$ and $0 < \gamma \le 1$ is completely monotone. Then, setting $\gamma = 1$, $c = (\frac{\|h\|}{\psi(u^2)^{1/2}})^\lambda$ and $t = \xi$, it is clear that $C_\xi(\mathbf{h}, u)$ belongs to the Gneiting class of covariance functions (Gneiting, 2002) for any $\xi > 0$. Also, the mapping $[0, \infty) \to \mathbb{R}: \xi \mapsto C_\xi(\mathbf{h}, u)$ satisfies the hypothesis of Lemma 1, for all fixed $(\mathbf{h}, u) \in \mathbb{R}^k \times \mathbb{R}$.

The rest of the proofs is split into two parts corresponding to each model. It consists in finding multivariate mixtures $\mathbf{m}(\xi) = [m_{ij}(\xi)]_{i,j=1}^p \in \mathbb{R}^p \times \mathbb{R}^p$ tailored to either the Matérn or the Cauchy class, that are symmetric and non-negative definite for all $\xi \ge 0$.
## Proof of Theorem 1

Let $\nu > 0$ and $r > 0$. From Eq. 3.471.9 in Gradshteyn and Ryzhik (2007), we have

$$ \mathcal{M}(\mathbf{h}; r, \nu) = \int_{0}^{\infty} \exp\{-\|\mathbf{h}\|^2 \xi\} m^{\mathcal{M}}(\xi; r, \nu) d\xi, $$

with

$$ m^{\mathcal{M}}(\xi; r, \nu) = \left(\frac{r^2}{4}\right)^{\nu} \frac{\xi^{-1-\nu}}{\Gamma(\nu)} \exp\left\{-\frac{r^2}{4\xi}\right\}, \quad \xi > 0. $$

Following the univariate mixture representation above, we set

$$ \mathbf{m}^{\mathcal{M}}(\xi) = [m_{ij}^{\mathcal{M}}(\xi)]_{i,j=1}^{p} = \left[ \rho_{ij} \left( \frac{r_{ij}^2}{4} \right)^{\nu_{ij}} \frac{\xi^{-1-\nu_{ij}}}{\Gamma(\nu_{ij})} \exp \left\{ -\frac{r_{ij}^2}{4\xi} \right\} \right]_{i,j=1}^{p}, \quad \xi > 0, $$

with

$$ r_{ij} = \left\{ \frac{(r_i^2 + r_j^2)}{2} \right\}^{1/2}, \quad \forall i, j = 1, \dots, p \text{ and } r_i > 0, \quad \forall i = 1, \dots, p, $$

$$ \nu_{ij} = (\nu_i + \nu_j)/2, \quad \forall i, j = 1, \dots, p, $$

$$ \rho_{ij} = \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2} \Gamma(\nu_j)^{1/2}} \frac{r_i^{\nu_i} r_j^{\nu_j}}{r_{ij}^{2\nu_{ij}}}, \quad \forall i, j = 1, \dots, p, $$

where $\beta = [\beta_{ij}]_{i,j=1}^p$ is a correlation matrix. We have, for each $i, j = 1, \dots, p$,

$$ \begin{aligned} m_{ij}^{\mathcal{M}}(\xi) &= \rho_{ij} \left(\frac{r_{ij}^2}{4}\right)^{\nu_{ij}} \frac{\xi^{-1-\nu_{ij}}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{r_{ij}^2}{4\xi}\right\} \\ &= \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2} \Gamma(\nu_j)^{1/2}} \frac{r_i^{\nu_i} r_j^{\nu_j}}{r_{ij}^{2\nu_{ij}}} \left(\frac{r_{ij}^2}{4}\right)^{\nu_{ij}} \frac{\xi^{-1-\nu_{ij}}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{r_{ij}^2}{4\xi}\right\} \\ &= \beta_{ij} \frac{r_i^{\nu_i}}{2^{\nu_i} \Gamma(\nu_i)^{1/2}} \xi^{-(1+\nu_i)/2} \exp\left\{-\frac{r_i^2}{8\xi}\right\} \cdot \frac{r_j^{\nu_j}}{2^{\nu_j} \Gamma(\nu_j)^{1/2}} \xi^{-(1+\nu_j)/2} \exp\left\{-\frac{r_j^2}{8\xi}\right\} \\ &= \beta_{ij} m_i^{\mathcal{M}}(\xi)^{1/2} m_j^{\mathcal{M}}(\xi)^{1/2}. \end{aligned} $$

Therefore $\mathbf{m}^{\mathcal{M}}(\xi)$ is the Hadamard product of the correlation matrix $\beta$ and the outer product $\bar{\mathbf{m}}^{\mathcal{M}}(\xi)$ by itself, with $\bar{\mathbf{m}}^{\mathcal{M}}(\xi) = (m_1^{\mathcal{M}}(\xi)^{1/2}, \dots, m_p^{\mathcal{M}}(\xi)^{1/2})^\top$. By Schur's Theorem (Horn and Johnson, 2012, p. 455), the matrix $\mathbf{m}^{\mathcal{M}}(\xi)$ is thus non-negative definite since it is the Hadamard product of two non-negative definite matrices.
Setting $\lambda = 2$ in Eq. (15) with $\mathbf{m}^M(\xi) = [m_{ij}^M(\xi)]_{i,j=1}^p$ as defined above leads thus to the valid $p$-variate matrix-valued covariance function $[C_{ij}^M(\mathbf{h}, u)]_{i,j=1}^p$ on $\mathbb{R}^k \times \mathbb{R}$ by application of Lemma 1:

$$ 
\begin{aligned}
C_{ij}^{M}(\mathbf{h}, u) &= \sigma_i \sigma_j \int_0^\infty \frac{1}{\psi(\|\mathbf{u}\|^2)^{k/2}} \exp\left\{-\xi \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^2\right\} m_{ij}^M(\xi) d\xi \\
&= \frac{\sigma_i \sigma_j}{\psi(u^2)^{k/2}} \rho_{ij} \mathcal{M}\left(\frac{\mathbf{h}}{\psi(\|\mathbf{u}\|^2)^{1/2}}; r_{ij}, \nu_{ij}\right)
\end{aligned}
$$

<figure><img src="image_9.png" /></figure>

## Proof of Theorem 2

It is well known that the Cauchy covariance function is the Laplace transform of a Gamma distribution. Therefore, we set

$$ \mathbf{m}^C(\xi) = [m_{ij}^C(\xi)]_{i,j=1}^p = \left[ \rho_{ij} \frac{1}{r_{ij}^{\nu_{ij}}} \frac{\xi^{\nu_{ij}-1}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{1}{r_{ij}}\xi\right\} \right]_{i,j=1}^p, \quad \xi > 0, $$

with

$$ r_{ij} = \{(r_i^{-1} + r_j^{-1})/2\}^{-1}, \forall i, j = 1, \dots, p \text{ and } r_i > 0, \forall i = 1, \dots, p, $$

$$ \nu_{ij} = (\nu_i + \nu_j)/2, \forall i, j = 1, \dots, p, $$

$$ \rho_{ij} = \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2} \Gamma(\nu_j)^{1/2}} \frac{r_{ij}^{\nu_{ij}}}{(r_i^{\nu_i} r_j^{\nu_j})^{1/2}}, \forall i, j = 1, \dots, p, $$

where $\beta = [\beta_{ij}]_{i,j=1}^p$ is a correlation matrix. This multivariate mixture is non-negative definite. Indeed, we have

$$ 
\begin{aligned}
m_{ij}^{C}(\xi) &= \rho_{ij} \frac{1}{r_{ij}^{\nu_{ij}}} \frac{\xi^{\nu_{ij}-1}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{1}{r_{ij}}\xi\right\} \\
&= \beta_{ij} \frac{\Gamma(\nu_{ij})}{\Gamma(\nu_i)^{1/2} \Gamma(\nu_j)^{1/2}} \frac{r_{ij}^{\nu_{ij}}}{(r_i^{\nu_i} r_j^{\nu_j})^{1/2}} \frac{1}{r_{ij}^{\nu_{ij}}} \frac{\xi^{\nu_{ij}-1}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{1}{r_{ij}}\xi\right\} \\
&= \beta_{ij} \frac{\xi^{\frac{\nu_i-1}{2}}}{(\Gamma(\nu_i)r_i^{\nu_i})^{1/2}} \exp\left\{-\frac{1}{2r_i}\xi\right\} \cdot \frac{\xi^{\frac{\nu_j-1}{2}}}{(\Gamma(\nu_j)r_j^{\nu_j})^{1/2}} \exp\left\{-\frac{1}{2r_j}\xi\right\} \\
&= \beta_{ij} m_i^C(\xi)^{1/2} m_j^C(\xi)^{1/2}.
\end{aligned}
$$

The same arguments as in Theorem 1 are used to conclude that $\mathbf{m}^C(\xi)$ is non-negative definite. The last step of the proof consists in applying Lemma 1 with $\mathbf{m}^C(\xi) = [m_{ij}^C(\xi)]_{i,j=1}^p$
as defined above and by setting $\lambda = 2$ in Eq. (15). The result is the valid $p$-variate matrix-valued covariance function on $\mathbb{R}^k \times \mathbb{R}$

$$ 
\begin{align*}
C_{ij}^{\mathcal{C}}(\mathbf{h}, u) &= \sigma_i \sigma_j \int_0^{\infty} \frac{1}{\psi(u^2)^{k/2}} \exp\left\{-\xi \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^{\lambda}\right\} m_{ij}^{\mathcal{C}}(\xi) d\xi \\
&= \frac{\sigma_i \sigma_j}{\psi(u^2)^{k/2}} \rho_{ij} \int_0^{\infty} \exp\left\{-\xi \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^{\lambda}\right\} \frac{1}{r_{ij}^{\nu_{ij}}} \frac{\xi^{\nu_{ij}-1}}{\Gamma(\nu_{ij})} \exp\left\{-\frac{1}{r_{ij}}\xi\right\} d\xi \\
&= \frac{\sigma_i \sigma_j}{\psi(u^2)^{k/2}} \rho_{ij} \left\{1 + r_{ij} \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^{\lambda}\right\}^{-\nu_{ij}} \\
&\quad \times \int_0^{\infty} \left\{\frac{1}{r_{ij}} + \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^{\lambda}\right\}^{\nu_{ij}} \frac{\xi^{\nu_{ij}-1}}{\Gamma(\nu_{ij})} \exp\left\{-\xi \left[\frac{1}{r_{ij}} + \left(\frac{\|\mathbf{h}\|}{\psi(u^2)^{1/2}}\right)^{\lambda}\right]\right\} d\xi \\
&= \frac{\sigma_i \sigma_j}{\psi(u^2)^{k/2}} \rho_{ij} C\left(\frac{\mathbf{h}}{\psi(u^2)^{1/2}; r_{ij}, \nu_{ii}, \lambda}\right) \square
\end{align*}
$$

## Acknowledgments

Two anonymous referees are gratefully acknowledged. The authors wish to thank Carlo Gaetan from the Ca' Foscari University in Venezia (Italy) for fruitful discussions and Liliane Bel from AgroParisTech in Paris (France) for her very careful reading of earlier versions of this manuscript. This work was supported by the metaprogramme Adaptation of Agriculture and Forests to Climate Change (AAFCC) of the French National Institute for Agricultural Research (INRA). Emilio Porcu has been supported by Proyecto Fondecyt n. 1130647 from Chilean Ministry of Education.

## References

Apanasovich, T. V. and Genton, M. G. (2010). Cross-covariance functions for multivariate random fields based on latent dimensions. Biometrika, 97(1):15–30.

Apanasovich, T. V., Genton, M. G., and Sun, Y. (2012). A valid Matérn class of cross-
covariance functions for multivariate random fields with any number of components. *Journal of the American Statistical Association*, 107(497):180-193.

Bevilacqua, M. and Gaetan, C. (2015). Comparing composite likelihood methods based on pairs for spatial Gaussian random fields. *Statistics and Computing*, 25(5):877-892.

Bevilacqua, M., Gaetan, C., Mateu, J., and Porcu, E. (2012). Estimating space and space-time covariance functions for large data sets: a weighted composite likelihood approach. *Journal of the American Statistical Association*, 107(497):268-280.

Castruccio, S., Huser, R., and Genton, M. (2015). High-order composite likelihood inference for max-stable distributions and processes. arXiv:1411.0086v3.

Cressie, N. (1993). *Statistics for Spatial Data: Wiley Series in Probability and Statistics*. Wiley-Interscience New York.

Cressie, N. and Wikle, C. K. (2011). *Statistics for spatio-temporal data*. John Wiley & Sons.

Daley, D., Porcu, E., and Bevilacqua, M. (2014). Classes of compactly supported covariance functions for multivariate random fields. *Stochastic Environmental Research and Risk Assessment*, 29(4):1249-1263.

De Iaco, S., Myers, D., Palma, M., and Posa, D. (2013). Using simultaneous diagonalization to identify a space-time linear coregionalization model. *Mathematical Geosciences*, 45(1):69-86.

DLMF. NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.10 of 2015-08-07.

Gelfand, A. E. and Banerjee, S. (2010). Multivariate spatial process models. In Gelfand, A. E., Diggle, P., Fuentes, M., and Guttorp, P., editors, *Handbook of Spatial Statistics*, pages 495-515. Chapman & Hall/CRC.
Genton, M. G. and Kleiber, W. (2015). Cross-covariance functions for multivariate geostatistics. *Statistical Science*, 30(2):147–163.

Gneiting, T. (2002). Nonseparable, stationary covariance functions for space-time data. *Journal of the American Statistical Association*, 97(458):590–600.

Gneiting, T., Genton, M., and Guttorp, P. (2007). Geostatistical space-time models, stationarity, separability and full symmetry. In Finkenstädt, B., Held, L., and Isham, V., editors, *Statistical Methods for Spatio-Temporal Systems*, pages 151–175. Chapman & Hall/CRC.

Gneiting, T., Kleiber, W., and Schlather, M. (2010). Matérn cross-covariance functions for multivariate random fields. *Journal of the American Statistical Association*, 105(491):1167–1177.

Gneiting, T. and Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and estimation. *Journal of the American Statistical Association*, 102(477):359–378.

Gneiting, T. and Schlather, M. (2004). Stochastic models that separate fractal dimension and the Hurst effect. *SIAM Review*, 46(2):269–282.

Goulard, M. and Voltz, M. (1992). Linear coregionalization model: Tools for estimation and choice of cross-variogram matrix. *Mathematical Geology*, 24(3):269–286.

Gradshteyn, I. and Ryzhik, I. (2007). *Table of integrals, series, and products, Seventh edition*. Academic Press.

Horn, R. A. and Johnson, C. R. (2012). *Matrix analysis*. Cambridge University Press.

Li, B., Genton, M. G., and Sherman, M. (2008). Testing the covariance structure of multivariate random fields. *Biometrika*, 95(4):813–829.

Lindsay, B. G. (1988). Composite likelihood methods. *Contemporary Mathematics*, 80(1):221–239.
Matérn, B. (1986). *Spatial Variation*. Lecture Notes in Statistics, 36, Springer.

Porcu, E. and Schilling, R. L. (2011). From Schoenberg to Pick-Nevanlinna: Toward a complete picture of the variogram class. *Bernoulli*, 17(1):441–455.

Porcu, E. and Zastavnyi, V. (2011). Characterization theorems for some classes of covariance functions associated to vector valued random fields. *Journal of Multivariate Analysis*, 102:1293–1301.

Reisert, M. and Burkhardt, H. (2007). Learning equivariant functions with matrix valued kernels. *The Journal of Machine Learning Research*, 8:385–408.

Schlather, M. (2010). Some covariance models based on normal scale mixtures. *Bernoulli*, 16(3):780–797.

Schoenberg, I. J. (1938). Metric spaces and completely monotone functions. *Annals of Mathematics*, 39(4):811–841.

Varin, C., Reid, N., and Firth, D. (2011). An overview of composite likelihood methods. *Statistica Sinica*, 21(1):5–42.

Wackernagel, H. (2003). *Multivariate geostatistics, Third Edition*. Springer, Berlin.

Zhang, H. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. *Journal of the American Statistical Association*, 99(465):250–261.