# «A GENERALIZED CONVOLUTION MODEL FOR MULTIVARIATE NONSTATIONARY SPATIAL PROCESSES Anandamayee Majumdar, Debashis Paul and Dianne Bautista Arizona ...»

Statistica Sinica 20 (2010), 675-695

## A GENERALIZED CONVOLUTION MODEL FOR

## MULTIVARIATE NONSTATIONARY SPATIAL PROCESSES

Anandamayee Majumdar, Debashis Paul and Dianne Bautista

Arizona State University, University of California, Davis

and Ohio State University

Abstract: We propose a ﬂexible class of nonstationary stochastic models for multivariate spatial data. The method is based on convolutions of spatially varying covariance kernels and produces mathematically valid covariance structures. This method generalizes the convolution approach suggested by Majumdar and Gelfand (2007) to extend multivariate spatial covariance functions to the nonstationary case.

A Bayesian method for estimation of the parameters in the covariance model based on a Gibbs sampler is proposed, then applied to simulated data. Model comparison is performed with the coregionalization model of Wackernagel (2003) that uses a stationary bivariate model. Based on posterior prediction results, the performance of our model appears to be considerably better.

Key words and phrases: Convolution, nonstationary process, posterior inference, predictive distribution, spatial statistics, spectral density.

1. Introduction Spatial modeling with ﬂexible classes of covariance functions has become a central topic of spatial statistics in recent years. One of the traditional approaches to modeling spatial stochastic processes is to consider parametric families of stationary processes, or processes that can be described through parametric classes of semi-variograms (Cressie (1993)). However, in spite of its simplicity, computational tractability, and interpretability, the stationarity assumption is often violated in practice, particularly when the data come from large, heterogeneous, regions. In various ﬁelds of applications, like soil science, environmental science, etc., it is often more reasonable to view the data as realizations of processes that only in a small neighborhood of a location behave like stationary processes. Also, it is often necessary to model two or more processes simultaneously and account for the possible correlation among various coordinate processes. For example, Majumdar and Gelfand (2007) consider an atmospheric pollution data consisting of 3 pollutants : CO, N O and N O2, whose concentrations in the atmosphere are correlated. A key question studied in this paper is modeling this correlation among the various coordinates while allowing for nonstationarity in space for the

## 676 ANANDAMAYEE MAJUMDAR, DEBASHIS PAUL AND DIANNE BAUTISTA

multivariate process. We propose a ﬂexible semiparametric model for multivariate nonstationary spatial processes. After reviewing the existing literature on nonstationary spatial modeling.A considerable amount of work over the last decade or so has focussed on modeling locally stationary processes (Fuentes (2002), Fuentes, Chen, Davis and Lackmann (2005), Gelfand, Schmidt, Banerjee and Sirmans (2004), Higdon (1997), Paciorek and Schervish (2006) and Nychka, Wikle, and Royle (2002)).

Dahlhaus (1996, 1997) gives a more formal treatment of locally stationary processes in the time series context in terms of evolutionary spectra of time series.

This research on the modeling of nonstationary processes might be thought of as the semi-parametric modeling of covariance functions. Higdon (2002) and Higdon, Swall, and Kern (1999) model the process as a convolution of a stationary process with a kernel of varying bandwidth. Thus, the observed process Y (s) ∫ is of the form Y (s) = Ks (x)Z(x)dx, where Z(x) is a stationary process, and the kernel Ks depends on the location s. Fuentes (2002) and Fuentes and Smith (2001) consider a convolution model in which the kernel has a ﬁxed bandwidth, while the process has a spatially varying parameter. Thus, ∫ K(s − x)Zθ(x) (s)dx, Y (s) = (1.1) D where {Zθ(x) (·) : x ∈ D} is a collection of independent stationary processes with covariance function parameterized by the function θ(·). Nychka, Wikle, and Royle (2002) consider a multiresolution analysis-based approach to model the spatial inhomogeneity that utilizes the smoothness of the process and its eﬀect on the covariances of the basis coeﬃcients, when the process is represented in a suitable wavelet-type basis.

