Covariance Functions for Multivariate Gaussian Fields Evolving Temporally over Planet Earth

TL;DR

This paper develops new flexible covariance functions for multivariate Gaussian fields evolving over time on spherical surfaces, crucial for modeling climate and ocean data across the globe.

Contribution

It introduces a novel family of matrix-valued covariance functions for spherical space-time data, extending Gneiting models and latent dimension approaches.

Findings

Models effectively capture space-time dependencies on the sphere.

Proposed covariance functions improve prediction accuracy.

Application to climate data demonstrates practical utility.

Abstract

The construction of valid and flexible cross-covariance functions is a fundamental task for modeling multivariate space-time data arising from climatological and oceanographical phenomena. Indeed, a suitable specification of the covariance structure allows to capture both the space-time dependencies between the observations and the development of accurate predictions. For data observed over large portions of planet Earth it is necessary to take into account the curvature of the planet. Hence the need for random field models defined over spheres across time. In particular, the associated covariance function should depend on the geodesic distance, which is the most natural metric over the spherical surface. In this work, we propose a flexible parametric family of matrix-valued covariance functions, with both marginal and cross structure being of the Gneiting type. We additionally introduce a different multivariate Gneiting model based on the adaptation of the latent dimension approach to the spherical context. Finally, we assess the performance of our models through the study of a bivariate space-time data set of surface air temperatures and precipitations.

Figures

Tables

Table 1. Table 3.1 : Some examples of completely monotone functions. Here, K ν subscript 𝐾 𝜈 K_{\nu} denotes the modified Bessel function of second kind of degree ν 𝜈 \nu .

Function	Parameters restriction
$g\left(t\right)=\exp \left(-c{t}^{\gamma }\right)$	$c>0$,
$g\left(t\right)={\left(2^{\nu -1}\Gamma \left(\nu \right)\right)}^{-1}{\left(c\sqrt{t}\right)}^{\nu }{K}_{\nu }\left(c\sqrt{t}\right)$	$c>0$, $\nu >0$
$g\left(t\right)={\left(1+c{t}^{\gamma }\right)}^{-\nu }$	$c>0$, , $\nu >0$
$g\left(t\right)=2^{\nu }{\left(\exp \left(c\sqrt{t}\right)+\exp \left(-c\sqrt{t}\right)\right)}^{-\nu }$	$c>0$, $\nu >0$

Table 2. Table 3.2 : Some examples of Bernstein functions.

Function	Parameters restriction
$f\left(t\right)={\left(a{t}^{\alpha }+1\right)}^{\beta }$	$a>0$, , $0\le \beta \le 1$
$f\left(t\right)=\ln \left(a{t}^{\alpha }+b\right)/\ln \left(b\right)$	$a>0$, $b>1$
$f\left(t\right)=\left(a{t}^{\alpha }+b\right)/\left(b\left(a{t}^{\alpha }+1\right)\right)$	$a>0$, ,

Table 3. Table 4.1 : CL estimates for Models A-C. For the separable Models A and B, there is a single spatial scale parameter: c S := c 11 = c 22 = c 12 assign subscript 𝑐 𝑆 subscript 𝑐 11 subscript 𝑐 22 subscript 𝑐 12 c_{S}:=c_{11}=c_{22}=c_{12} . For Model C, the cross scale parameter c 12 subscript 𝑐 12 c_{12} is the maximum between the marginal spatial scales c 11 subscript 𝑐 11 c_{11} and c 22 subscript 𝑐 22 c_{22} .

		${\hat{\sigma }}_{11}^2$	${\hat{\sigma }}_{22}^2$	${\hat{\rho }}_{12}$	${\hat{c}}_{11}$	${\hat{c}}_{22}$	${\hat{c}}_{12}$	${\hat{c}}_T$
Model A		1.85	$4.71\times {10}^{-5}$	0.28	102	102	102	11.88
Model B		1.84	$4.72\times {10}^{-5}$	0.28	2602	2602	2602	22.58
Model C		1.85	$4.69\times {10}^{-5}$	0.28	2901	2245	2901	22.92

Table 4. Table 4.2 : CL estimates for Model D.

		${\hat{a}}_{11}$	${\hat{a}}_{21}$	${\hat{a}}_{22}$	${\hat{c}}_{S,1}$	${\hat{c}}_{S,2}$	${\hat{c}}_{T,1}$	${\hat{c}}_{T,2}$	${\hat{\zeta }}_{12}$
Model D		1.38	$3.98\times {10}^{-3}$	$5.58\times {10}^{-3}$	47133	30805	41.36	4.03	1.66

Table 5. Table 4.3 : Comparison between Models A-D in terms of Log-CL and cross-validation scores. Here, MSE i subscript MSE 𝑖 {\rm MSE}_{i} and CRPS i subscript CRPS 𝑖 {\rm CRPS}_{i} denote respectively the MSE and CRPS associated to Variable i 𝑖 i , for i = 1 , 2 𝑖 1 2 i=1,2 .

	Log-CL $\times {10}^{-5}$	${\mathrm{MSE}}_1$	${\mathrm{MSE}}_2\times {10}^5$	${\mathrm{CRPS}}_1$	${\mathrm{CRPS}}_2$
Model A	5.79	0.27	2.07	1.21	1.52
Model B	5.80	0.13	1.62	0.67	1.02
Model C	5.81	0.11	1.61	0.62	0.93
Model D	5.81	0.11	1.84	0.63	0.62

Table 6. Table A.1 : Summary of the notation used along the Appendices. Throughout, in the univariate case ( m = 1 𝑚 1 m=1 ) we omit the super index: 𝒫 ​ ( ℰ ) 𝒫 ℰ \mathcal{P}(\mathcal{E}) , Φ d , 𝒮 subscript Φ 𝑑 𝒮 \Phi_{d,\mathcal{S}} , Φ d , ℐ subscript Φ 𝑑 ℐ \Phi_{d,\mathcal{I}} , Ψ d , ℐ subscript Ψ 𝑑 ℐ \Psi_{d,\mathcal{I}} and Υ d , k subscript Υ 𝑑 𝑘 \Upsilon_{d,k} are used instead of 𝒫 1 ​ ( ℰ ) superscript 𝒫 1 ℰ \mathcal{P}^{1}(\mathcal{E}) , Φ d , 𝒮 1 superscript subscript Φ 𝑑 𝒮 1 \Phi_{d,\mathcal{S}}^{1} , Φ d , ℐ 1 superscript subscript Φ 𝑑 ℐ 1 \Phi_{d,\mathcal{I}}^{1} , Ψ d , ℐ 1 superscript subscript Ψ 𝑑 ℐ 1 \Psi_{d,\mathcal{I}}^{1} and Υ d , k 1 superscript subscript Υ 𝑑 𝑘 1 \Upsilon_{d,k}^{1} , respectively.

Notation	Description
${𝒫}^m\left(ℰ\right)$	Positive definite matrix-valued ($m\times m$) functions on $ℰ\times ℰ$.
${\Phi }_{d,𝒮}^m$	Continuous, bounded and stationary elements of ${𝒫}^m\left({ℝ}^d\right)$.
${\Phi }_{d,ℐ}^m$	Continuous, bounded and Euclidean isotropic elements of ${𝒫}^m\left({ℝ}^d\right)$.
${\Psi }_{d,ℐ}^m$	Continuous, bounded and geodesically isotropic elements of ${𝒫}^m\left({𝕊}^d\right)$.
${{\rm Y}}_{d,k}^m$	Elements in ${𝒫}^m\left({𝕊}^d\times {ℝ}^k\right)$ being, continuous, bounded, geodesically
	isotropic in the spherical variable and stationary in the Euclidean variable.

Equations

$$\sum _{\ell =1}^n\sum _{r=1}^n{\mathbf{a}}_{\ell }^{\top }\text{C}(\theta ({\mathbf{x}}_{\ell },{\mathbf{x}}_r),t_{\ell }-t_r){\mathbf{a}}_r\ge 0,$$

