Asymmetric Matrix-Valued Covariances for Multivariate Random Fields on Spheres

TL;DR

This paper introduces a novel approach for constructing asymmetric matrix-valued covariance functions on spheres, improving the modeling of multivariate spatial data over large geographic areas.

Contribution

It proposes a rotation-based method to generate asymmetric covariances on spheres, addressing the limitations of symmetric models in geostatistics.

Findings

Enhanced predictive performance over symmetric models

Effective in real data analysis and simulations

Addresses the gap in asymmetric covariance modeling on spheres

Abstract

Matrix-valued covariance functions are crucial to geostatistical modeling of multivariate spatial data. The classical assumption of symmetry of a multivariate covariance function is overlay restrictive and has been considered as unrealistic for most of real data applications. Despite of that, the literature on asymmetric covariance functions has been very sparse. In particular, there is some work related to asymmetric covariances on Euclidean spaces, depending on the Euclidean distance. However, for data collected over large portions of planet Earth, the most natural spatial domain is a sphere, with the corresponding geodesic distance being the natural metric. In this work, we propose a strategy based on spatial rotations to generate asymmetric covariances for multivariate random fields on the $d$-dimensional unit sphere. We illustrate through simulations as well as real data analysis that our proposal allows to achieve improvements in the predictive performance in comparison to the symmetric counterpart.

Figures

Tables

Table 1. Table 1: Cross-validation scores for models (M1)–(M3), separable and non-separable, using symmetric (S) and asymmetric (A) versions, with ρ 12 = 0.25 subscript 𝜌 12 0.25 \rho_{12}=0.25 .

			Separable					Non-separable
			$\eta =0.1$		$\eta =0.6$			$\eta =0.1$		$\eta =0.6$
			S	A	S	A		S	A	S	A
(M1)		MSPE	0.983	0.968	0.964	0.935		0.964	0.953	0.954	0.924
		LSCORE	1.410	1.402	1.398	1.383		1.396	1.386	1.387	1.373
(M2)		MSPE	0.969	0.956	0.961	0.930		0.924	0.910	0.929	0.893
		LSCORE	1.398	1.393	1.395	1.381		1.357	1.348	1.352	1.327
(M3)		MSPE	0.981	0.966	0.980	0.959		0.993	0.988	0.972	0.945
		LSCORE	1.409	1.399	1.407	1.396		1.414	1.409	1.405	1.388

Table 2. Table 2: Cross-validation scores for models (M1)–(M3), separable and non-separable, using symmetric (S) and asymmetric (A) versions, with ρ 12 = 0.5 subscript 𝜌 12 0.5 \rho_{12}=0.5 .

			Separable					Non-separable
			$\eta =0.1$		$\eta =0.6$			$\eta =0.1$		$\eta =0.6$
			S	A	S	A		S	A	S	A
(M1)		MSPE	0.969	0.945	0.979	0.919		0.940	0.919	0.950	0.882
		LSCORE	1.401	1.389	1.405	1.369		1.383	1.372	1.387	1.345
(M2)		MSPE	0.971	0.950	0.975	0.908		0.938	0.904	0.924	0.856
		LSCORE	1.400	1.390	1.406	1.366		1.356	1.339	1.358	1.326
(M3)		MSPE	0.992	0.969	0.993	0.911		0.981	0.951	0.970	0.897
		LSCORE	1.415	1.405	1.414	1.364		1.408	1.393	1.403	1.366

Table 3. Table 3: Composite Likelihood estimates for Models 1-4.

	${\hat{\sigma }}_{TT}^2$	${\hat{\sigma }}_{PP}^2$	${\hat{\rho }}_{TP}$	${\hat{c}}_{TT}$	${\hat{c}}_{PP}$	${\hat{c}}_{TP}$	$\hat{\eta }$	${\hat{\alpha }}_1$	${\hat{\alpha }}_2$
Model 1	34.643	$9.262\times {10}^{-5}$	0.342	5102	5102	5102	-	-	-
Model 2	34.275	$9.363\times {10}^{-5}$	0.342	4789	5414	5414	-	-	-
Model 3	34.524	$9.235\times {10}^{-5}$	0.342	5089	5089	5089	0.043	1.512	2.678
Model 4	34.098	$9.439\times {10}^{-5}$	0.342	4815	5427	5427	0.047	1.504	2.236