One of the central themes of the various modeling schemes described above is that a process may be represented in the spectral domain locally as a superposition of Fourier frequencies with suitable (possibly spatially varying) weight functions. Recent work of Pintore and Holmes (2006) provides a solid mathematical foundation to this approach. Paciorek and Schervish (2006) derive an explicit representation for the covariance function for Higdon’s model when the kernel is multivariate Gaussian and use it to deﬁne a nonstationary version of the Mat´rn covariance function by utilizing the Gaussian scale mixture repree sentation of positive deﬁnite functions. Also, there are works on a diﬀerent type of nonstationary modeling through spatial deformations (see e.g., Sampson and Guttorp (1992)), but they do not concern us here.

The modeling approaches mentioned so far focus primarily on one-dimensional processes. In this paper, our main focus is on modeling nonstationary, multidimensional spatial processes. Existing approaches to modeling the multivariate

## GENERALIZED CONVOLUTION MODEL FOR SPATIAL PROCESSES 677

where ρj (·) are stationary covariance functions, and Tj are positive semideﬁnite matrices of rank 1. Christensen and Amemiya (2002) consider a diﬀerent class of multivariate processes that depend on a latent shifted-factor model structure.

Our work can be viewed as a generalization of the convolution model for correlated Gaussian processes proposed by Majumdar and Gelfand (2007). We extend their model to nonstationary settings. A key motivation is the assertion that when spatial inhomogeneity in the process is well-understood in terms of dependence on geographical locations, it makes sense to use that information directly in the speciﬁcation of the covariance kernel. For example, soil concentrations of Nitrogen, Carbon, and other nutrients and/or pollutants, that are spatially distributed, are relatively homogenous across similar land-use types (e.g., agricultural, urban, desert, transportation - and so on), but are non-homogeneous across spatial locations with diﬀerent land-use types. Usually the land-use types and their boundaries are clearly known (typically from satellite imagery). This is then an instance when nonstationary models are clearly advantageous compared to stationary models. Another example concerns land-values and diﬀerent economic indicators in a spatial area. Usually land-values are higher around (possibly multiple) business centers, and such information may be incorporated in the model as the known centers of the kernels at (3.1). It is also important for modeling multidimensional processes that the degree of correlations among the coordinate processes across diﬀerent spatial scales is allowed to vary. Keeping these goals in mind, we present a class of models that behave locally like stationary processes, but are globally nonstationary. The main contributions of this paper are: (i) speciﬁcation of the multivariate spatial cross-covariance function in terms of Fourier transforms of spatially varying spectra; (ii) incorporation of correlations among coordinate processes that vary with both frequency and location; (iii) derivation of precise mathematical conditions under which the process is nonsingular; and (iv) the provision for including local information about the process (e.g., smoothness, scale of variability, gradient of spatial correlation along a given direction) directly into the covariance model. The last goal is achieved by expressing the spatially varying coordinate spectra fj (s, ω) (as in (2.6)) as a sum of kernel-weighted stationary spectra, where the kernels have known shapes and diﬀerent (possibly pre-speciﬁed) centers, bandwidths and orientations. We also

## 678 ANANDAMAYEE MAJUMDAR, DEBASHIS PAUL AND DIANNE BAUTISTA

present a Bayesian estimation procedure based on Gibbs sampling for estimating a speciﬁc parametric covariance function and study its performance through simulation studies.The paper is organized as follows. We specify the model and discuss its properties in Section 2. In Section 3, we propose a special parametric subclass that is computationally easier to deal with. Also, we discuss various aspects of the model, such as parameter identiﬁability and the relation to some existing models, by focussing attention on a special bivariate model. In Section 4, we give an outline of a simulation study that illustrates the characteristics of the various processes generated by our model in the two-dimensional setting. In Section 5, we present a Bayesian estimation procedure and conduct a simulation study to demonstrate its eﬀectiveness. In Section 6, we discuss some related research directions. Some technical details and a detailed outline of the Gibbs sampling procedure for posterior inference are given in the supplementary material.

2. Construction of Covariances Through Convolution We consider a real-valued point-referenced univariate spatial process, Y (s), associated with locations s ∈ Rd. In this section, we construct a Gaussian spatial process model for an arbitrary ﬁnite set of locations in a region D ⊂ Rd by generalizing the construction of Majumdar and Gelfand (2007), and then extend it to whole of Rd.

and deﬁne A(ω) to be the N k × N k matrix with (l, l )th block All (ω), for 1 ≤ l, l ≤ k. Then A(ω) = F(ω)[(e(ω) R0 (ω)) ⊗ R(ω)]F(ω), where denotes Schur (or Hadamard) product, i.e., coordinate-wise product of two matrices of same dimension, and ⊗ denotes the Kronecker product.