$$\text{C}(\theta ,u)={\text{C}}_{\mathtt{S}}(\theta )\circ {\text{C}}_{\mathtt{T}}(u),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$$\text{C}(\theta ,u)=\text{A}{\text{C}}_{\mathtt{S}}(\theta ){\text{C}}_{\mathtt{T}}(u),(\theta ,u)\in [0,\pi ]\times \mathbb{R}.$$

$$\text{C}(\theta ,u)=\text{A}\text{C}(\theta ,u),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$$\text{C}(\theta ,u)=\sum _{k=1}^q{\text{A}}_k{\text{C}}_k(\theta ,u),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$$\mathbf{Y}(\mathbf{x},t):=\mathbf{Z}\left({\mathcal{R}}^t\mathbf{x}\right),(\mathbf{x},t)\in {\mathbb{S}}^d\times \mathbb{R}.$$

$${\text{C}}_{ij}(\theta ,u)={\sigma }_{ii}{\sigma }_{jj}{\rho }_{ij}\text{C}(\theta ,u;{\mathbf{\lambda }}_{ij}),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$${\int }_{\Omega }\text{C}(\theta ,u,\xi )\mu (\mathrm{d}\xi )$$

$$\text{C}(\theta ,u)=\frac{1}{{\left\{f(\theta ){|}_{[0,\pi ]}\right\}}^{1/2}}g\left(\frac{|u{|}^2}{f(\theta ){|}_{[0,\pi ]}}\right),\theta \in [0,\pi ],u\in \mathbb{R},$$

$$\text{C}(\theta ,u)=\frac{1}{f(|u|{)}^{n+2}}g\left(\theta f(|u|)\right),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$$\sum _{i\ne j}|{\rho }_{ij}|(c_{ii}/c_{ij}{)}^{n+1}\le 1,$$

$${\text{C}}_{ij}(\theta ,u)=\frac{{\sigma }_{ii}{\sigma }_{jj}{\rho }_{ij}}{f(|u|{)}^{n+2}}g\left(\frac{\theta f(|u|)}{c_{ij}}\right),(\theta ,u)\in [0,\pi ]\times \mathbb{R},$$

$${\text{C}}_{ij}(\theta ,u)=\frac{{\sigma }_{ii}{\sigma }_{jj}{\rho }_{ij}}{{\left(1+\frac{1.7|u|}{c_T}\right)}^3}\exp \left\{-\frac{3\theta \left(1+\frac{1.7|u|}{c_T}\right)}{c_{ij}}\right\},(\theta ,u)\in [0,\pi ]\times \mathbb{R}.$$

$$Z_i(\mathbf{x},t)=Y(\mathbf{x},t,{\mathbf{\xi }}_i),\mathbf{x}\in {\mathbb{S}}^d,t\in \mathbb{R},{\mathbf{\xi }}_i\in {\mathbb{R}}^k,$$

$${\text{C}}_{ij}(\theta ,u)=\text{K}(\theta ,u,{\mathbf{\xi }}_i-{\mathbf{\xi }}_j),\theta \in [0,\pi ],u\in \mathbb{R}.$$

$$\text{K}(\theta ,\mathbf{u},\mathbf{v})=\frac{1}{\{f_2(\theta ){|}_{[0,\pi ]}{\}}^{l/2}{\left\{f_1\left[\frac{\Vert \mathbf{v}{\Vert }^2}{f_2(\theta ){|}_{[0,\pi ]}}\right]\right\}}^{k/2}}g\left(\frac{\Vert \mathbf{u}{\Vert }^2}{f_1\left[\frac{\Vert \mathbf{v}{\Vert }^2}{f_2(\theta ){|}_{[0,\pi ]}}\right]}\right),(\theta ,\mathbf{u},\mathbf{v})\in [0,\pi ]\times {\mathbb{R}}^l\times {\mathbb{R}}^k,$$

$$\text{K}(\theta ,u,{\zeta }_{ij})=\frac{1}{{\left(\frac{400\theta }{c_{S,1}}+1\right)}^{1/2}{\left(\frac{{\zeta }_{ij}}{\left(\frac{400\theta }{c_{S,1}}+1\right)}+1\right)}^{1/2}}\exp \left\{-\frac{3|u|/c_{T,1}}{{\left(\frac{{\zeta }_{ij}}{\left(\frac{400\theta }{c_{S,1}}+1\right)}+1\right)}^{1/2}}\right\},$$

$$\text{R}(\theta ,u)=\frac{1}{{\left(\frac{400\theta }{c_{S,2}}+1\right)}^{1/2}}\exp \left\{-\frac{3|u|/c_{T,2}}{{\left(\frac{400\theta }{c_{S,2}}+1\right)}^{1/2}}\right\}.$$

$${\text{C}}_{11}(\theta ,u)$$

$${\text{C}}_{22}(\theta ,u)$$

$${\text{C}}_{12}(\theta ,u)$$

$${\text{C}}_{ij}(\theta ,u)=\frac{{\sigma }_{ii}{\sigma }_{jj}{\rho }_{ij}}{{\left(\frac{400\theta }{c_S}+1\right)}^{1/2}}\exp \left\{\frac{-3|u|/c_T}{{\left(\frac{400\theta }{c_S}+1\right)}^{1/2}}\right\},$$

$$\sum _{\ell =1}^n\sum _{r=1}^n{\mathbf{a}}_{\ell }^{\top }\text{F}(e_{\ell },e_r){\mathbf{a}}_r\ge 0.$$

$$\text{F}(\mathbf{x},\mathbf{y})=\tilde{\mathbf{\phi }}(\mathbf{x}-\mathbf{y})=[{\tilde{\phi }}_{ij}(\mathbf{x}-\mathbf{y}){]}_{i,j=1}^m,\mathbf{x},\mathbf{y}\in {\mathbb{R}}^d.$$

$$\tilde{\mathbf{\phi }}(\mathbf{h})={\int }_{{\mathbb{R}}^d}\exp \{ı{\mathbf{h}}^{\top }\mathbf{\omega }\}\text{d}{\tilde{\mathbf{\Lambda }}}_d(\mathbf{\omega }),\mathbf{h}\in {\mathbb{R}}^d,$$

$$\tilde{\mathbf{\phi }}(\mathbf{h})={\int }_{{\mathbb{R}}^d}\exp \{ı{\mathbf{h}}^{\top }\mathbf{\omega }\}{\tilde{\mathbf{\lambda }}}_d(\mathbf{\omega })\text{d}\mathbf{\omega },\mathbf{h}\in {\mathbb{R}}^d,$$

$${\tilde{\mathbf{\lambda }}}_d(\mathbf{\omega })=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}\exp \{-ı{\mathbf{h}}^{\top }\mathbf{\omega }\}\tilde{\mathbf{\phi }}(\mathbf{h})\text{d}\mathbf{h},\mathbf{\omega }\in {\mathbb{R}}^d.$$

$$|{\tilde{\phi }}_{ij}(\mathbf{h}){|}^2\le {\tilde{\phi }}_{ii}(0){\tilde{\phi }}_{jj}(0),\text{ for all }\mathbf{h}\in {\mathbb{R}}^d.$$

$$\text{F}(\mathbf{x},\mathbf{y})=\mathbf{\phi }(\Vert \mathbf{x}-\mathbf{y}\Vert ),\mathbf{x},\mathbf{y}\in {\mathbb{R}}^d.$$

$$\mathbf{\phi }(r)={\int }_{[0,\infty )}{\Omega }_d(r\omega )\text{d}{\mathbf{\Lambda }}_d(\omega ),r\ge 0,$$

$${\Omega }_d(z)=\Gamma (d/2)(z/2{)}^{-(d-2)/2}{J}_{(d-2)/2}(z),z\ge 0,$$

Full text