Table 4. Table 4: Comparison of Log-CL values and cross-validation scores for Models 1-4.

	$\#$ parameters	Log-CL	MSPE	LSCORE
Model 1	4	20460	1.309	$-1.235$
Model 2	5	20485	1.306	$-1.236$
Model 3	7	20500	1.157	$-1.271$
Model 4	8	20518	1.139	$-1.277$

Equations

$$F_{ij}(\mathbf{x},\mathbf{y}):=\text{cov}\{Z_i(\mathbf{x}),Z_j(\mathbf{y})\},\mathbf{x},\mathbf{y}\in {\mathbb{S}}^d.$$

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

$$\frac{1}{(1+r^2-2r\mu {)}^{\lambda }}=\sum _{k=0}^{\infty }{r}^k{\mathcal{P}}_k^{\lambda }(\mu ),\mu \in [-1,1],r\in (-1,1).$$

$$\mathbf{C}(\theta )=\sum _{k=0}^{\infty }{\mathbf{B}}_k{\mathcal{P}}_k^{(d-1)/2}(\cos \theta ),\theta \in [0,\pi ],$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}\frac{2^{1-\nu }}{\Gamma (\nu )}{\left(\frac{\theta }{c_{ij}}\right)}^{\nu }{K}_{\nu }\left(\frac{\theta }{c_{ij}}\right),\theta \in [0,\pi ],$$

$$|{\rho }_{12}|\le {\left(\frac{c_{11}{c}_{22}}{c_{12}^2}\right)}^{\nu }.$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}{\left(1+\frac{{\theta }^{\gamma }}{c_{ij}}\right)}^{-\nu },\theta \in [0,\pi ],$$

$${\rho }_{12}^2\le {\left(\frac{c_{11}{c}_{22}}{c_{12}^2}\right)}^{\nu }.$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}{\left(1-\frac{\theta }{c_{ij}}\right)}_{+}^{\nu }\left(1+\nu \frac{\theta }{c_{ij}}\right),\theta \in [0,\pi ],$$

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

$$Z_i^a(\mathbf{x})=Z_i({\mathbf{R}}_i\mathbf{x}),\mathbf{x}\in {\mathbb{S}}^d,i=1,\dots ,p.$$

$$F_{ij}^a(\mathbf{x},\mathbf{y})$$

$$\mathbf{R}(\delta ):=\left[\begin{array}{cc}\cos \delta & \sin \delta \\ -\sin \delta & \cos \delta \end{array}\right],-\pi <\delta \le \pi ,$$

$${\delta }_1+\cdots +{\delta }_p=0.$$

$${\mathbf{F}}^a(\mathbf{x},\mathbf{y})=\sum _{k=0}^{\infty }{\mathbf{B}}_k\circ {\mathbf{P}}_k(\mathbf{x},\mathbf{y}),$$

$${\mathbf{P}}_k(\mathbf{x},\mathbf{y})=\left(\begin{array}{cc}{\mathcal{P}}_k^0(\cos \theta )& {\mathcal{P}}_k^0(\cos \theta \cos \eta +(x_1{y}_2-x_2{y}_1)\sin \eta )\\ {\mathcal{P}}_k^0(\cos \theta \cos \eta -(x_1{y}_2-x_2{y}_1)\sin \eta )& {\mathcal{P}}_k^0(\cos \theta )\end{array}\right),$$

$${\mathbf{R}}_{\mathbf{\omega }}(\delta )={\mathbf{I}}_3+(\sin \delta )\mathbf{\Omega }+(1-\cos \delta ){\mathbf{\Omega }}^2,$$

$$\mathbf{\Omega }=\left[\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right].$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}\exp \left(-\frac{3\theta }{c_{ij}}\right),\theta \in [0,\pi ],i,j=1,2,$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}{\left(1+\frac{19\theta }{c_{ij}}\right)}^{-1},\theta \in [0,\pi ],i,j=1,2,$$

$$C_{ij}(\theta )={\sigma }_i{\sigma }_j{\rho }_{ij}{\left(1-\frac{\theta }{c_{ij}}\right)}_{+}^4\left(1+4\frac{\theta }{c_{ij}}\right),\theta \in [0,\pi ],i,j=1,2,$$

Full text

1 Introduction

