Covariance Functions for Multivariate Gaussian Fields evolving temporally over Planet Earth

5.2 MB·

by  Alfredo Alegría

Covariance Functions for Multivariate Gaussian Fields evolving
temporally over Planet Earth
Alfredo Alegr´ıa*, Emilio Porcu*, Reinhard Furrer**, and Jorge Mateu***
7
1
0
* Department of Mathematics, University Federico Santa Maria, Valparaiso, Chile.
2

n **Department of Mathematics and Department of Computer Science, University of Zurich,
a
J Zurich, Switzerland.

1
2 ***Department of Mathematics, University Jaume I, Castello´n, Spain.

]
T
S January 24, 2017
.
h
t
a
Abstract
m
[
Theconstructionofvalidandflexiblecross-covariancefunctionsisafundamentaltaskformodelling

1
multivariate space-time data arising from climatological and oceanographical phenomena. Indeed,
v
0
asuitablespecificationofthecovariancestructureallowstocaptureboththespace-timedependen-
1
0 cies between the observations and the development of accurate predictions. For data observed over
6
0 largeportionsofplanetEarthitisnecessarytotakeintoaccountthecurvatureoftheplanet. Hence
.
1 the need for random field models defined over spheres across time. In particular, the associated
0
7 covariance function should depend on the geodesic distance, which is the most natural metric over
1
: thesphericalsurface. Inthiswork,weprovideacompletecharacterisationforthespace-timecross-
v
i covariances, that are continuous, geodesically isotropic in the spatial variable and stationary in the
X
r temporal one, and discuss some fundaments and construction principles of matrix-valued covari-
a
ances on spheres. 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 simulation experi-
ments, and we analyse a bivariate space-time data set of surface air temperatures and atmospheric
pressures.
Keywords: Atmospheric pressure; Geodesic; Gneiting classes; Latent dimensions; Space-time;
Sphere; Temperature.
1
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
modelling the observations as a partial realisation of a space-time multivariate random field (RF),
denoted as
{Z(x,t) = (Z (x,t),...,Z (x,t))(cid:62) : (x,t) ∈ D×T},
1 m
where(cid:62)isthetransposeoperatorandmisapositiveintegerrepresentingthenumberofcomponentsof
the field. If m = 1, we say that Z is a univariate (or scalar-valued) random field, whereas for m > 1, Z
is called an m-variate (or vector-valued) random field. Here, D and T denote the spatial and temporal
domains, respectively. Throughout, we assume that Z is Gaussian, so that a suitable specification
of the covariance structure of Z is crucial to develop both accurate inferences and predictions over
unobserved sites (Cressie, 1993).
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 the whole planet Earth, the approach based on projections
generates apparent 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 random fields defined spatially on a sphere (see
Marinucci and Peccati, 2011).
Let d be a positive integer. The d-dimensional unit sphere is denoted as Sd := {x ∈ Rd+1 : (cid:107)x(cid:107) = 1},
where (cid:107)·(cid:107) 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
2
distance is defined as the mapping θ : Sd×Sd → [0,π] given by
θ := θ(x,y) = arccos(x(cid:62)y).
The construction of valid and flexible 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, 2016). 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 works of Berg and Porcu (2016) and Porcu et al.
(2016) being notable exceptions.
In this work, we propose some progress in this area. First, we give a complete characterisation for the
space-time matrix-valued covariances, being continuous, geodesically isotropic in the spatial variable
andstationaryinthetemporalone. Ourcharacterisationgeneralisestothemultivariatecasetheresult
given by Berg and Porcu (2016). Also, we discuss several interesting aspects, such as different notions
of separability and some construction principles, of cross-covariances on Sd×R.
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). We extend to
the multivariate scenario the modified Gneiting class introduced by Porcu et al. (2016). Furthermore,
we adapt to the spherical context the latent dimension approach (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.
Ourfindingsareillustratedthroughsimulationstudiesaswellasarealdataapplication. Inparticular,
we analyse a bivariate space-time data set of surface air temperatures and atmospheric pressures.
These variables have been generated from the Community Climate System Model (CCSM) provided
by NCAR (National Center for Atmospheric Research). We compare the performance of the proposed
models to the traditional linear model of coregionalisation.
3
The remainder of the article is organised as follows. Section 2 introduces the main notation and some
preliminaryresultsonpositivedefinitefunctions. InSection3,wegiveacharacterisationforthespace-
time matrix-valued covariance functions, defined spatially in terms of the geodesic distance. Some
fundaments of space-time cross-covariances on spheres are additionally discussed. Section 4 provides
some multivariate Gneiting type families. In order to illustrate the properties of the proposed models,
Section 5 contains simulation experiments and a real data example of surface air temperatures and
atmospheric pressures.
2 Background
This section introduces the main notation and preliminaries of the paper. We start with positive
definite functions, that arise in statistics as the covariances of Gaussian random fields as well as the
characteristic functions of probability distributions.
Definition 2.1. Let E be a non-empty set and m ∈ N. We say that the matrix-valued function
F : E×E → Rm×m is positive definite if for all integer n ≥ 1, {e ,...,e } ⊂ E and {a ,...,a } ⊂ Rm,
1 n 1 n
the following inequality holds:
n n
(cid:88)(cid:88)
a(cid:62)F(e ,e )a ≥ 0. (2.1)
(cid:96) (cid:96) r r
(cid:96)=1r=1
We denote as Pm(E) the class of such mappings F satisfying (2.1).
Next, we focus on the cases where E is either Rd, Sd or Sd×Rk, for d,k ∈ N. For a clear presentation
of the results, Table 2.1 summarises the notation introduced along the paper.
2.1 Matrix-valued positive definite functions on Euclidean spaces
This section provides a brief review about matrix-valued positive definite functions on the Euclidean
spaceE = Rd. Specifically,weexposesomecharacterisationsforthestationaryandEuclideanisotropic
members of the class Pm(Rd).
4
Table2.1: Summaryofthenotationusedalongthepaper. Throughout, intheunivariatecase(m = 1)
we omit the super index: P(E), Φ , Φ , Ψ and Υ are used instead of P1(E), Φ1 , Φ1 , Ψ1
d,S d,I d,I d,k d,S d,I d,I
and Υ1 , respectively.
d,k
Notation Description
Pm(E) Positive definite matrix-valued (m×m) functions on E ×E.
Φm Continuous, bounded and stationary elements of Pm(Rd).
d,S
Φm Continuous, bounded and Euclidean isotropic elements of Pm(Rd).
d,I
Ψm Continuous, bounded and geodesically isotropic elements of Pm(Sd).
d,I
Υm Elements in Pm(Sd×Rk) being, continuous, bounded, geodesically
d,k
isotropic in the spherical variable and stationary in the Euclidean variable.
Stationarity: the class Φm
d,S
We say that F ∈ Pm(Rd) is stationary if there exists a mapping ϕ˜ : Rd → Rm×m such that
F(x,y) = ϕ˜(x−y) = [ϕ˜ (x−y)]m , x,y ∈ Rd. (2.2)
ij i,j=1
We call Φm the class of continuous mappings ϕ˜ such that F in (2.2) is positive definite. Cram´er’s
d,S
Theorem (Cramer, 1940) establishes that ϕ˜ ∈ Φm if and only if it can be represented through
d,S
(cid:90)
ϕ˜(h) = exp{ıh(cid:62)ω}dΛ˜ (ω), h ∈ Rd, (2.3)
d
Rd
√
where ı = −1 ∈ C and Λ˜ : Rd → Cm×m is a matrix-valued mapping, with increments being
d
Hermitianandpositivedefinitematrices,andwhoseelements,Λ˜ (·),fori,j = 1,...,m,arefunctions
d,ij
ofboundedvariation(seeWackernagel,2003). Inparticular,thediagonalterms,Λ˜ (ω),arereal,non-
d,ii
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 Λ˜ (·) are absolutely continuous, then Equation (2.3) simplifies to
d
(cid:90)
ϕ˜(h) = exp{ıh(cid:62)ω}λ˜ (ω)dω, h ∈ Rd,
d
Rd
with λ˜ (ω) = [λ˜ (ω)]m being Hermitian and positive definite, for any ω ∈ Rd. The mapping
d d,ij i,j=1
λ˜ (ω) is known as the matrix-valued spectral density and classical Fourier inversion yields
d
(cid:90)
1
λ˜ (ω) = exp{−ıh(cid:62)ω}ϕ˜(h)dh, ω ∈ Rd.
d (2π)d
Rd
5
Finally, the following inequality between the elements of ϕ˜ is true
|ϕ˜ (h)|2 ≤ ϕ˜ (0)ϕ˜ (0), for all h ∈ Rd.
ij ii jj
However, the maximum value of the mapping ϕ˜ (h), with i (cid:54)= j, is not necessarily reached at h = 0.
ij
In general, ϕ˜ is not itself a scalar-valued positive definite function when i (cid:54)= j.
ij
Euclidean isotropy: the class Φm
d,I
Consider an element F in Pm(Rd) and suppose that there exists a continuous and bounded mapping
ϕ : R → Rm×m such that
+
F(x,y) = ϕ((cid:107)x−y(cid:107)), x,y ∈ Rd.
Then, F is called stationary and Euclidean isotropic (or radial). We denote as Φm the class of
d,I
bounded, continuous, stationary and Euclidean isotropic mappings ϕ(·) = [ϕ (·)]m .
ij i,j=1
When m = 1, characterisation of the class Φ was provided through the celebrated paper by Schoen-
d,I
berg (1938). Alonso-Malaver et al. (2015) characterise the class Φm through the continuous members
d,I
ϕ having representation
(cid:90)
ϕ(r) = Ω (rω)dΛ (ω), r ≥ 0,
d d
[0,∞)
where Λ : [0,∞) → Rm×m is a matrix-valued mapping, with increments being positive definite
d
matrices, and elements Λ (·) of bounded variation, for each i,j = 1,...,m. Here, the function Ω (·)
d,ij d
is defined as
Ω (z) = Γ(d/2)(z/2)−(d−2)/2J (z), z ≥ 0,
d (d−2)/2
with Γ being the Gamma function and J the Bessel function of the first kind of degree ν (see
ν
Abramowitz and Stegun, 1970). If the elements of Λ (·) are absolutely continuous, then we have an
d
associated spectral density λ : [0,∞) → Rm×m as in the stationary case, which is called, following
d
Daley and Porcu (2014), a d-Schoenberg matrix.
The classes Φm are non-increasing in d, and the following inclusion relations are strict
d,I
∞
(cid:92)
Φm := Φm ⊂ ··· ⊂ Φm ⊂ Φm .
∞,I d,I 2,I 1,I
d=1
6
The elements ϕ in the class Φm can be represented as
∞,I
(cid:90)
ϕ(r) = exp(−r2ω2)dΛ(ω), r ≥ 0,
[0,∞)
where Λ is a matrix-valued mapping with similar properties as Λ .
d
2.2 Matrix-valued positive definite functions on Sd and the class Ψm
d,I
In this section, we pay attention to matrix-valued positive definite functions on the unit sphere.
Consider F = [F ]m ∈ Pm(Sd). We say that F is geodesically isotropic if there exists a bounded
ij i,j=1
and continuous mapping ψ : [0,π] → Rm×m such that
F(x,y) = ψ(θ(x,y)), x,y ∈ Sd.
The continuous mappings ψ are the elements of the class Ψm and the following inclusion relations
d,I
are true:
∞
(cid:92)
Ψm = Ψm ⊂ ··· ⊂ Ψm ⊂ Ψm , (2.4)
∞,I d,I 2,I 1,I
d=1
where Ψm is the class of geodesically isotropic positive definite functions being valid on the Hilbert
∞,I
sphere S∞ = {(xn)n∈N ∈ RN : (cid:80)n∈Nx2n = 1}.
The elements of the class Ψm have an explicit connection with Gegenbauer (or ultraspherical) poly-
d,I
nomials (Abramowitz and Stegun, 1970). Here, Gλ denotes the λ-Gegenbauer polynomial of degree n,
n
which is defined implicitly through the expression
∞
1 (cid:88)
= rnGλ(cosθ), θ ∈ [0,π], r ∈ (−1,1).
(1+r2−2rcosθ)λ n
n=0
In particular, T := G0 and P := G1/2 are respectively the Chebyshev and Legendre polynomials of
n n n n
degree n.
The following result (Hannan, 2009; Yaglom, 1987) offers a complete characterisation of the classes
Ψm and Ψm , and corresponds to the multivariate version of Schoenberg’s Theorem (Schoenberg,
d,I ∞,I
1942). Equalities and summability conditions for matrices must be understood in a componentwise
sense.
Theorem 2.1. Let d and m be positive integers.
7
(1) The mapping ψ is a member of the class Ψm if and only if it admits the representation
d,I
∞ (d−1)/2
(cid:88) Gn (cosθ)
ψ(θ) = B , θ ∈ [0,π], (2.5)
n,d (d−1)/2
G (1)
n=0 n
where {B }∞ is a sequence of symmetric, positive definite and summable matrices.
n,d n=0
(2) The mapping ψ is a member of the class Ψm if and only if it can be represented as
∞,I
∞
(cid:88)
ψ(θ) = B (cosθ)n, θ ∈ [0,π],
n
n=0
where {B }∞ is a sequence of symmetric, positive definite and summable matrices.
n n=0
UsingorthogonalitypropertiesofGegenbauerpolynomials(AbramowitzandStegun,1970)andthrough
classical Fourier inversion we can prove that
1 (cid:90) π
B = ψ(θ)dθ,
0,1
π
0
2 (cid:90) π
B = cos(nθ)ψ(θ)dθ, for n ≥ 1, (2.6)
n,1
π
0
whereas for d ≥ 2, we have
(2n+d−1)[Γ((d−1)/2)]2 (cid:90) π
B = G(d−1)/2(cosθ)(sinθ)d−1ψ(θ)dθ, n ≥ 0, (2.7)
n,d 23−dπ Γ(d−1) n
0
where integration is taken componentwise. The matrices {B }∞ are called Schoenberg’s matrices.
n,d n=0
For the case m = 1, such result is reported by Gneiting (2013).
2.3 Scalar-valued positive definite functions on Sd ×Rk and the class Υ
d,k
Now, we consider the class of scalar-valued positive definite functions on Sd×Rk, for d,k ∈ N. Partic-
ularly, we focus on those continuous elements being geodesically isotropic in the spherical component
and stationary in the Euclidean one. Suppose that F is a member of the class P(Sd ×Rk) and there
exists a bounded and continuous mapping C : [0,π]×Rk → R such that
F((x,t),(y,s)) = C(θ(x,y),t−s), x,y ∈ Sd,t,s ∈ Rk. (2.8)
8
The continuous mappings C such that F in (2.8) is positive definite are the elements of the class Υ .
d,k
These classes are non-increasing in d and we have the inclusions
∞
(cid:92)
Υ := Υ ⊂ ··· ⊂ Υ ⊂ Υ .
∞,k d,k 2,k 1,k
d=1
The following theorem given by Berg and Porcu (2016) characterises completely the classes Υ and
d,k
Υ .
∞,k
Theorem 2.2. Let d and k be positive integers and C : [0,π]×Rk → R a continuous mapping, with
C(0,0) < ∞.
(1) The mapping C belongs to the class Υ if and only if it can be represented as
d,k
∞ (d−1)/2
C(θ,u) = (cid:88)ϕ˜ (u)Gn (cosθ), (θ,u) ∈ [0,π]×Rk, (2.9)
n,d
(d−1)/2
G (1)
n=0 n
where {ϕ˜ (·)}∞ is a sequence of elements in Φ , such that (cid:80)∞ ϕ˜ (0) < ∞.
n,d n=0 k,S n=0 n,d
(2) The mapping C belongs to the class Υ if and only if it admits the representation
∞,k
∞
(cid:88)
C(θ,u) = ϕ˜ (u)(cosθ)n, (θ,u) ∈ [0,π]×Rk, (2.10)
n
n=0
where {ϕ˜ (·)}∞ is a sequence of elements in Φ , with (cid:80)∞ ϕ˜ (0) < ∞.
n n=0 k,S n=0 n
Remark 2.1. The result given by Berg and Porcu (2016) is more general, since Rk can be replaced
for any locally compact group.
The orthogonality of Gegenbauer polynomials (Abramowitz and Stegun, 1970) implies that
1 (cid:90) π
ϕ˜ (u) = C(θ,u)dθ,
0,1
π
0
2 (cid:90) π
ϕ˜ (u) = cos(nθ)C(θ,u)dθ, for n ≥ 1, (2.11)
n,1
π
0
whereas for d ≥ 2,
(2n+d−1)[Γ((d−1)/2)]2 (cid:90) π
ϕ˜ (u) = G(d−1)/2(cosθ)(sinθ)d−1C(θ,u)dθ, n ≥ 0. (2.12)
n,d 23−dπ Γ(d−1) n
0
Berg and Porcu (2016) call the sequence of mappings {ϕ˜ (·)}∞ as Schoenberg’s functions.
n,d n=0
9
3 Space-time cross-covariances for RFs defined on spheres across
time
Let d, k and m be positive integers. We now study the class of matrix-valued positive definite
functions on Sd ×Rk, being bounded, continuous, geodesically isotropic in the spherical component
and stationary in the Euclidean one. The case k = 1 is particularly important, since Pm(Sd×R) can
be interpreted as the class of admissible space-time covariances for multivariate Gaussian RFs, with
spatial locations on the unit sphere. However, the main result of this section, given in Theorem 3.1 is
provided for an arbitrary k ∈ N, since we require this general framework in the coming sections, as it
will become apparent subsequently.
3.1 Matrix-valued positive definite functions on Sd × Rk: the class Υm and its
d,k
characterisation
Consider F ∈ Pm(Sd × Rk) and suppose that there exists a bounded and continuous mapping C :
[0,π]×Rk → Rm×m such that
F((x,t),(y,s)) = C(θ(x,y),t−s), x,y ∈ Sd,t,s ∈ Rk. (3.1)
Such mappings C are the elements of the class Υm . These classes are non-increasing in d and we
d,k
have the inclusions
∞
(cid:92)
Υm := Υm ⊂ ··· ⊂ Υm ⊂ Υm .
∞,k d,k 2,k 1,k
d=1
Now, we propose the generalisation of Theorem 2.2 to the multivariate case. Theorem 3.1 offers a
completecharacterisationoftheclassΥm andΥm ,foranym ≥ 1. Again,equalitiesandsummability
d,k ∞,k
conditions must be understood in a componentwise sense.
Theorem 3.1. Let d, k and m be positive integers and C : [0,π] × Rk → Rm×m a continuous
matrix-valued mapping, with C (0,0) < ∞, for all i = 1,...,m.
ii
(1) The mapping C belongs to the class Υm if and only if
d,k
∞ (d−1)/2
C(θ,u) = (cid:88)ϕ˜ (u)Gn (cosθ), (θ,u) ∈ [0,π]×Rk, (3.2)
n,d
(d−1)/2
G (1)
n=0 n
10