Keywords: Geodesic; Gneiting classes; Latent dimensions; Precipitations; Space-time; Sphere; Temperature.

1 Introduction

Monitoring several georeferenced variables is a common practice in a wide range of disciplines such as climatology and oceanography. The phenomena under study are often observed over large portions of the Earth and across several time instants. Since there is only a finite sample from the involved variables, geostatistical models are a useful tool to capture both the spatial and temporal interactions between the observed data, as well as the uncertainty associated to the limited available information (Cressie, , 1993; Wackernagel, , 2003; Gneiting et al., , 2007). The geostatistical approach consists in modeling the observations as a partial realization of a space-time multivariate random field (RF), denoted as $\{\bm{Z}(\bm{x},t)=(Z_{1}(\bm{x},t),\ldots,Z_{m}(\bm{x},t))^{\top}:(\bm{x},t)\in\mathcal{D}\times\mathcal{T}\}$, where $\top$ is the transpose operator and $m$ is a positive integer representing the number of components of the field. If $m=1$, we say that $Z(\bm{x},t)$ is a univariate (or scalar-valued) RF, whereas for $m>1$, $\bm{Z}(\bm{x},t)$ is called an $m$-variate (or vector-valued) RF. Here, $\mathcal{D}$ and $\mathcal{T}$ denote the spatial and temporal domains, respectively. Throughout, we assume that $\bm{Z}(\bm{x},t)$ is a zero mean Gaussian field, so that a suitable specification of its covariance structure is crucial to develop both accurate inferences and predictions over unobserved sites (Cressie, , 1993).

Parametric families of matrix-valued covariance functions are typically given in terms of Euclidean distances. The literature for this case is extensive and we refer the reader to the review by Genton and Kleiber, (2015) with the references therein. The main motivation to consider the Euclidean metric is the existence of several methods for projecting the geographical coordinates, longitude and latitude, onto the plane. However, when a phenomenon is observed over large portions of planet Earth, the approach based on projections generates distortions in the distances associated to distant locations on the globe. The reader is referred to Banerjee, (2005), where the impact of the different types of projections with respect to spatial inference is discussed.

Indeed, the geometry of the Earth must be considered. Thus, it is more realistic to work under the framework of RFs defined spatially on a sphere (see Marinucci and Peccati, , 2011). Let $d$ be a positive integer. The $d$-dimensional unit sphere is denoted as $\mathbb{S}^{d}:=\{\bm{x}\in\mathbb{R}^{d+1}:\|\bm{x}\|=1\}$, where $\|\cdot\|$ represents the Euclidean norm. The most accurate metric in the spherical scenario is the geodesic (or great circle) distance, which roughly speaking corresponds to the arc joining any two points located on the sphere, measured along a path on the spherical surface. Formally, the geodesic distance is defined as the mapping $\theta:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow[0,\pi]$ given by $\theta:=\theta(\bm{x},\bm{y})=\arccos(\bm{x}^{\top}\bm{y})$.

The construction of valid and flexible parametric covariance functions in terms of the geodesic distance is a challenging problem and requires the application of the theory of positive definite functions on spheres (Schoenberg, , 1942; Yaglom, , 1987; Hannan, , 2009; Berg and Porcu, , 2017). In the univariate and merely spatial case, Huang et al., (2011) study the validity of some specific covariance functions. The essay by Gneiting, (2013) contains a wealth of results related to the validity of a wide range of covariance families. Other related works are the study of star-shaped random particles (Hansen et al., , 2011) and convolution roots (Ziegel, , 2014). However, the spatial and spatio-temporal covariances in the multivariate case are still unexplored, with the work of Porcu et al., (2016) being a notable exception.

The Gneiting class (Gneiting, , 2002) is one of the most popular space-time covariance families and some adaptations in terms of geodesic distance have been given by Porcu et al., (2016). In this paper, we extend to the multivariate scenario the modified Gneiting class introduced by Porcu et al., (2016). Furthermore, we adapt the latent dimension approach to the spherical context (Apanasovich and Genton, , 2010) and we then generate additional Gneiting type matrix-valued covariances. The proposed models are non-separable with respect to the components of the field nor with respect to the space-time interactions. To obtain these results, we have demonstrated several technical results that can be useful to develop new research in this area. Our findings are illustrated through a real data application. In particular, we analyze a bivariate space-time data set of surface air temperatures and precipitations. These variables have been obtained from the Community Climate System Model (CCSM) provided by NCAR (National Center for Atmospheric Research).

The remainder of the article is organized as follows. Section 2 introduces some concepts and construction principles of cross-covariances on spheres across time. In Section 3 we propose some multivariate Gneiting type covariance families. Section 4 contains a real data example of surface air temperatures and precipitations. Section 5 concludes the paper with a discussion.

2 Fundaments and Principles of Space-time Matrix-valued Covariances

This section is devoted to introduce the basic material related to multivariate fields on spheres across time. Let $\bm{Z}(\bm{x},t)$, for $\bm{x}\in\mathbb{S}^{d}$ and $t\in\mathbb{R}$, be an $m$-variate space-time field and let $\textbf{\text{C}}:[0,\pi]\times\mathbb{R}\rightarrow\mathbb{R}^{m\times m}$ be a continuous matrix-valued mapping, whose elements are defined as $\text{C}_{ij}(\theta,u)=\text{cov}\{Z_{i}(\bm{x},t+u),Z_{j}(\bm{y},t)\}$, where $\theta=\theta(\bm{x},\bm{y})$. Following Porcu et al., (2016), we say that C is geodesically isotropic in space and stationary in time. Throughout, the diagonal elements $\text{C}_{ii}$ are called marginal covariances, whereas the off-diagonal members $\text{C}_{ij}$ are called cross-covariances. Any parametric representation of C must respect the positive definite condition, i.e. it must satisfy

[TABLE]

for all positive integer $n$, $\{(\bm{x}_{1},t_{1}),\ldots,(\bm{x}_{n},t_{n})\}\subset\mathbb{S}^{d}\times\mathbb{R}$ and $\{\bm{a}_{1},\ldots,\bm{a}_{n}\}\subset\mathbb{R}^{m}$. Appendix A contains a complete background material on matrix-valued positive definite functions on Euclidean and spherical domains.

We call the mapping C space-time $m$-separable if there exists two mappings $\textbf{\text{C}}_{{\tt S}}:[0,\pi]\rightarrow\mathbb{R}^{m\times m}$ and $\textbf{\text{C}}_{{\tt T}}:\mathbb{R}\rightarrow\mathbb{R}^{m\times m}$, being merely spatial and temporal matrix-valued covariances, respectively, such that

[TABLE]

where $\circ$ denotes the Hadamard product. We call the space-time $m$-separability property complete if there exists a symmetric, positive definite matrix $\textbf{\text{A}}\in\mathbb{R}^{m\times m}$, a univariate spatial covariance $\text{C}_{{\tt S}}:[0,\pi]\rightarrow\mathbb{R}$ , and a univariate temporal covariance $\text{C}_{{\tt T}}:\mathbb{R}\rightarrow\mathbb{R}$, such that

[TABLE]

We finally call $m$-separable the mapping C for which there exists a univariate space-time covariance $\text{C}:[0,\pi]\times\mathbb{R}\rightarrow\mathbb{R}$, and a matrix A, as previously defined, such that

[TABLE]

and of course the special case $\text{C}(\theta,u)=\text{C}_{{\tt S}}(\theta)\text{C}_{{\tt T}}(u)$ offers complete space-time $m$-separability as previously discussed.

Separability is a very useful property for both modeling and estimation purposes, because the related covariance matrices admit nice factorizations, with consequent alleviation of the computational burdens. At the same time, it is a very unrealistic assumption and the literature has been focussed on how to develop non-separable models. How to escape from separability is a major deal, and we list some strategies that can be adapted from others proposed in Euclidean spaces. To the knowledge of the authors, none of these strategies have ever been implemented on spheres or spheres across time.