In recent years, the use of multivariate geostatistical models has become increasingly popular in disciplines such as climatology, environmental sciences and oceanography. Specifically, the georeferenced observations are assumed as a partial realization of a multivariate Gaussian random field (MGRF) and the spatial dependencies between the observed data are modeled through the matrix-valued covariance function of the field. Covariance functions are the main tool for spatial interpolation methods and the development of accurate predictions. The best linear unbiased predictor, under a mean squared error criterion, is known in multivariate geostatistics as the co-kriging predictor. We refer the reader to [18] for a detailed study on MGRFs and [8] for a review on parametric families of matrix-valued covariance functions.

A typical assumption in geostatistics is that of spatial isotropy, which means that the multivariate covariance depends solely on spatial distance. The practical and theoretical contributions in spatial statistics have been mainly developed on Euclidean spaces with the corresponding Euclidean metric. This choice is supported by the existence of several methods to project the geographical locations onto the plane. However, for data recorded over large portions of planet Earth, such projection generate significant distortions in the distances. The problem of an appropriate metric has been discussed by [4], [9] and [17]. Indeed, if we agree that a sphere is a good approximation of planet Earth, it is more realistic to work under the framework of random fields indexed on a sphere [16]. Moreover, the most natural metric to describe the spatial interactions is the geodesic distance, which corresponds to the shortest path joining two points on the spherical surface.

The construction of flexible families of matrix-valued covariance functions on spheres is a challenging problem, since the corresponding covariance matrices must be positive definite. The literature in this direction has been very sparse with the work of [17], who derive conditions for the validity of the multivariate Matérn and Cauchy models, being a notable exception. On the other hand, [15] shows that compactly supported models, introduced on Euclidean spaces, can be adopted on specific low dimensional spheres. Additional approaches based on latent dimensions are explored by [2].

The parametric covariance families mentioned above depend solely on the geodesic distance (geodesic isotropy). In particular, spatial isotropy implies symmetry of the matrix-valued covariance structure. However, in real data applications, asymmetrical behaviors are usually observed. Indeed, we need to go beyond the class of geodesically isotropic covariances to generate asymmetric models. In this paper, we propose a strategy based on spatial rotations to generate asymmetric covariances from symmetric models. We show that the inclusion of asymmetry produce more flexibility in the statistical modeling, and particularly, that the improvements in terms of predictive performance are significant. Our proposal is the spherical counterpart of the work developed by [14] on Euclidean spaces.

The remainder of the article is as follows. Section 2 contains the background material on MGRFs on spheres as well as a brief summary of some symmetric covariance families, which are the building block to construct asymmetric models. Section 3 introduces a general approach to produce asymmetric covariances. In Section 4, we show through simulation experiments that asymmetric models can provide better predictive results in comparison to the symmetric counterparts. A bivariate data set is analyzed in Section 5. Finally, Section 6 contains some conclusions.

2 Background

We briefly review the literature on MGRFs on spheres and the parametric families of matrix-valued covariance functions introduced in the literature. Although the models considered in this section are symmetric, they are the building block to construct explicit asymmetric models in the subsequent sections.

2.1 MGRFs on Spheres

Throughout the article, we consider MGRFs, $\{\bm{Z}(\bm{x})=(Z_{1}(\bm{x}),\ldots,Z_{p}(\bm{x}))^{\top}:\bm{x}\in\mathbb{S}^{d}\}$, defined on the $d$-dimensional unit sphere $\mathbb{S}^{d}=\{\bm{x}\in\mathbb{R}^{d+1}:\|\bm{x}\|=1\}$, $d\in\mathbb{N}$, where $p\in\mathbb{N}$ denotes the number of components of the field, ⊤ the transpose operator and $\|\cdot\|$ the Euclidean norm in $\mathbb{R}^{d+1}$. We call $\bm{Z}(\bm{x})$ a $p$-variate random field. The natural metric on the spherical surface is the geodesic distance defined through 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})$.

We assume, without loss of generality, that $\bm{Z}(\bm{x})$ has zero mean, since we are essentially interested on its covariance structure. Let $i,j\in\{1,\ldots,p\}$ and consider the matrix-valued mapping ${\bm{F}}:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow\mathbb{R}^{p\times p}$, with $(i,j)$-th element defined through

[TABLE]