Note that, for an arbitrary a ∈ Ck, a∗ (e(ω) R0 (ω))a = b∗ R0 (ω)b, where bl = al e−iω sl, l = 1,..., k. Therefore, if R0 (ω) is positive deﬁnite, then so is the T

2.2. Construction of nonstationary covariances on Rd We now generalize the construction of the nonstationary N × N covariance function C from the set {s1,..., sk } to the entire space Rd. Since a Gaussian process is determined entirely by its mean and covariance, given points s1,..., sk ∈ Rd, we can ﬁnd a zero mean Gaussian random vector (Yjl : 1 ≤ j ≤ N, 1 ≤ l ≤ k) with covariance matrix given by C. Moreover, this vector can be viewed as the realization of an N -dimensional random process Y(s) = (Y1 (s),..., YN (s)) at the points s1,..., sk, if we deﬁne Yjl = Yj (sl ). The next theorem states that an extension of the process Y(s) to arbitrary domains in Rd is possible.

## GENERALIZED CONVOLUTION MODEL FOR SPATIAL PROCESSES 681

The function fj (s, ω) can be interpreted as a location-dependent spectral density of a locally stationary stochastic process. If fj (s, ω) = fj (ω) for all j = 1,..., N, and ρ0 (s, s, ω) = 1, then C as in Theorem 2 becomes a covariance function of an N -dimensional stationary process on Rd.

2.3. Suﬃcient conditions for positive deﬁniteness In this subsection, we present a suﬃcient condition on the Fourier transforms of the cross-correlation functions, namely {ρjj }j=j, that guarantee the positivedeﬁniteness of the covariance function in the convolution model presented in Section 2.2 when the number of variables N is at most 4.

** Theorem 3. When N ≤ 4, suﬃcient conditions for positive deﬁniteness of R(ω) are the following.**

(i) 1 |ρjj |2 for all 1 ≤ j j ≤ N.

(ii) 1 |ρjl |2 + |ρlm |2 + |ρmj |2 − 2Re(ρjl ρlm ρmj ), for all 1 ≤ j l m ≤ N.

(iii) If 1 ≤ j = l = m = n ≤ N, then

Equality in place of any of the inequalities implies singularity of the matrix R(ω).

2.4. A general model A general formulation for the nonstationary covariance kernels comes from introducing some structure to the correlation function ρ0 (s, t, ω). One proposal

## 682 ANANDAMAYEE MAJUMDAR, DEBASHIS PAUL AND DIANNE BAUTISTA

and we are assuming that all the measurability conditions needed on the processes to deﬁne the stochastic integral are satisﬁed. The most manageable case from a practical point of view though, in our opinion, is when ρl (s, t) = ρ(s − t; θl ) for some parametric correlation function ρ(·; θ).

3.1. Speciﬁcation of the parametric spectral density and correlation We now give a complete description of a model that maintains a balance between ﬂexibility, and computational cost and interpretability. We choose ρ1 (·; τ ) to be an arbitrary parametric stationary correlation function on Rd with parameter τ. We assume that fj (ω; θjl ) is of the form cjl γ(ω; θjl ) for some scale parameter cjl 0 (note that we express θjl = (cjl, θjl )), and a parametric class of spectral densities γ(·; θ) that is closed under product. The latter means that given any m ≥ 1, there exists a function γ(·; ·) and functions cγ (· · · ), dγ (· · · ) of m variables such that, given parameters θ1,..., θm, ∏ m γ(·; θi ) = dγ (θ1, · · ·, θm )γ(·; cγ (θ1, · · ·, θm )).

i=1

where F −1 γ denotes the inverse Fourier transform of γ, i.e., the covariance function whose spectral density is γ. Also, for j = j, Gjj (s; θjl, θj l, νjj, κ) = cjl cjl dγ (θjl, θjl )(F −1 γ)(s; cγ (θjl, θjl )). (3.6)

3.3. Comparison with other nonstationary models Here we compare our model with other well-known models for nonstationary spatial covariances. For brevity, we focus on the univariate process speciﬁed by (1.1). Assuming that the set D is ﬁnite, say D = {x1,..., xM }, the covariance kernel for the process Y (·) is