Linear Models of Coregionalization

Let $q$ be a positive integer. Given a collection of matrices $\textbf{\text{A}}_{k}$, $k=1,\ldots,q$, and a collection of univariate space-time covariances $\text{C}_{k}:[0,\pi]\times\mathbb{R}\rightarrow\mathbb{R}$, the linear model of coregionalization (LMC) has expression

[TABLE]

where a simplification of the type $\text{C}_{k}(\theta,u)=\text{C}_{k,{\tt S}}(\theta)\text{C}_{{\tt T}}(u)$ might be imposed. This model has several drawbacks that have been discussed in Gneiting et al., (2010) as well as in Daley et al., (2015).

Lagrangian Frameworks

Let $\bm{Z}$ be an $m$-variate Gaussian field on the sphere with covariance $\textbf{\text{C}}_{{\tt S}}:[0,\pi]\rightarrow\mathbb{R}^{m\times m}$. Let ${\cal R}$ be a random orthogonal matrix with a given probability law. Let

[TABLE]

Then, $\bm{Y}$ is a RF with transport effect over the sphere. Clearly, $\bm{Y}$ is not Gaussian and evaluation of the corresponding covariance might be cumbersome, as shown in Alegría and Porcu, (2017).

Multivariate Parametric Adaptation

Let $p$ be a positive integer. Let $\text{C}(\cdot,\cdot;\bm{\lambda})$, for $\bm{\lambda}\in\mathbb{R}^{p}$, be a univariate space-time covariance. Let $\bm{\lambda}_{ij}\in\mathbb{R}^{p}$, for $i,j=1,\ldots,m$. For $|\rho_{ij}|\leq 1$ and $\rho_{ii}=1$, and $\sigma_{ii}>0$, find the parametric restrictions such that $\textbf{\text{C}}:[0,\pi]\times\mathbb{R}\rightarrow\mathbb{R}^{m\times m}$ defined through

[TABLE]

is a valid covariance. Sometimes the restriction on the parameters can be very severe, in particular when $m$ is bigger than 2. In Euclidean spaces this strategy has been adopted by Gneiting et al., (2010) and Apanasovich et al., (2012) for the Matérn model, and by Daley et al., (2015) for models with compact support.

Scale Mixtures

Let $(\Omega,{\cal A},\mu)$ be a measure space. Let $\textbf{\text{C}}:[0,\pi]\times\mathbb{R}\times\Omega\rightarrow\mathbb{R}^{m\times m}$ such that

1. 1.

$\textbf{\text{C}}(\cdot,\cdot,\xi)$ is a valid covariance for all $\xi$ in $\Omega$; 2. 2.

$\textbf{\text{C}}(\theta,u,\cdot)\in L_{1}\left(\Omega,{\cal A},\mu\right)$ for all $(\theta,u)\in[0,\pi]\times\mathbb{R}$.

Then,

[TABLE]

is still a valid covariance (Porcu and Zastavnyi, , 2011). Of course, simple strategies can be very effective. For instance, one might assume that $\textbf{\text{C}}(\cdot,\cdot,\xi)=\text{C}(\cdot,\cdot)\textbf{\text{A}}(\xi)$, with C a univariate covariance and $\textbf{\text{A}}(\xi)\in\mathbb{R}^{m\times m}$ being a positive definite matrix for any $\xi\in\Omega$ and such that the hypothesis of integrability above is satisfied.

3 Matrix-valued Covariance Functions of Gneiting Type

This section provides some general results for the construction of multivariate covariances for fields defined on $\mathbb{S}^{d}\times\mathbb{R}$. The most important feature of the proposed models is that they are non-separable mappings, allowing more flexibility in the study of space-time phenomena. In particular, we focus on covariances with a Gneiting’s structure (Gneiting, , 2002). See also Bourotte et al., (2016) for multivariate Gneiting classes on Euclidean spaces.

Recently, Porcu et al., (2016) have proposed some versions of the Gneiting model for RFs with spatial locations on the unit sphere. Let us provide a brief review. The mapping $\text{C}:[0,\pi]\times\mathbb{R}\rightarrow\mathbb{R}$ defined as

[TABLE]

is called Gneiting model. Arguments in Porcu et al., (2016) show that C is a valid univariate covariance, for all $d\in\mathbb{N}$, for $g:[0,\infty)\rightarrow[0,\infty)$ being a completely monotone function, i.e. $g$ is infinitely differentiable on $(0,\infty)$ and $(-1)^{n}g^{(n)}(t)\geq 0$, for all $n\in\mathbb{N}$ and $t>0$, and $f:[0,\infty)\rightarrow(0,\infty)$ is a positive function with completely monotone derivative. Such functions $f$ are called Bernstein functions (Porcu and Schilling, , 2011). Here, $f(\theta)|_{[0,\pi]}$ denotes the restriction of the mapping $f$ to the interval $[0,\pi]$. Tables 3.1 and 3.2 contain some examples of completely monotone and Bernstein functions, respectively. Additional properties about these functions are studied in Porcu and Schilling, (2011).

We also pay attention to the modified Gneiting class (Porcu et al., , 2016) defined through

[TABLE]

where $n\leq 3$ is a positive integer, $g:[0,\infty)\rightarrow[0,\infty)$ is a completely monotone function and $f:[0,\infty)\rightarrow(0,\infty)$ is a strictly increasing and concave function on the positive real line. The mapping (3.2) is a valid covariance for any $d\leq 2n+1$.

3.1 Multivariate Modified Gneiting Class on Spheres Acros Time

We are now able to illustrate the main result within this section. Specifically, we propose a multivariate space-time class generating Gneiting functions with different scale parameters.

Theorem 3.1**.**

Let $m\geq 2$ and $n\leq 3$ be positive integers. Let $g:[0,\infty)\to[0,\infty)$ be a completely monotone function. Consider $f:[0,\infty)\to(0,\infty)$ being strictly increasing and concave. Let $\sigma_{ii}>0$, $|\rho_{ij}|\leq 1$ and $c_{ij}>0$, for $i,j=1,\ldots,m$, be constants yielding the additional condition

[TABLE]

for each $i=1,\ldots,m$. Then, the mapping C, whose members $\text{C}_{ij}:[0,\pi]\times\mathbb{R}\to\mathbb{R}$ are defined through

[TABLE]

is a matrix-valued covariance for any $d\leq 2n+1$.

The condition (3.1) comes from the arguments in Daley et al., (2015). The proof of Theorem 3.1 is deferred to Appendix B, coupled with the preliminary notation and results introduced in Appendix A.

For example, taking $n=1$, $g(t)=\exp(-3t)$ and $f(t)=1+1.7t/c_{T}$, for $c_{T}>0$, we can use Equation (3.2), and the restrictions of Theorem 3.1, to generate a model of the form

[TABLE]

The special parameterization used in the covariance (3.3) ensures that $\text{C}_{ij}(\theta,0)/(\sigma_{ii}\sigma_{jj})<0.05$ and $\text{C}_{ij}(0,|u|)/(\sigma_{ii}\sigma_{jj})<0.05$, for $\theta>c_{ij}$ and $|u|>c_{T}$, respectively.

3.2 A Multivariate Gneiting Family Based on Latent Dimension Approaches

We now propose families of matrix-valued covariances obtained from univariate models. The method used is known as latent dimension approach and has been studied in the Euclidean case by Porcu et al., (2006), Apanasovich and Genton, (2010) and Porcu and Zastavnyi, (2011). Next, we illustrate the spirit of this approach.