The diagonal elements, ${F}_{ii}:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow\mathbb{R}$, are called marginal covariances, whereas the off-diagonal elements, ${F}_{ij}:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow\mathbb{R}$, for $i\neq j$, are called cross-covariances. The mapping ${\bm{F}}$ is (semi) positive definite, i.e., for all positive integer $n$, $\{\bm{x}_{1},\ldots,\bm{x}_{n}\}\subset\mathbb{S}^{d}$ and $\{\bm{a}_{1},\ldots,\bm{a}_{n}\}\subset\mathbb{R}^{p}$, the following inequality holds

[TABLE]

Note that, in general, the cross-covariance function is not symmetric in the sense that $F_{ij}(\bm{x},\bm{y})\neq F_{ji}(\bm{x},\bm{y})$ and $F_{ij}(\bm{x},\bm{y})\neq F_{ij}(\bm{y},\bm{x})$. However, when this covariance depends on $\bm{x}$ and $\bm{y}$ only through its geodesic distance (so-called geodesic isotropy), such a covariance function is clearly symmetric. The rest of this section focuses on geodesically isotropic models.

Following [17], we call $\Psi_{d}^{p}$ the class of continuous mappings $\bm{C}:[0,\pi]\rightarrow\mathbb{R}^{p\times p}$, being continuous and such that, for a positive definite function ${\bm{F}}:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow\mathbb{R}^{p\times p}$, we have ${\bm{F}}(\bm{x},\bm{y})=\bm{C}(\theta(\bm{x},\bm{y}))$. The mapping $\bm{C}$ is called radial part of a geodesically isotropic covariance ${\bm{F}}$. Moreover, the inclusion $\Psi_{d+1}^{p}\subset\Psi_{d}^{p}$, for all positive integer $d$, is strict.

Throughout, $\mathcal{P}_{k}^{\lambda}(\cdot)$ denotes the $\lambda$-ultraspherical polynomial of degree $k$ [1], which is defined implicitly as

[TABLE]

The following result [11, 19] characterizes completely the class $\Psi_{d}^{p}$ in terms of ultraspherical polynomials. The equalities and summability conditions must be understood in a componentwise sense.

Theorem 2.1**.**

Let $d$ and $p$ be positive integers. Let $\bm{C}:[0,\pi]\rightarrow\mathbb{R}^{p\times p}$ be a continuous mapping with elements $C_{ij}(\cdot)$, such that $C_{ii}(0)<+\infty$, for all $i=1,\ldots,p$. Then, $\bm{C}$ is a member of the class $\Psi_{d}^{p}$ if and only if it admits the representation

[TABLE]

where $\{\bm{B}_{k}\}_{k=0}^{\infty}$ is a sequence of symmetric, positive definite and summable matrices of order $p\times p$.

2.2 Parametric Families in $\Psi_{d}^{p}$

We now provide a review on parametric families of geodesically isotropic matrix-valued covariance functions. We focus on the Matérn and Cauchy families, as well as on compactly supported models. In particular, the Matérn and Cauchy families are valid on spheres of any dimension, whereas the compactly supported models require additional restrictions on the dimension.

Matérn model

Let $K_{\nu}$ be the modified Bessel function of second kind and $\Gamma$ be the Gamma function [1]. The multivariate Matérn model on the sphere [17] is given by

[TABLE]

provided that $\sigma_{i}$ and $c_{ij}$ are positive, $0<\nu\leq 1/2$, $|\rho_{ij}|\leq 1$, with $\rho_{ii}=1$, for all $i,j=1,\ldots,p$. Here, $\sigma_{i}$ is the standard deviation of the $i$-th component of the field, $c_{ij}$ are scale parameters, $\rho_{ij}$ is the collocated correlation coefficient between the components $i$ and $j$, and $\nu$ is a smoothness parameter. Some additional conditions for the parameters are required (see [17]). In the following sections, we focus on parsimonious bivariate ($p=2$) models and such conditions are relaxed. For instance, the parameterization $c_{12}=\max\{c_{11},c_{22}\}$ has associated condition

[TABLE]

For a more detailed study of the parametric conditions, we refer the reader to [10], [3] and [17].

Cauchy model

The conditions for matrix-valued covariance functions of generalized Cauchy type on spheres have been developed by [17]. The generalized Cauchy model is given by

[TABLE]

with $\sigma_{i}$, $c_{ij}$ and $\nu$ being positive, $0<\gamma\leq 1$, $|\rho_{ij}|\leq 1$ and $\rho_{ii}=1$, for all $i,j=1,\ldots,p$. Here, the interpretation of the parameters is similar to the Matérn model. Additional conditions must be imposed in order to obtain a positive definite mapping (see [17]). For instance, in a bivariate case with $c_{12}=(c_{11}+c_{22})/2$, we have the restriction

[TABLE]

Moreover, [17] provide explicit conditions on the parameters for the $p$-variate case, with $p>2$.

Compactly supported models

[15] shows that a radial matrix-valued positive definite function on $\mathbb{R}^{3}$, with elements having compact support less than $\pi$, can be coupled with the geodesic distance to obtain a geodesically isotropic positive definite function on $\mathbb{S}^{3}$.

Formally, consider the mapping $\bm{\varphi}:[0,+\infty)\rightarrow\mathbb{R}^{p\times p}$, such that $(\bm{x},\bm{y})\mapsto\bm{\varphi}(\|\bm{x}-\bm{y}\|)$, with elements $\varphi_{ij}(\cdot)$, is positive definite on $\mathbb{R}^{3}$, with $\varphi_{ij}(r)=0$, for all $r>\pi$ and $i,j=1,\ldots,p$. Therefore, the mapping $\bm{C}:[0,\pi]\rightarrow\mathbb{R}^{p\times p}$ defined by $\bm{C}(\theta)=\bm{\varphi}(\theta)|_{[0,\pi]}$ belongs to $\Psi_{3}^{p}$, where $f|_{A}$ denotes the restriction of $f$ to a Borel set $A\subset\mathbb{R}$.

In particular, the inclusion between the classes $\Psi_{d}^{p}$ imply that, under this construction principle, we can obtain models on spheres of dimension 1, 2 and 3. Examples can be generated from Table 1 in [6]. For instance, the Wendland model is defined as

[TABLE]

where $\nu\geq 2$, $0<c_{ij}\leq\pi$, $\sigma_{i}>0$, $|\rho_{ij}|\leq 1$ and $\rho_{ii}=1$, for all $i,j=1,\ldots,p$. Here $(a)_{+}=\max\{a,0\}$ and the interpretation of the parameters is similar to the previous examples. In addition, we have the condition (see [6])

[TABLE]

In the subsequent sections, we consider a bivariate Wendland model with the parsimonious choice $c_{12}=\max\{c_{11},c_{22}\}$.

The characterization of the cross scale parameter $c_{ij}$ as a function of the marginal ones, $c_{ii}$ and $c_{jj}$, is a very useful strategy to avoid an excessive number of parameters in the estimation procedure. In some cases, this parsimonious strategy can provide better predictive results than the full model (see [10]).

Additional parametric models can be constructed using latent dimensions [2], or by extending the spherical convolution approach [21] as well as covariance approximation methods [12] to the multivariate case.

3 Asymmetric Covariances on $\mathbb{S}^{d}$

A strategy to construct asymmetric models is now provided. We start with a $p$-variate zero mean Gaussian field, $\bm{Z}(\bm{x})$, with geodesically isotropic, and thus symmetric, covariance structure $\bm{F}(\bm{x},\bm{y})=\bm{C}(\theta(\bm{x},\bm{y}))\in\Psi_{d}^{p}$.

Let $\{{\bm{R}}_{1},\ldots,{\bm{R}}_{p}\}\subset\mathbb{R}^{(d+1)\times(d+1)}$ be a collection of rotation matrices, i.e., ${\bm{R}}_{i}$ is an orthogonal matrix with determinant one, for each $i=1,\ldots,p$. Such matrices satisfy the relation ${\bm{R}}_{i}^{\top}={\bm{R}}_{i}^{-1}$. Consider a new $p$-variate random field $\{\bm{Z}^{a}(\bm{x})=(Z_{1}^{a}(\bm{x}),\ldots,Z_{p}^{a}(\bm{x}))^{\top}:\bm{x}\in\mathbb{S}^{d}\}$ defined through

[TABLE]

The field introduced in Equation (4) is Gaussian and it has zero mean. Furthermore, the covariance function associated to $\bm{Z}^{a}(\bm{x})$, denoted as $\bm{F}^{a}:\mathbb{S}^{d}\times\mathbb{S}^{d}\rightarrow\mathbb{R}^{p\times p}$, is given by

[TABLE]