Consider a univariate Gaussian RF defined on the product space $\mathbb{S}^{d}\times\mathbb{R}\times\mathbb{R}^{k}$, for some positive integers $d$ and $k$, namely $\{Y(\bm{x},t,\bm{\xi}):(\bm{x},t,\bm{\xi})\in\mathbb{S}^{d}\times\mathbb{R}\times\mathbb{R}^{k}\}$. Suppose that there exists a mapping $\text{K}:[0,\pi]\times\mathbb{R}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ such that $\text{cov}\{Y(\bm{x},t,\bm{\xi}_{1}),Y(\bm{y},s,\bm{\xi}_{2})\}=\text{K}(\theta(\bm{x},\bm{y}),t-s,\bm{\xi}_{1}-\bm{\xi}_{2})$, for all $\bm{x},\bm{y}\in\mathbb{S}^{d}$, $t,s\in\mathbb{R}$ and $\bm{\xi}_{1},\bm{\xi}_{2}\in\mathbb{R}^{k}$. The idea is to define the components of an $m$-variate RF on $\mathbb{S}^{d}\times\mathbb{R}$ as

[TABLE]

for $i=1,\ldots,m$. Thus, the resulting covariance, $\textbf{\text{C}}(\cdot)=[\text{C}_{ij}(\cdot)]_{i,j=1}^{m}$, associated to $\bm{Z}(\bm{x},t)$ is given by

[TABLE]

Here, the vectors $\bm{\xi}_{i}$ are handled as additional parameters of the model. The following theorem allows to construct different versions of the Gneiting model, using the latent dimension technique.

Theorem 3.2**.**

Consider the positive integers $d$, $k$ and $l$. Let $g$ be a completely monotone function and $f_{i}$, $i=1,2$, Bernstein functions. Then,

[TABLE]

is a univariate covariance for fields defined on $\mathbb{S}^{d}\times\mathbb{R}^{l}\times\mathbb{R}^{k}$.

The proof of Theorem 3.2 requires technical lemmas and is deferred to Appendix C, coupled with the preliminary results introduced in Appendix A.

In order to avoid an excessive number of parameters, we follow the parsimonious strategy proposed by Apanasovich and Genton, (2010). Consider $k=1$ and the scalars $\{\xi_{1},\ldots,\xi_{m}\}$. We can consider the parameterization $\zeta_{ij}=|\xi_{i}-\xi_{j}|^{2}$, with $\zeta_{ij}=\zeta_{ji}>0$ and $\zeta_{ii}=0$, for all $i$ and $j$.

Following Apanasovich and Genton, (2010), a LMC based on the latent dimension approach can be used to achieve different marginal structures. Indeed, suppose that the components of a bivariate field are given by $Z_{1}(\bm{x},t)=a_{11}Y(\bm{x},t,\xi_{1})$ and $Z_{2}(\bm{x},t)=a_{21}Y(\bm{x},t,\xi_{2})+a_{22}W(\bm{x},t)$, where $Y(\bm{x},t,\xi)$ is a RF with covariance

[TABLE]

generated from Equation (3.1) and the first entries in Tables 3.1 and 3.2. The field $W(\bm{x},t)$ is independent of $Y(\bm{x},t,\xi)$, with covariance function of Gneiting type as in Equation (3.1),

[TABLE]

Thus, the covariance function of $\bm{Z}(\bm{x},t)$, with ${\rm K}$ and ${\rm R}$ defined by (3.2) and (3.3), respectively, is given by

[TABLE]

Note that $\text{R}(\theta,0)<0.05$ and $\text{R}(0,|u|)<0.05$, for $\theta>c_{S,2}$ and $|u|>c_{T,2}$, respectively. Moreover, since $\text{C}_{12}(0,0)=a_{11}a_{21}/(\zeta_{12}+1)^{1/2}$, the parameter $\zeta_{12}$ is related to the cross-scale as well as the collocated correlation coefficient between the fields.

4 Data Example: Temperatures and Precipitations

We illustrate the use of the proposed models on a climatological data set. We consider a bivariate space-time data set of surface air temperatures (Variable 1) and precipitations (Variable 2) obtained from a climate model provided by NCAR (National Center for Atmospheric Research), Boulder, CO, USA. Specifically, the data come from the Community Climate System Model (CCSM4.0) (see Gent et al., , 2011). Temperatures and precipitations are physically related variables and their joint modeling is crucial in climatological disciplines.

In order to attenuate the skewness of precipitations, we work with its cubic root. The units are degrees Kelvin for temperatures and centimeters for the transformed precipitations. The spatial grid is formed by longitudes and latitudes with $2.5\times 2.5$ degrees of spatial resolution (10368 grid points). We model planet Earth as a sphere of radius 6378 kilometers. We focus on December, between the years 2006 and 2015 (10 time instants), and the region with longitudes between 50 and 150 degrees, and latitudes between $-50$ and 0 degrees. The maximum distance between the points inside this region is approximately half of the maximum distance between any pair of points on planet Earth. Indeed, we are considering a large portion of the globe. For each variable and time instant, we use splines to remove the cyclic patterns along longitudes and latitudes. The residuals are approximately Gaussian distributed with zero mean (see Figure 4.1). We attenuate the computational burden by taking a sub grid with 10 degrees of resolution in the longitud direction. We thus work with only 200 spatial sites and the resulting data set consists of 4000 observations. The empirical variances are $2.033$ and $4.809\times 10^{-5}$, for Variables 1 and 2, respectively. In addition, these variables are moderately positively correlated. The empirical correlation between the components is $0.279$.

We are interested in showing that in real applications a non-separable model can produce significant improvements with respect to a separable one. For our purposes, we consider the following models:

• (A)

An $m$-separable covariance based on the univariate Gneiting class (3.1):

[TABLE]

where the vector of parameters is given by $(\sigma_{11}^{2},\sigma_{22}^{2},\rho_{12},c_{S},c_{T})^{\top}$.

• (B)

An $m$-separable version of the modified Gneiting class (3.3) with the parameterization $c_{S}:=c_{11}=c_{22}=c_{12}$. The vector of parameters is the same as in Model A.

• (C)

A non-separable version of the modified Gneiting class (3.3), with the parsimonious parameterization $c_{12}=\max\{c_{11},c_{22}\}$. The vector of parameters is given by $(\sigma_{11}^{2},\sigma_{22}^{2},\rho_{12},c_{11},c_{22},c_{T})^{\top}$.

• (D)

The non-separable LMC (3.4) based on the latent dimension approach. Here, the vector of parameters is given by $(a_{11},a_{21},a_{22},c_{S,1},c_{S,2},c_{T,1},c_{T,2},\zeta_{12})^{\top}$.

We estimate the models parameters using the pairwise composite likelihood (CL) method developed for multivariate RFs by Bevilacqua et al., (2016). This method is a likelihood approximation and offers a trade-off between statistical and computational efficiency. We only consider pairs of observations with spatial distances less than $1275.6$ kilometers (equivalent to 0.2 radians on the unit sphere) and temporal distance less than $4$ years.

The CL estimates for Models A-C are reported in Table 4.1, whereas Table 4.2 reports the results for Model D. In addition, Table 4.3 contains the Log-CL values attained at the optimum for each model. As expected, the non-separable covariances C and D exhibit the highest values of the Log-CL. In Figure 4.2, we compare the empirical spatial semi-variograms, for the temporal lags $|u|=0$ and $|u|=1$, versus the theoretical models C and D. Both models seem to capture the behavior of the empirical variograms at the limit values (sills). However, as noted by Gneiting et al., (2010), disagreements between empirical and theoretical fits are typically observed in practice, and it can be associated to biases in the empirical estimators.