for each $i,j=1,\ldots,p$. Obviously, $\bm{Z}^{a}(\bm{x})$ preserves the marginal structure associated to $\bm{Z}(\bm{x})$, since ${F}^{a}_{ii}(\bm{x},\bm{y})=C_{ii}(\theta(\bm{x},\bm{y}))$, for each $i=1,\ldots,p$. However, the cross-covariances ${F}_{ij}^{a}$, for $i\neq j$, are not necessarily geodesically isotropic. Thus, $\bm{F}^{a}$ does not belong to the class $\Psi_{d}^{p}$. Moreover, we have that ${F}_{ij}^{a}(\bm{x},\bm{y})\neq{F}^{a}_{ji}(\bm{x},\bm{y})$, for $\bm{x},\bm{y}\in\mathbb{S}^{d}$.

Here, the rotation matrices are handled as additional parameters of the model. In order to avoid identifiability problems in the estimation of such parameters, we need to add some constraints in a similar fashion to the work of [14]. Since only the product ${\bm{R}}_{i}^{-1}{\bm{R}}_{j}$ is relevant, we propose the condition ${\bm{R}}_{1}\cdots{\bm{R}}_{p}={\bm{I}}_{d+1}$, where ${\bm{I}}_{d+1}$ denotes the identity matrix of order $(d+1)\times(d+1)$. We illustrate these considerations for the cases $d=1$ and $d=2$.

Example 3.1**.**

In the unit circle $\mathbb{S}^{1}$, we have the general rotation matrix

[TABLE]

where $\delta$ denotes the rotation angle. For each $i=1,\ldots,p$, we define ${\bm{R}}_{i}:={\bm{R}}(\delta_{i})$. Thus, the constraint ${\bm{R}}_{1}\cdots{\bm{R}}_{p}={\bm{R}}(\delta_{1}+\cdots+\delta_{p})={\bm{I}}_{2}$ is equivalent to the condition $(\delta_{1}+\cdots+\delta_{p})\text{ mod }2\pi=0$. In particular, we use the following identifiability condition

[TABLE]

Direct inspection shows that ${\bm{R}}_{i}^{-1}{\bm{R}}_{j}={\bm{R}}(\delta_{j}-\delta_{i})$ and $\bm{x}^{\top}{\bm{R}}_{i}^{-1}{\bm{R}}_{j}\bm{y}=\cos(\theta)\cos(\delta_{j}-\delta_{i})+(x_{1}y_{2}-x_{2}y_{1})\sin(\delta_{j}-\delta_{i})$, where $\theta:=\theta(\bm{x},\bm{y})$ and $\bm{x}=(x_{1},x_{2})^{\top},\bm{y}=(y_{1},y_{2})^{\top}\in\mathbb{S}^{1}$.

For instance, consider a bivariate scenario with $\delta_{1}=\eta/2$ and $\delta_{2}=-\eta/2$, for $-\pi<\eta<\pi$. In this case, we have the additional asymmetry parameter $\eta$. The general expression for the bivariate asymmetric covariance functions on $\mathbb{S}^{1}$, generated from our approach, is given by

[TABLE]

where

[TABLE]

with $\circ$ denoting the Hadamard product. Recall that $\{\bm{B}_{k}\}_{k=0}^{\infty}$ is a summable sequence of symmetric positive definite matrices. In this case, the ultraspherical polynomials reduces to the Tchebyshev polynomials $\mathcal{P}_{k}^{0}(\cos\theta)=\cos(k\theta)$.

Example 3.2**.**

Consider now the unit sphere $\mathbb{S}^{2}$. The Rodrigues’s rotation formula provides the following representation for a rotation matrix, with respect to an axis $\bm{\omega}=(\omega_{1},\omega_{2},\omega_{3})^{\top}\in\mathbb{S}^{2}$ by an angle $\delta\in(-\pi,\pi]$,

[TABLE]

where

[TABLE]

Here, we take the collection of matrices ${\bm{R}}_{i}={\bm{R}}_{\bm{\omega}}(\delta_{i})$. Note that we have the relation ${\bm{R}}_{i}^{-1}{\bm{R}}_{j}={\bm{R}}_{\bm{\omega}}(\delta_{j}-\delta_{i})$. For the bivariate case ($p=2$), considering the identifiability condition (5), $\delta_{1}=\eta/2$ and $\delta_{2}=-\eta/2$, we have the additional parameters $(\eta,\bm{\omega}^{\top})^{\top}\in(-\pi,\pi)\times\mathbb{S}^{2}$. Since the pairs $(\eta,\bm{\omega}^{\top})$ and $(-\eta,-\bm{\omega}^{\top})$ define the same rotation, we must to restrict the range of $\eta$ to the interval $[0,\pi)$. The axis of rotation can be parameterized in terms of spherical coordinates: $\bm{\omega}=(\cos\alpha_{1}\sin\alpha_{2},\sin\alpha_{1}\sin\alpha_{2},\cos\alpha_{2})^{\top}$, with $0\leq\alpha_{1}<2\pi$ and $0\leq\alpha_{2}\leq\pi$. In conclusion, we have the vector of additional parameters $(\eta,\alpha_{1},\alpha_{2})^{\top}$.

4 Simulation Study

We illustrate simulation experiments related to bivariate Gaussian fields on $\mathbb{S}^{2}$ with asymmetric covariances. We use initial covariance functions either of Matérn, Cauchy and Wendland type, as introduced in the previous sections. In particular, we consider the mentioned models with the following parameterizations:

• (M1)

Matérn model with $\nu=1/2$, which reduces to the exponential model

[TABLE]

where $c_{12}=\max\{c_{11},c_{22}\}$.

• (M2)

Cauchy model with $\nu=1$,

[TABLE]

where $c_{12}=(c_{11}+c_{22})/2$.

• (M3)

Wendland model with $\nu=4$,

[TABLE]

where $c_{12}=\max\{c_{11},c_{22}\}$.

The parameterization of (M1) and (M2) ensures that $(C_{ij}(\theta)/C_{ij}(0))<0.05$ for $\theta>c_{ij}$. For model (M3), $(C_{ij}(\theta)/C_{ij}(0))$ is exactly equal to 0 beyond the cut-off distance $c_{ij}$. The large number of parameters considered in our simulation experiments supports the parsimonious parameterization used for $c_{12}$. In particular, for the generalized Cauchy model (M2) we consider $c_{12}$ as the average of the marginal scale parameters, since under this assumption we have the condition (2), which can be checked easily. The same happens for the other models (see conditions (1) and (3)). Figure 1 shows these covariances under a specific parameter setting: $\sigma_{1}^{2}=\sigma_{2}^{2}=1$, $\rho_{12}=0.5$, $c_{11}=0.1$, $c_{22}=0.2$.

In the special case where $c_{ij}=c>0$, for all $i,j=1,2$, we say that the covariances are separable, otherwise they are non-separable. We make simulation experiments for both cases. For each non-separable model, after adding asymmetry, the vector of parameters is given by $(\sigma_{1}^{2},\sigma_{2}^{2},\rho_{12},c_{11},c_{22},\eta,\alpha_{1},\alpha_{2})^{\top}$. In the separable case, on the other hand, the vector of parameters reduces to $(\sigma_{1}^{2},\sigma_{2}^{2},\rho_{12},c,\eta,\alpha_{1},\alpha_{2})^{\top}$. In our studies, we consider 225 spatial sites in a grid generated combining 15 equispaced azimuthal and polar angles. Specifically, the polar angles are taken as $2\pi(k-1)/15$, whereas the azimuthal angles as $\pi(k-1)/15$, for $k=1,\ldots,15$.

We assess the typical variance of the likelihood-based estimates when the asymmetric effect is added. In particular, we consider the pairwise composite likelihood (CL) approach developed by [5], which is a fast likelihood approximation method. We set $\sigma_{1}^{2}=\sigma_{2}^{2}=1$, $\rho_{12}=0.5$, $c_{11}=0.1$, $c_{22}=0.2$, $\alpha_{1}=\alpha_{2}=\pi/2$ and $\eta=0.1$. Under the separable scenario, we set $c=0.1$. Figure 2 reports the bias of the CL estimates for each model on the basis of 500 independent simulated data sets. We have used a special case of CL method, by considering only pairs of data whose spatial distance is less than 1 radian, although in our experience, smaller cut-off distances can also ensure stability of the method. Even though we are dealing with a large number of parameters, the boxplots show that the direct maximization of the CL objective function provides quite reasonable results. The empirical $95\%$ confidence interval for $\eta$ is $[0.055,0.190]$ in the non-separable model (M1), and similar intervals are obtained for the other models. Another appealing estimation procedure is to consider a two steps maximization routine (see [14]), which is useful in terms of computational complexity.