Finally, we compare the predictive performance of the covariances through a cross-validation study based on the cokriging predictor. We use a drop-one prediction strategy and quantify the discrepancy between the real and predicted values for each variable through the Mean Squared Error (MSE) and the Continuous Ranked Probability Score (CRPS) (Zhang and Wang, , 2010). Smaller values of these indicators imply better predictions. Table 4.3 displays the results, where $\text{MSE}_{i}$ denotes the MSE associated to the variable $i$, for $i=1,2$. The interpretation of $\text{CRPS}_{i}$ is similar. Apparently, the performance of Model A is quite poor in comparison to the other models. Moreover, the non-separable models, C and D, have better results than the $m$-separable model B. Note that although Model B has a smaller $\text{MSE}_{2}$ value than Model D, D is globally superior. If we choose C instead of B, we have a significant improvement in the prediction of Variable 1, since the corresponding ${\rm MSE}_{1}$ ratio is approximately 0.84. Similarly, the ${\rm CRPS}_{1}$ ratio is 0.93. Also, direct inspection of ${\rm MSE}_{2}$ and ${\rm CRPS}_{2}$ show that Model C provides improvements in the prediction of Variable 2 in comparison to Model B. Indeed, for this specific data set, Model C shows the best results in terms of MSE, outperforming even Model D. A nice property of the multivariate modified Gneiting class C is that it has physically interpretable parameters. Naturally, the proposed models can be combined to provide more flexible models.

5 Discussion

In the paper, we have discussed several construction principles for multivariate space-time covariances, which allow to escape from separability. In particular, we have proposed Gneiting type families of cross-covariances and their properties have been illustrated through a real data example of surface air temperatures and precipitations. The proposed models have shown a good performance using separable models as a benchmark. We believe that the methodology used to prove our theoretical results can be adapted to find additional flexible classes of matrix-valued covariances and to develop new applications on univariate or multivariate global phenomena evolving temporally. Specifically, our proposals can be used as a building block for the construction of more sophisticated models, such as non-stationary RFs.

Acknowledgments

Alfredo Alegría’s work was initiated when he was a PhD student at Universidad Técnica Federico Santa María. Alfredo Alegría was supported by Beca CONICYT-PCHA/Doctorado Nacional/2016-21160371. Emilio Porcu is supported by Proyecto Fondecyt Regular number 1170290. Reinhard Furrer acknowledges support of the Swiss National Science Foundation SNSF-144973. Jorge Mateu is partially supported by grant MTM2016-78917-R. We additionally acknowledge the World Climate Research Programme’s Working Group on Coupled modeling, which is responsible for Coupled Model Intercomparison Project (CMIP).

Appendices

A Background on Positive Definite Functions

We start by defining positive definite functions, that arise in statistics as the covariances of Gaussian RFs as well as the characteristic functions of probability distributions.

Definition A.1**.**

Let $\mathcal{E}$ be a non-empty set and $m\in\mathbb{N}$. We say that the matrix-valued function $\textbf{\text{F}}:\mathcal{E}\times\mathcal{E}\rightarrow\mathbb{R}^{m\times m}$ is positive definite if for all integer $n\geq 1$, $\{e_{1},\ldots,e_{n}\}\subset\mathcal{E}$ and $\{\bm{a}_{1},\ldots,\bm{a}_{n}\}\subset\mathbb{R}^{m}$, the following inequality holds:

[TABLE]

We denote as $\mathcal{P}^{m}(\mathcal{E})$ the class of such mappings F satisfying Equation (A.1).

Next, we focus on the cases where $\mathcal{E}$ is either $\mathbb{R}^{d}$, $\mathbb{S}^{d}$ or $\mathbb{S}^{d}\times\mathbb{R}^{k}$, for $d,k\in\mathbb{N}$. For a clear presentation of the results, Table A.1 summarizes the notation introduced along this Appendix.

A.1 Matrix-valued Positive Definite Functions on Euclidean Spaces: The Classes $\Phi^{m}_{d,\mathcal{S}}$ and $\Phi^{m}_{d,\mathcal{I}}$

This section provides a brief review about matrix-valued positive definite functions on the Euclidean space $\mathcal{E}=\mathbb{R}^{d}$. Specifically, we expose some characterizations for the stationary and Euclidean isotropic members of the class $\mathcal{P}^{m}(\mathbb{R}^{d})$.

We say that $\textbf{\text{F}}\in\mathcal{P}^{m}(\mathbb{R}^{d})$ is stationary if there exists a mapping $\tilde{\bm{\varphi}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m\times m}$ such that

[TABLE]

We call $\Phi_{d,\mathcal{S}}^{m}$ the class of continuous mappings $\tilde{\bm{\varphi}}$ such that F in (A.2) is positive definite. Cramér’s Theorem (Cramer, , 1940) establishes that $\tilde{\bm{\varphi}}\in\Phi_{d,\mathcal{S}}^{m}$ if and only if it can be represented through

[TABLE]

where $\imath=\sqrt{-1}\in\mathbb{C}$ and $\tilde{\bm{\Lambda}}_{d}:\mathbb{R}^{d}\rightarrow\mathbb{C}^{m\times m}$ is a matrix-valued mapping, with increments being Hermitian and positive definite matrices, and whose elements, $\tilde{\Lambda}_{d,ij}(\cdot)$, for $i,j=1,\ldots,m$, are functions of bounded variation (see Wackernagel, , 2003). In particular, the diagonal terms, $\tilde{\Lambda}_{d,ii}(\bm{\omega})$, are real, non-decreasing and bounded, whereas the off-diagonal elements are generally complex-valued. Cramer’s Theorem is the multivariate version of the celebrated Bochner’s Theorem (Bochner, , 1955). If the elements of $\tilde{\bm{\Lambda}}_{d}(\cdot)$ are absolutely continuous, then Equation (A.3) simplifies to

[TABLE]

with $\tilde{\bm{\lambda}}_{d}(\bm{\omega})=[\tilde{\lambda}_{d,ij}(\bm{\omega})]_{i,j=1}^{m}$ being Hermitian and positive definite, for any $\bm{\omega}\in\mathbb{R}^{d}$. The mapping $\tilde{\bm{\lambda}}_{d}(\bm{\omega})$ is known as the matrix-valued spectral density and classical Fourier inversion yields

[TABLE]

Finally, the following inequality between the elements of $\tilde{\bm{\varphi}}$ is true

[TABLE]

However, the maximum value of the mapping $\tilde{\varphi}_{ij}(\bm{h})$, with $i\neq j$, is not necessarily reached at $\bm{h}=\bm{0}$. In general, $\tilde{\varphi}_{ij}$ is not itself a scalar-valued positive definite function when $i\neq j$.

Consider an element F in $\mathcal{P}^{m}(\mathbb{R}^{d})$ and suppose that there exists a continuous and bounded mapping $\bm{\varphi}:\mathbb{R}_{+}\rightarrow\mathbb{R}^{m\times m}$ such that

[TABLE]

Then, F is called stationary and Euclidean isotropic (or radial). We denote as $\Phi^{m}_{d,\mathcal{I}}$ the class of bounded, continuous, stationary and Euclidean isotropic mappings $\bm{\varphi}(\cdot)=[\varphi_{ij}(\cdot)]_{i,j=1}^{m}$.

When $m=1$, characterization of the class $\Phi_{d,\mathcal{I}}$ was provided through the celebrated paper by Schoenberg, (1938). Alonso-Malaver et al., (2015) characterize the class $\Phi^{m}_{d,\mathcal{I}}$ through the continuous members $\bm{\varphi}$ having representation

[TABLE]

where $\bm{\Lambda}_{d}:[0,\infty)\rightarrow\mathbb{R}^{m\times m}$ is a matrix-valued mapping, with increments being positive definite matrices, and elements $\Lambda_{d,ij}(\cdot)$ of bounded variation, for each $i,j=1,\ldots,m$. Here, the function $\Omega_{d}(\cdot)$ is defined as

[TABLE]

with $\Gamma$ being the Gamma function and $J_{\nu}$ the Bessel function of the first kind of degree $\nu$ (see Abramowitz and Stegun, , 1970). If the elements of $\bm{\Lambda}_{d}(\cdot)$ are absolutely continuous, then we have an associated spectral density $\bm{\lambda}_{d}:[0,\infty)\rightarrow\mathbb{R}^{m\times m}$ as in the stationary case, which is called, following Daley and Porcu, (2014), a $d$-Schoenberg matrix.

The classes $\Phi^{m}_{d,\mathcal{I}}$ are non-increasing in $d$, and the following inclusion relations are strict

[TABLE]

The elements $\bm{\varphi}$  in the class $\Phi^{m}_{\infty,\mathcal{I}}$ can be represented as

[TABLE]

where $\bm{\Lambda}$ is a matrix-valued mapping with similar properties as $\bm{\Lambda}_{d}$.

A.2 Matrix-valued Positive Definite Functions on $\mathbb{S}^{d}$: The Class $\Psi_{d,\mathcal{I}}^{m}$

In this section, we pay attention to matrix-valued positive definite functions on the unit sphere. Consider $\textbf{\text{F}}=[\text{F}_{ij}]_{i,j=1}^{m}\in\mathcal{P}^{m}(\mathbb{S}^{d})$. We say that F is geodesically isotropic if there exists a bounded and continuous mapping $\bm{\psi}:[0,\pi]\rightarrow\mathbb{R}^{m\times m}$ such that

[TABLE]

The continuous mappings $\bm{\psi}$ are the elements of the class $\Psi_{d,\mathcal{I}}^{m}$ and the following inclusion relations are true:

[TABLE]

where $\Psi_{\infty,\mathcal{I}}^{m}$ is the class of geodesically isotropic positive definite functions being valid on the Hilbert sphere $\mathbb{S}^{\infty}=\{(x_{n})_{n\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}:\sum_{n\in\mathbb{N}}x_{n}^{2}=1\}$.

The elements of the class $\Psi_{d,\mathcal{I}}^{m}$ have an explicit connection with Gegenbauer (or ultraspherical) polynomials (Abramowitz and Stegun, , 1970). Here, $\mathcal{G}_{n}^{\lambda}$ denotes the $\lambda$-Gegenbauer polynomial of degree $n$, which is defined implicitly through the expression

[TABLE]

In particular, $\mathcal{T}_{n}:=\mathcal{G}_{n}^{0}$ and $\mathcal{P}_{n}:=\mathcal{G}_{n}^{1/2}$ are respectively the Chebyshev and Legendre polynomials of degree $n$.

The following result (Hannan, , 2009; Yaglom, , 1987) offers a complete characterization of the classes $\Psi_{d,\mathcal{I}}^{m}$ and $\Psi_{\infty,\mathcal{I}}^{m}$, and corresponds to the multivariate version of Schoenberg’s Theorem (Schoenberg, , 1942). Equalities and summability conditions for matrices must be understood in a componentwise sense.

Theorem A.1**.**

Let $d$ and $m$ be positive integers.

• (1)

The mapping $\bm{\psi}$ is a member of the class $\Psi_{d,\mathcal{I}}^{m}$ if and only if it admits the representation

[TABLE]

where $\{\textbf{\text{B}}_{n,d}\}_{n=0}^{\infty}$ is a sequence of symmetric, positive definite and summable matrices.

• (2)

The mapping $\bm{\psi}$ is a member of the class $\Psi_{\infty,\mathcal{I}}^{m}$ if and only if it can be represented as

[TABLE]

where $\{\textbf{\text{B}}_{n}\}_{n=0}^{\infty}$ is a sequence of symmetric, positive definite and summable matrices.

Using orthogonality properties of Gegenbauer polynomials (Abramowitz and Stegun, , 1970) and through classical Fourier inversion we can prove that

[TABLE]

whereas for $d\geq 2$, we have

[TABLE]

where integration is taken componentwise. The matrices $\{\textbf{\text{B}}_{n,d}\}_{n=0}^{\infty}$ are called Schoenberg’s matrices. For the case $m=1$, such result is reported by Gneiting, (2013).

A.3 Matrix-valued Positive Definite Functions on $\mathbb{S}^{d}\times\mathbb{R}^{k}$: The Class $\Upsilon_{d,k}^{m}$

Let $d$, $k$ and $m$ be positive integers. We now focus on the class of matrix-valued positive definite functions on $\mathbb{S}^{d}\times\mathbb{R}^{k}$, being bounded, continuous, geodesically isotropic in the spherical component and stationary in the Euclidean one. The case $k=1$ is particularly important, since $\mathcal{P}^{m}(\mathbb{S}^{d}\times\mathbb{R})$ can be interpreted as the class of admissible space-time covariances for multivariate Gaussian RFs, with spatial locations on the unit sphere.

Consider $\textbf{\text{F}}\in\mathcal{P}^{m}(\mathbb{S}^{d}\times\mathbb{R}^{k})$ and suppose that there exists a bounded and continuous mapping $\textbf{\text{C}}:[0,\pi]\times\mathbb{R}^{k}\rightarrow\mathbb{R}^{m\times m}$ such that

[TABLE]

Such mappings C are the elements of the class $\Upsilon_{d,k}^{m}$. These classes are non-increasing in $d$ and we have the inclusions

[TABLE]

Ma, (2016) proposes the generalization of Theorem A.1 to the space-time case. Theorem A.2 below offers a complete characterization of the class $\Upsilon_{d,k}^{m}$ and $\Upsilon_{\infty,k}^{m}$, for any $m\geq 1$. Again, equalities and summability conditions must be understood in a componentwise sense.

Theorem A.2**.**

Let $d$, $k$ and $m$ be positive integers and $\textbf{\text{C}}:[0,\pi]\times\mathbb{R}^{k}\to\mathbb{R}^{m\times m}$ a continuous matrix-valued mapping, with $\text{C}_{ii}(0,\bm{0})<\infty$, for all $i=1,\ldots,m$.

• (1)

The mapping C belongs to the class $\Upsilon_{d,k}^{m}$ if and only if

[TABLE]

with $\left\{\tilde{\bm{\varphi}}_{n,d}(\cdot)\right\}_{n=0}^{\infty}$ being a sequence of members of the class $\Phi_{k,\mathcal{S}}^{m}$, with the additional requirement that the sequence of matrices $\{\tilde{\bm{\varphi}}_{n,d}(\bm{0})\}_{n=0}^{\infty}$ is summable.

• (2)

The mapping C belongs to the class $\Upsilon_{\infty,k}^{m}$ if and only if

[TABLE]

with $\left\{\tilde{\bm{\varphi}}_{n}(\cdot)\right\}_{n=0}^{\infty}$ being a sequence of members of the class $\Phi_{k,\mathcal{S}}^{m}$, with the additional requirement that the sequence of matrices $\{\tilde{\bm{\varphi}}_{n}(\bm{0})\}_{n=0}^{\infty}$ is summable.

Again, using orthogonality arguments, we have

[TABLE]

whereas for $d\geq 2$,

[TABLE]

B Proof of Theorem 3.1

In order to illustrate the results following subsequently, a technical Lemma will be useful. We do not provide a proof because it is obtained following the same arguments as in Porcu and Zastavnyi, (2011).

Lemma .1**.**

Let $m$, $d$ and $k$ be strictly positive integers. Let $(X,{\cal B},\mu)$ be a measure space, for $X\subset\mathbb{R}$ and $\cal B$ being the Borel sigma algebra. Let $\bm{\psi}:[0,\pi]\times X\to\mathbb{R}$ and $\bm{\varphi}:[0,\infty)\times X\to\mathbb{R}$ be continuous mappings satisfying

$\bm{\psi}(\cdot,\xi)\in\Psi_{d,{\cal I}}^{m}$ a.e. $\xi\in X$; 2. 2.

$\bm{\psi}(\theta,\cdot)\in L_{1}(X,{\cal B},\mu)$ for any $\theta\in[0,\pi]$; 3. 3.

$\bm{\varphi}(\cdot,\xi)\in\Phi_{k,{\cal I}}^{m}$ a.e. $\xi\in X$; 4. 4.

$\bm{\varphi}(u,\cdot)\in L_{1}(X,{\cal B},\mu)$ for any $u\in[0,\infty)$.