Now, we study the predictive performance of asymmetric models through a cross-validation analysis. We simulate bivariate fields from the asymmetric versions of models (M1)–(M3), with both separable and non-separable parameterizations. We set $\sigma_{1}^{2}=\sigma_{2}^{2}=1$, $c_{11}=0.1$, $c_{22}=0.2$, $\alpha_{1}=\alpha_{2}=\pi/2$, as in the previous example, and consider two scenarios for the collocated correlation coefficient: $\rho_{12}=0.25,0.5$, and two scenarios for the asymmetry parameter: $\eta=0.1,0.6$. Again, in the separable case, we consider a similar experiment with $c=0.1$. We quantify the accuracy of the co-kriging predictor, with both the symmetric and asymmetric covariances, in terms of the mean squared prediction error (MSPE) and the log-score (LSCORE) (see [20]). These indicators are evaluated using a drop-one prediction strategy. Tables 1 and 2 report the results for $\rho_{12}=0.25$ and $\rho_{12}=0.5$, respectively. It is clearly visible that the asymmetric models produce better results. Note that if the asymmetry and the correlation between the fields are strong, then the improvements are more significant.

5 A Real Data Example

This section assesses the statistical performance of asymmetric models on a real data example. We analyze a bivariate data set of Temperatures (T) and Precipitations (P). These data outputs have been obtained from the Community Climate System Model (CCSM4.0, [7]) provided by National Center for Atmospheric Research (NCAR), CO, USA. The joint modeling of these variables is of major interest in climatological disciplines.

The spatial resolution is 2.5 degrees in latitude and longitude. We assume planet Earth as a sphere of radius 6378 kilometers. We pay attention to the geographical region delimited by the longitudes $[0,100]$ and latitudes $[-40,60]$ degrees (see Figure 3) and consider December 2015. The maximum geodesic distance between the points in this region is 70% of the maximum distance on planet Earth. We work with the cubic root of precipitations in order to attenuate the skewness of this variable. The units are Kelvin degrees for temperatures and centimeters for transformed precipitations. In order to achieve geodesic isotropy, we remove through splines the trend and cyclic patterns of the variables, using longitude and latitude as covariates. The resulting residuals are approximately Gaussian with zero mean.

We consider the bivariate exponential model (M1) with the following features:

• •

Model 1. Separable exponential model without asymmetry.

• •

Model 2. Non-separable exponential model without asymmetry.

• •

Model 3. Separable exponential model after adding asymmetry.

• •

Model 4. Non-separable exponential model after adding asymmetry.

Table 3 illustrates the pairwise CL estimates for each model. We have considered CL method with pairs of data whose geodesic distance is less than 1000 kilometers. The estimates of the scale parameters are given in kilometers. We compare the models in terms of the Log-Composite Likelihood (Log-CL) value at the optimum, as well as through their predictive performance. The accuracy of the prediction, for each model, is evaluated in terms of MSPE and LSCORE, using a drop-one prediction strategy (see Table 4). It is clear that the asymmetric versions provide better results with respect to the symmetric ones. Note that for this specific data set, Model 3 outperforms Model 2, which means that the inclusion of asymmetry can produce more flexibility in comparison to the inclusion of non-separability.

6 Conclusions

We have provided a strategy to construct asymmetric matrix-valued covariance functions on spheres. We show that our proposal can produce significant improvements in the predictive performance of the model. Our developments have been illustrated through simulated data, under several models and parametric settings. We have also analyzed a real data set of temperatures and precipitations on a large portion of planet Earth, where the asymmetric models produce better results.

An interesting research direction is the construction of parametric tests, based on the asymmetry parameter $\eta$, to contrast the hypothesis of symmetry in a given data set. In particular, a challenging problem is the adaptation of the methodology used by [13], where the authors evaluate several types of common assumptions on multivariate covariance functions, to the spherical context.

Acknowledgments

Alfredo Alegría is supported by Beca CONICYT-PCHA/Doctorado Nacional/2016-21160371. Emilio Porcu is supported by Proyecto Fondecyt Regular number 1130647. Reinhard Furrer acknowledges support of the Swiss National Science Foundation SNSF-144973. We additionally acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for Coupled Model Intercomparison Project (CMIP).