Let $\textbf{\text{C}}:[0,\pi]\times[0,\infty)\rightarrow\mathbb{R}^{m\times m}$ be the mapping defined through

[TABLE]

Then, C is continuous and bounded. Further, C belongs to the class $\Upsilon_{d,k}^{m}$.

Of course, Lemma .1 is a particular case of the scale mixtures introduced in Section 2.

Proof of Theorem 3.1 Let $(X,{\cal B},\mu)$ as in Lemma .1 and consider $X=\mathbb{R}_{+}$ with $\mu$ the Lebesgue measure. We offer a proof of the constructive type. Let us define the function $\bm{\psi}(\theta,\xi)$ with members $\psi_{ij}(\cdot,\cdot)$ defined through

[TABLE]

where, as asserted, the constants $\sigma_{ii}$, $\rho_{ij}$ and $c_{ij}$ are determined according to condition (3.1). Let us now define the mapping $(u,\xi)\mapsto\varphi(u,\xi):=\xi^{n+1}(1-\xi f(u))_{+}^{\ell}$, with $(u,\xi)\in[0,\infty)\times X$. It can be verified that both $\bm{\psi}$ and $\varphi$ satisfy requirements 1–4 in Lemma .1. In particular, Condition 1 yields thanks to Lemmas 3 and 4 in Gneiting, (2013), as well as Theorem 1 in Daley et al., (2015). Also, arguments in Porcu et al., (2016) show that Condition 3 holds for any $\ell\geq 1$. We can now apply Lemma .1, so that we have that

[TABLE]

is a member of the class $\Upsilon_{2n+1,1}^{m}$ for any $\ell\geq 1$. Pointwise application of an elegant scale mixture argument as in Proposition 1 of Porcu et al., (2016) shows that

[TABLE]

where $\mathcal{B}$ denotes the Beta function (Abramowitz and Stegun, , 1970). We now omit the factor $\mathcal{B}(n+2,\ell+1)$ since it does not affect positive definiteness. Now, standard convergence arguments show that

[TABLE]

with the convergence being uniform in any compact set. The proof is then completed in view of Bernstein’s theorem (Feller, , 1966).

C Proof of Theorem 3.2

Before we state the proof of Theorem 3.2, we need to introduce three auxiliary lemmas.

Lemma .1**.**

Let $\text{C}:[0,\pi]\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ be a continuous, bounded and integrable function, for some positive integer $k$. Then $\text{C}\in\Upsilon_{d,k}$, for $d\geq 1$, if and only if the mapping $\psi_{\bm{\omega}}:[0,\pi]\rightarrow\mathbb{R}$ defined as

[TABLE]

belongs to the class $\Psi_{d,\mathcal{I}}$, for all $\bm{\omega}\in\mathbb{R}^{k}$.

We do not report the proof of Lemma .1 since the arguments are the same as in the proof of Lemma .2 below. Note that this lemma is a spherical version of the result given by Cressie and Huang, (1999).

Lemma .2**.**

Let $\text{C}:[0,\pi]\times\mathbb{R}^{l}\times\mathbb{R}^{k}\rightarrow\mathbb{R}$ be a continuous, bounded and integrable function, for some positive integers $l$ and $k$. Then, $\text{C}\in\Upsilon_{d,k+l}$, with $d\geq 1$, if and only if the mapping $\text{C}_{\bm{\omega}}:[0,\pi]\times\mathbb{R}^{l}\rightarrow\mathbb{R}$ defined as

[TABLE]

belongs to the class $\Upsilon_{d,l}$, for all $\bm{\omega}\in\mathbb{R}^{k}$.

Proof of Lemma .2 Suppose that $\text{C}\in\Upsilon_{d,k+l}$, then the characterization of Berg and Porcu, (2017) implies that

[TABLE]

where $\{\tilde{\varphi}_{n,d}(\cdot,\cdot)\}_{n=0}^{\infty}$ is a sequence of functions in $\Phi_{k+l,\mathcal{S}}$, with $\sum_{n=0}^{\infty}\tilde{\varphi}_{n,d}(\bm{0},\bm{0})<\infty$. Therefore,

[TABLE]

where the last step is justified by dominated convergence. We need to prove that for each fixed $\bm{\omega}\in\mathbb{R}^{k}$, the sequence of functions

[TABLE]

belongs to the class $\Phi_{l,\mathcal{S}}$, a.e. $\bm{\omega}\in\mathbb{R}^{k}$. In fact, we have that

[TABLE]

Since $\tilde{\varphi}_{n,d}(\cdot,\cdot)$ belongs to $\Phi_{k+l,\mathcal{S}}$, Bochner’s Theorem implies that the right side in Equation (C.4) is non-negative everywhere. This implies that $\tilde{\lambda}_{n,d}(\cdot;\bm{\omega})$ belongs to the class $\Phi_{l,\mathcal{S}}$. Also, direct inspection shows that $\sum_{n=0}^{\infty}\tilde{\lambda}_{n,d}(\bm{0};\bm{\omega})<\infty$, for all $\bm{\omega}\in\mathbb{R}^{k}$. The necessary part is completed.

On the other hand, suppose that for each $\bm{\omega}\in\mathbb{R}^{k}$ the function $\text{C}_{\bm{\omega}}(\theta,\bm{u})$ belongs to the class $\Upsilon_{d,l}$, then there exists a sequence of mappings $\{\tilde{\lambda}_{n,d}(\cdot;\bm{\omega})\}_{n=0}^{\infty}$ in $\Phi_{l,\mathcal{S}}$ for each $\bm{\omega}\in\mathbb{R}^{k}$, such that

[TABLE]

Thus,

[TABLE]

We conclude the proof by invoking again Bochner’s Theorem and the result of Berg and Porcu, (2017).

Lemma .3**.**

Let $d$ and $k$ be two positive integers. Consider $g$ and $f$ be completely monotone and Bernstein functions, respectively. Then,

[TABLE]

belongs to the class $\Upsilon_{d,k}$, for any positive integer $d$.

Proof of Lemma .3 By Lemma .1, we must show that $\psi_{\bm{\omega}}$, defined through Equation (C.1), belongs to the class $\Psi_{d,\mathcal{I}}$, for all $\bm{\omega}\in\mathbb{R}^{k}$. In fact, we can assume that C is integrable, since the general case is obtained with the same arguments given by Gneiting, (2002). Bernstein’s Theorem establishes that $g$ can be represented as the Laplace transform of a bounded measure $G$, then

[TABLE]

where the last equality follows from Fubini’s Theorem and $\text{d}G(r)=r^{k/2}\text{d}\tilde{G}(r)$. In addition, the composition between a negative exponential and a Bernstein function is completely monotone on the real line (Feller, , 1966). Then, for any $\bm{\omega}$ and $r$, the mapping $\theta\mapsto\exp\{-\|\bm{\omega}\|^{2}f(\theta)|_{[0,\pi]}/(4r)\}$ is the restriction of a completely monotone function to the interval $[0,\pi]$. Theorem 7 in Gneiting, (2013) implies that such mapping, and thus $\psi_{\bm{\omega}}$, belongs to the class $\Psi_{d,\mathcal{I}}$, for any $d\in\mathbb{N}$ and $\bm{\omega}\in\mathbb{R}^{k}$.

Proof of Theorem 3.2 By Lemma .2, we must show that (C.2) belongs to the class $\Upsilon_{d,l}$, for all $d\in\mathbb{N}$. In fact, assuming again that C is integrable and invoking Bernstein’s Theorem we have

[TABLE]

where the last equality follows from Fubini’s Theorem and $\text{d}G(r)=r^{k/2}\text{d}\tilde{G}(r)$. In addition, for any $\bm{\omega}$, the mapping

[TABLE]

is completely monotone (Feller, , 1966). Therefore,

[TABLE]

and by Lemma .3, we have that $\text{C}_{\bm{\omega}}\in\Upsilon_{d,l}$, for all $d\in\mathbb{N}$.