# A flexible Clayton-like spatial copula with application to bounded support
data

Moreno Bevilacqua  Eloy Alvarado  Christian Caamaño-Carrillo  Facultad de
Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Viña del Mar, Chile
Dipartimento di Scienze Ambientali, Informatica e Statistica, Ca’ Foscari
University of Venice, Venice, Italy  Departamento de Industrias, Universidad
Técnica Federico Santa María, Santiago, Chile  Departamento de Estadística,
Universidad del Bío-Bío, Concepción, Chile

###### Abstract

The Gaussian copula is a powerful tool that has been widely used to model
spatial and/or temporal correlated data with arbitrary marginal distributions.
However, this kind of model can potentially be too restrictive since it
expresses a reflection symmetric dependence. In this paper, we propose a new
spatial copula model that makes it possible to obtain random fields with
arbitrary marginal distributions with a type of dependence that can be
reflection symmetric or not.

Particularly, we propose a new random field with uniform marginal
distributions that can be viewed as a spatial generalization of the classical
Clayton copula model. It is obtained through a power transformation of a
specific instance of a beta random field which in turn is obtained using a
transformation of two independent Gamma random fields.

For the proposed random field, we study the second-order properties and we
provide analytic expressions for the bivariate distribution and its
correlation. Finally, in the reflection symmetric case, we study the
associated geometrical properties. As an application of the proposed model we
focus on spatial modeling of data with bounded support. Specifically, we focus
on spatial regression models with marginal distribution of the beta type. In a
simulation study, we investigate the use of the weighted pairwise composite
likelihood method for the estimation of this model. Finally, the effectiveness
of our methodology is illustrated by analyzing point-referenced vegetation
index data using the Gaussian copula as benchmark. Our developments have been
implemented in an open-source package for the R statistical environment.

###### keywords:

Archimedean Copula, Beta random fields, Composite likelihood, Reflection
Asymmetry.

###### MSC:

[2020] Primary 62H11 , Secondary 62M30

††journal: Journal of Multivariate Analysis

##  1 Introduction

Many applications of statistics across a wide range of disciplines such as
climatology, environmental sciences and engineering, to name a few, show an
increasing interest in the statistical analysis of geo-referenced spatial
data. In order to model the inherent uncertainty of the data, Gaussian random
fields play a fundamental role; see [65, 2, 18], for instance. Indeed, the
Gaussian random fields can be completely specified in terms of mean and
correlation function and flexible correlation models such as the Matérn [65]
or the Generalized Wendland correlation models [10] can be specified to model
the geometrical properties of the random field.

Unfortunately, the Gaussian assumption is rarely met in practice. Indeed, in
many geostatistical applications, including climatology, oceanography, the
environment and the study of natural resources, the Gaussian framework is
unrealistic because the observed data have different features such as
asymmetry, heavy tails, positivity and bounded support.

The spatial generalized mixed models proposed by Diggle et al. [20] has been
widely used to model non-Gaussian spatial data [53]. Under this framework,
non-Gaussian models are specified using a specific link function and a latent
Gaussian random field through a conditional independence assumption. However,
the conditional independence assumption underlying models of this kind leads
to random fields with non-standard marginal distributions and with a “forced”
nugget effect that implies no mean square continuity [24, 53]. This can be a
potential problem when modeling data exhibiting continuity and/or
differentiability.

More flexible models can be obtained through a suitable transformation of a
Gaussian random field or independent copies thereof. This is an appealing
approach because the correlation of the transformed random field depends on
the correlation of the Gaussian random field and very often the geometrical
properties of the non-Gaussian random field are inherited from the Gaussian
random field. Notable examples can be found in [19] for log-Gaussian random
fields, [56] for Gaussian-log-Gaussian random fields, [73] for Skew-Gaussian
random fields, [72] for Tukey $g$-$h$ random fields, [7] for $t$ random
fields, [8] for Weibull random fields, [14] for sinh-arcsinh random fields,
and [53] for Poisson random fields just to mention a few.

A general powerful modeling tool to obtain non-Gaussian random fields with
arbitrary marginal distributions can be obtained under the copula framework
[37]. Unfortunately, the adaptation of the copula framework to the spatial
continuous setting is not a trivial task. Hereafter, we consider random fields
defined on $A$ a subset of ${\rm I\\!R}^{d}$ even if the proposed methodology
can be easily extended to more complicated spaces such as the continuous
space-time framework [27] or the spherical space [28, 59] or a linear network
[1]. A random field $\\{S(\bm{s}),\bm{s}\in A\\}$ with arbitrary marginal cdf
$F_{S}$ can be obtained through the transformation

$S(\bm{s})=F_{S}^{-1}\\{U(\bm{s})\\},$

where $U$ is a random field with uniform marginals and $F_{S}^{-1}$ is the
generalized inverse distribution function. This approach makes it possible to
model the marginal distributions and the dependence structure separately. In
particular the dependence is driven by the dependence of the random field $U$.

Inference and prediction of the non-Gaussian random field $S$ can be a
challenging task in particular when considering only one realization from $S$
which is the typical case for spatial data. It is important to highlight that
when considering independent realizations from $S$, the inference generally
simplifies and the set of copula models that can be used is larger; see for
instance [43] and the references therein.

However, in this paper we assume that only one realization of $S$ is available
which is the typical case for spatial data. In this case, a popular choice is
$U(\bm{s})=\Phi\\{Z(\bm{s})\\}$ where $\Phi$ is the cdf of the Gaussian
distribution and $\\{Z(\bm{s}),\bm{s}\in A\\}$ is a standard Gaussian random
field. This is the so-called Gaussian copula random field that has been widely
applied in the analysis of spatial and/or temporal correlated data; see [5,
38, 50, 32, 34, 19], just to mention a few. The main reason for success of the
Gaussian copula is that the dependence structure is indexed by a correlation
function which is helpful in many applications such as spatial statistics,
time series or longitudinal data.

However, Gaussian copula random fields suffer from some limitations due to the
kind of dependence structure they can model. In particular, the resulting type
of dependence is reflection symmetric, which basically means that high values
exhibit a spatial dependence similar to low ones. These restrictions are often
violated when analyzing Gaussian or non Gaussian real data. For instance
reflection asymmetry has been observed in rainfall or wind speed data [5, 67,
8]. To overcome these restrictions, a random field $U$ with uniform marginals
and more complex dependence structures than $\Phi\\{Z(\bm{s})\\}$ is needed.

A first flexible copula model for spatial data has been proposed by [5] and
deeply studied in [61] and [60]. The so-called Chi-squared copula has been
obtained by considering the cdf associated with a squared Gaussian random
field and the associated uniform random field takes the form
$U(\bm{s})=|1-2\Phi\\{Z(\bm{s})\\}|$. However this kind of model offers only a
specific type of bivariate dependence. As a more flexible alternative, a vine
copula approach can be performed [22, 32], but the computational time explodes
as the dimension increases because the number of possible configurations
becomes very high.

In this paper, we propose a new flexible random field with uniform marginals.
The main features of this model is that, unlike the Gaussian copula random
field, the type of dependence can be reflection symmetric or not. In addition,
as in the Gaussian copula or Chi-squared copula case, the dependence structure
of the proposed random field is indexed by a correlation function with an
additional parameter that characterizes the type of reflection (a)symmetry.

Specifically, we first define a general class of random fields with uniform
marginal distributions that we call Archimedean random field because it is
obtained generalizing the stochastic representation associated to the
Archimedean copula [25, 37]. Then we focus on a specific instance of the
Archimedean class that we call Clayton random field because it can be viewed
as a spatial generalization of the Clayton copula. It turns out that the
Clayton random field can be written as a power transformation of a specific
random field with beta marginals which is obtained through a transformation of
Gamma random fields obtained as sum of squared Gaussian random fields sharing
a common underlying correlation function.

For the Clayton random field, we study the second-order properties and we
provide analytic expressions for the bivariate pdf and cdf (the associated
bivariate copula density and copula function, respectively) and the
correlation function. It turns out that the bivariate pdf, the bivariate cdf
and the correlation function depend on some special functions, namely the
Gaussian hypergeometric function, the Appell function of the fourth type, and
the Kampé de Fériet function [31]. Finally, in a special case, we study the
associated geometrical properties. Particularly, we show that nice properties
such as mean-square continuity and degrees of mean-square differentiability
can be inherited from the underlying Gaussian random field.

As an application of the Clayton random field, we consider a spatial random
field with beta marginal distribution with a parametrization that allows to
specify a spatial mean regression model [23]. Estimating the (reparametrized)
beta random field or more generally random fields with arbitrary marginal
distributions obtained from the Clayton random field is not an easy task. In
fact, it must be said that the likelihood function associated with the
proposed Clayton random field involves an analytically intractable form. As a
consequence, likelihood-based estimation methods are unfeasible. Exploiting
the results on the bivariate distribution of the Clayton random field, we
investigate the use of the weighted pairwise composite likelihood [46, 68] as
method of estimation. In particular, in our simulation study, we study the
performance of the weighted pairwise composite likelihood when estimating the
parameters of the (reparametrized) beta random field. The methods proposed in
this paper are implemented in the R package GeoModels [13] and the R code for
reproducing the work is available as an online supplement.

The remainder of the article is organized as follows. In Section 2, we first
study an auxiliary random field with beta marginals which is a key tool when
defining the Clayton random field. In Section 3, we first define the
Archimedean random field and then we study a specific case, the Clayton random
field, and in particular we provide the associated bivariate pdf and cdf and
the correlation function. In addition we study geometrical properties of the
Clayton random field in the special reflection symmetric case. In Section 4,
we present a simulation study in order to investigate the performance of the
weighted pairwise composite likelihood method when estimating the mean
regression parameters and the dependence parameters of the (reparametrized)
beta random field obtained from the Clayton random field. In Section 5, we
illustrate the application of the proposed methodology by analyzing spatially
point-referenced vegetation indexes data using a random field with beta
marginal distribution. In particular we compare the proposed model using the
Gaussian copula model as a benchmark. All the proofs are deferred to the
Appendix.

##  2 A random field with beta marginals

Henceforth, given a weakly stationary random field $\\{Q(\bm{s}),\bm{s}\in
A\\}$ with $E\\{Q(\bm{s})\\}=\mu(\bm{s})$ and ${\rm
var}\left\\{Q(\bm{s})\right\\}=\sigma^{2}$, we denote by
$\rho_{Q}(\bm{h})={\rm corr}\\{Q(\bm{s}_{i}),Q(\bm{s}_{j})\\}$ its correlation
function, where $\bm{h}=\bm{s}_{i}-\bm{s}_{j}\in A$ is the lag separation
vector. For any set of distinct points
$(\bm{s}_{1},\ldots,\bm{s}_{n})^{\top}$, $n\in\mathbb{N}$, we denote by
$\bm{Q}_{ij}=(Q(\bm{s}_{i}),Q(\bm{s}_{j}))^{\top}$, $i\neq j$, the bivariate
random vector and by $\bm{Q}=(Q(\bm{s}_{1}),\ldots,Q(\bm{s}_{n}))^{\top}$ the
multivariate random vector. Moreover, we denote by $f_{Q}$ and $F_{Q}$ the
marginal probability density function (pdf) and cumulative distribution
function (cdf) of $Q(\bm{s})$, respectively, with $f_{\bm{Q}_{ij}}$
$F_{\bm{Q}_{ij}}$ the pdf (probability density function) and cdf (cumulative
distribution function) of $\bm{Q}_{ij}$ and with $f_{\bm{Q}}$ the pdf of
$\bm{Q}$ with $F_{\bm{Q}}$ the cdf of $\bm{Q}$ . Additionally, we denote by
$\Phi$ and $\phi$ the cdf and pdf of the standard Gaussian random variable,
respectively. Finally, if a random field has a $X$ marginal distribution, then
we call it $X$ random field.

Let $\\{Z(\bm{s}),\bm{s}\in A\\}$, a standard Gaussian random field with
correlation function $\rho_{Z}(\bm{h})$ and with some abuse of notation we set
$\rho(\bm{h})=\rho_{Z}(\bm{h})$. Henceforth, we call $Z$ and $\rho(\bm{h})$ as
the underlying Gaussian random field and correlation function respectively.

To build our random field with beta marginals we first consider a Gamma random
field $\\{G_{\psi}(\bm{s}),\bm{s}\in A\\}$ with marginal distribution
$G_{\psi}(\bm{s})\sim\mathcal{G}(\psi/2,1)$, defined as

$G_{\psi}(\bm{s})=\sum_{i=1}^{\psi}Z_{i}(\bm{s})^{2}/2,$

where $Z_{i}$, $i\in\\{1,\ldots,\psi\\}$ are mutually independent copies of
$Z$ with $\rho_{G_{\psi}}(\bm{h})=\rho^{2}(\bm{h})$. The associated
multivariate density $f_{\bm{G}_{\psi}}$ was discussed earlier by
Krishnamoorthy and Parthasarathy [42], and its properties have been studied
since then by several authors; see [41, 62]. It turns out that the analytical
expressions of the multivariate density can be derived only in some special
case. However, in the bivariate case, the pdf has general expression given by
[40]

$f_{\bm{G}_{\psi;ij}}(g_{i},g_{j})=\frac{(g_{i}g_{j})^{\psi/2-1}e^{-\frac{(g_{i}+g_{j})}{1-\rho^{2}(\bm{h})}}}{\Gamma\left(\psi/2\right)\\{1-\rho^{2}(\bm{h})\\}^{\nu/2}}\left[\frac{\sqrt{\rho^{2}(\bm{h})g_{i}g_{j}}}{1-\rho^{2}(\bm{h})}\right]^{1-\psi/2}I_{\psi/2-1}\left[\frac{2\sqrt{\rho^{2}(\bm{h})g_{i}g_{j}}}{1-\rho^{2}(\bm{h})}\right],$
(1)

where $I_{a}(x)$ is the modified Bessel function of the first kind of order a.
Now consider the random field $\\{Y_{\nu,\alpha}(\bm{s}),\bm{s}\in A\\}$
defined as:

$Y_{\nu,\alpha}(\bm{s})=\dfrac{H_{\nu}(\bm{s})}{H_{\nu}(\bm{s})+N_{\alpha}(\bm{s})}.$
(2)

where $H_{\nu}$, and $N_{\alpha}$ are two independent copies of the Gamma
random field $G_{\psi}$, $\psi\in\\{1,2,\ldots\\}$ sharing an underlying
common correlation function $\rho(\bm{h})$. By construction
$Y_{\nu,\alpha}(\bm{s})$ has beta marginal distribution
$\mathcal{B}(\nu/2,\alpha/2)$ with mean and variance given by ${\rm
E}\left\\{Y_{\nu,\alpha}(\bm{s})\right\\}=\nu/(\nu+\alpha)$ and ${\rm
var}\left\\{Y_{\nu,\alpha}(\bm{s})\right\\}=2\nu\alpha/\\{(\nu+\alpha)^{2}(\nu+\alpha+2)\\}$,
respectively. Hereafter, we call $Y_{\nu,\alpha}$ the auxiliary beta random
field with underlying correlation function $\rho(\bm{h})$.

A possible drawback for the gamma random fields used in (2) is that it is a
limited model due to the restrictions to the half-integers for the shape
parameter. Actually, in some special cases, it can assume any positive value
greater than zero. This feature is intimately related to the infinite
divisibility of the squared of the underlying Gaussian random fields $Z^{2}$
as shown in [41]. Characterization of the infinite divisibility of $Z^{2}$ has
been studied in [4, 21, 33, 70]. In particular, Bapat [4] provides a
characterization based on $\Omega_{n}$, the correlation matrix associated with
$\rho(\bm{h})$. Specifically, $\psi>0$ if and only if there exists a matrix
$S_{n}$ such that $S_{n}\Omega_{n}^{-1}S_{n}$ is an $M$-matrix [58], where
$S_{n}$ is a signature matrix, i.e., a diagonal matrix of size $n$ with
entries either $1$ or $-1$. This condition is satisfied, for instance, by a
stationary Gaussian random process $Z$ defined on $A=\mathbb{R}$ with an
exponential correlation function. However, it is not guaranteed to be
satisfied for any covariance matrix in any dimension.

These sorts of restrictions on the shape parameter of the Gamma random field
$G_{\nu}$ are clearly inherited by the auxiliary beta random field, which
means that $Y_{\nu,\alpha}$ is well defined for
$\nu,\alpha\in\\{1,2,\ldots\\}$. This certainly restricts the flexibility of
this kind of model.

In what follows, we make use of the generalized hypergeometric function
defined in Gradshteyn and Ryzhik [31] by:

${}_{p}F_{q}(a_{1},a_{2},\ldots,a_{p};b_{1},b_{2},\ldots,b_{q};x)=\sum\limits_{k=0}^{\infty}\frac{(a_{1})_{k},(a_{2})_{k},\ldots,(a_{p})_{k}}{(b_{1})_{k},(b_{2})_{k},\ldots,(b_{q})_{k}}\frac{x^{k}}{k!}\;\;\;\text{for}\;\;\;p,q\in\\{0,1,\ldots\\},$
(3)

where $(a)_{k}=\Gamma(a+k)/\Gamma(a)$, for $k\in\mathbb{N}\cup\\{0\\}$, is the
Pochhammer symbol. In particular we focus on the Gaussian hypergeometric
function

${}_{2}F_{1}(a,b,c;x)=\sum\limits_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}}{(c)_{k}}\frac{x^{k}}{k!},\;\;\;|x|<1,$
(4)

and we also consider the Appell hypergeometric function of the fourth type
[31] defined as

$F_{4}(a,b;c,c^{\prime};w,z)=\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{(a)_{k+m}(b)_{k+m}w^{k}z^{m}}{k!m!(c)_{k}(c^{\prime})_{m}},\;\;\;\;|\sqrt{w}|+|\sqrt{z}|<1,$
(5)

which is a special function of two variables. The special functions (4) and
(5) are related through the identity

$F_{4}(a,b;c,c^{\prime};w,z)=\sum\limits_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}z^{k}}{k!(c^{\prime})_{k}}{}_{2}F_{1}(a+k,b+k;c;w),\;\;\;\;|\sqrt{w}|+|\sqrt{z}|<1;$
(6)

see [15].

Equation (6) is useful for the computation of the $F_{4}$ function given that
efficient numerical computation of the ${}_{2}F_{1}$ function can be found in
different libraries of statistical softwares including R, MATLAB and Python.
Furthermore, we introduce a two-variable power series generalization called
the Kampé de Feriet function defined as

$F_{E;G;H}^{A;B;C}\left[\begin{array}[]{c}(a);(b_{B});(c_{C})\\\
(e_{E});(g_{G});(h_{H})\end{array}\middle|x,y\right]=\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\prod\limits_{j=1}^{A}(a_{j})_{k+m}\prod\limits_{j=1}^{B}(b_{j})_{k}\prod\limits_{j=1}^{C}(c_{j})_{m}x^{k}y^{m}}{k!m!\prod\limits_{j=1}^{E}(e_{j})_{k+m}\prod\limits_{j=1}^{G}(g_{j})_{k}\prod\limits_{j=1}^{R}(h_{j})_{m}}.$
(7)

This special function can be viewed as a generalization of the $F_{4}$ and
other types of Appell functions. For instance, it can be shown that
$F_{2;0;0}^{0;1;1}=F_{4}$ [64].

The following theorem provide the pdf of the bivariate distribution of the
auxiliary beta random field defined in (2) in terms of the Appell $F_{4}$
function.

######  Theorem 1.

Let $Y_{\nu,\alpha}$ the auxiliary beta random field defined in (2) with
underlying correlation function $\rho(\bm{h})$. Then the pdf of
$\bm{Y}_{\nu,\alpha;ij}$ is given by

$f_{\bm{Y}_{\nu,\alpha;ij}}(\bm{y}_{ij})=\frac{(y_{i}y_{j})^{\nu/2-1}\\{(1-y_{i})(1-y_{j})\\}^{\alpha/2-1}\Gamma^{2}\left(c\right)}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{-c}}F_{4}\left(c,c;\nu/2,\alpha/2;\rho^{2}(\bm{h})y_{i}y_{j},\rho^{2}(\bm{h})(1-y_{i})(1-y_{j})\right),$

with $c=(\nu+\alpha)/2$.

Note that $f_{\bm{Y}_{\nu,\alpha;ij}}$ is well defined for $\nu,\alpha>0$
irrespectively of the correlation function, as it is obtained from a bivariate
Gamma distribution. Moreover, when $\rho(\bm{h})=0$, it can be easily shown
that
$f_{\bm{Y}_{\nu,\alpha;ij}}(\bm{y}_{ij})=f_{Y_{\nu,\alpha}}(y_{i})f_{Y_{\nu,\alpha}}(y_{j})$.
This implies that zero pairwise correlation implies pairwise independence.

The following theorem provides the correlation function of the auxiliary beta
random field $Y_{\nu,\alpha}$ in terms of the Kampé de Fériet function defined
in (7).

######  Theorem 2.

Let $Y_{\nu,\alpha}$ the auxiliary beta random field defined in (2) with
underlying correlation $\rho(\bm{h})$. Then:

$\rho_{Y_{\nu,\alpha}}(\bm{h})=\frac{\nu(c+1)}{\alpha}\left[\\{1-\rho^{2}(\bm{h})\\}^{c}A-1\right],$
(8)

where $c=(\nu+\alpha)/2$ and

$A=F_{2;0;1}^{2;1;2}\left[\begin{array}[]{c}c;\frac{\alpha}{2};\frac{\nu}{2}+1\\\
c+1;-;\frac{\nu}{2}\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right].$

##  3 Clayton random fields

The goal of this section is to introduce what we call the Clayton random
field, which is a random field with standard uniform marginal distributions
with a flexible type of dependence allowing reflection (a)symmetry. We call
the proposed random field ‘Clayton’ because it can be viewed as a
generalization, of the classical Clayton copula.

To clarify the link between our construction and the classical Clayton copula,
let us first consider a random vector
${\mathbf{U}}=(U_{1},\dots,U_{n})^{\top}$, where

$\mathbf{U}=\varphi\left(\mathbf{E}/M\right).$ (9)

Here ${\bf E}=(E_{1},\dots,E_{n})^{\top}$ is a vector of i.i.d. random
variables with standard exponential distribution and $M$ an independent
positive random variable. The function $\varphi:{\rm I\\!R}^{+}\to[0,1]$ in
(9) is applied pointwise and it is the Laplace transform of the random
variable $M$, i.e., a completely monotone function [52]. Then, $\mathbf{U}$ is
a vector with uniform marginals $U_{i}$ [49, 47]. The cdf associated to the
random vector $\mathbf{U}$ has been widely studied and it belongs to a sub-
class of the so-called Archimedean copulas.

An $n$-dimensional Archimedean copula is defined with a generator function
$\varphi$ through:

$F_{\mathbf{U}}(\bm{u})=\varphi\\{\varphi^{-1}(u_{1})+\dots+\varphi^{-1}(u_{n})\\}$
(10)

and it has been shown by McNeil and Nešlehová [51] that a necessary and
sufficient condition for $\varphi$ to define a $n$-dimensional copula is that
$\varphi$ is a $n$-monotone function on $[0,\infty]$. For the precise
definition of $n$-monotonicity we refer the interested reader to [48, 51, 71].
If $D_{p}$ is the set of all $p$-monotone functions then it is known that
$D_{2}\supsetneq D_{3}\supsetneq\cdots\supsetneq D_{\infty}$, where
$D_{\infty}$ denotes the set of all completely monotone functions.

This implies that any function $\varphi\in D_{\infty}$ can be used to define
an Archimedean copula in arbitrary dimension $n\geq 2$ and it has been shown
in [49] that the random vector $\mathbf{U}$ in $(\ref{pppo})$ has cdf of type
$(\ref{pppo2})$ when $\varphi\in D_{\infty}$. For instance, choosing
$M\sim\mathcal{G}(1/a,1)$ with Laplace transform given by
$\psi(x)=(1+x)^{-1/a}$, $a>0$ with inverse $\psi^{-1}(y)=y^{-a}-1$, then the
Clayton copula

$F_{\mathbf{U}}(\bm{u})=\left(u_{1}^{-a}+\dots+u_{n}^{-a}-n+1\right)^{1/a},$

is obtained as a special case.

Inspired by the stochastic representation in (9) we now propose a new random
field with uniform marginals. Specifically, our idea is to obtain a random
field with uniform marginals generalizing the stochastic representation in (9)
to the spatial setting. To achieve this task, we first need to relax the
independence assumption of the exponential random variables. This basically
prevents the generation of a random field with “a forced” nugget effect. In
addition we need to replace the positive random variable $M$ with a positive
random field $\\{M(\bm{s}),\bm{s}\in A\\}$. This allows to avoid
identifiability problems (similar to the ones described in [26] and [7]) when
estimating with just one realization of the random field (the typical setting
for spatial data). Note that when independent replicates are available, our
construction in principle could be simplified assuming $M$ as a sequence of
independent random variables of gamma type.

Following these arguments, we generalize the model in (9) by defining the
Archimedean random field, i.e., a new class of random fields
$\\{U(\bm{s}),\bm{s}\in A\\}$ with uniform marginals as

$U(\bm{s})=\psi\left\\{E(\bm{s})/M(\bm{s})\right\\},$ (11)

where $E$ is a random field with standard exponential marginal distribution,
$M$ is a positive random field and $\psi$ is the Laplace transform of the
random variable $M(\bm{s})$.

In this work we focus on the special case where $E\equiv G_{2}$ and $M$ is a
copy of the Gamma random field defined in Section 2, that is $M\equiv G_{\nu}$
with marginal distribution $\mathcal{G}(\nu/2,1)$. Since the associated
Laplace transform is given by $\psi(t)=(1+t)^{-\nu/2}$, the random field in
(11) assumes the form

$U_{\nu}(\bm{s})=\left\\{\frac{G_{\nu}(\bm{s})}{G_{\nu}(\bm{s})+G_{2}(\bm{s})}\right\\}^{\nu/2},\quad\nu\in\\{1,2,\ldots\\}.$

It turns out that, under this specific choice, $U_{\nu}$ can be written as a
power transformation of a special case of the auxiliary beta random field
introduced in Section 2 that is

$U_{\nu}(\bm{s})=\\{Y_{\nu,2}(\bm{s})\\}^{\nu/2},\quad\nu\in\\{1,2,\ldots\\},$

where $Y_{\nu,2}$ has marginal distribution $\mathcal{B}(\nu/2,1)$. Since we
choose the Laplace transform of a Gamma marginal distribution in analogy with
the classical Clayton copula, we call $U_{\nu}$ a Clayton random field.

The multivariate pdf of $U_{\nu}$ is not explicitly known. However, the pdf of
the bivariate random vector
$\bm{U}_{\nu;ij}=(U_{\nu}(\bm{s}_{i}),U_{\nu}(\bm{s}_{j}))^{\top}$ (bivariate
copula density) can be easily obtained using Theorem 1, viz.

$\displaystyle
f_{\bm{U}_{\nu;ij}}(\bm{u}_{ij})=\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}F_{4}\left(\nu/2+1,\nu/2+1;\nu/2,1;\rho^{2}(\bm{h})(u_{i}u_{j})^{2/\nu},\rho^{2}(\bm{h})(1-u^{2/\nu}_{i})(1-u^{2/\nu}_{j})\right).$
(12)

If the bivariate random vector has uniform marginals, reflection symmetry can
be easily checked using the associated bivariate density through
$f_{\bm{U}_{\nu;ij}}(u_{i},u_{j})=f_{\bm{U}_{\nu;ij}}(1-u_{i},1-u_{j})$; see
[37]. Recalling that the $F_{4}$ function is symmetric, i.e.,
$F_{4}(\cdot,\cdot,\cdot,x,y)=F_{4}(\cdot,\cdot,\cdot,y,x)$, it is apparent
that $f_{\bm{U}_{\nu;ij}}(u_{j},u_{j})=f_{\bm{U}_{\nu;ij}}(1-u_{i},1-u_{j})$
only when $\nu=2$, that is reflection asymmetry is obtained when $\nu\neq 2$.

Note that the parameter $\nu$ only affects the bivariate distribution. That
is, it can be viewed as dependence parameter. In addition, the bivariate cdf
of $U_{\nu}$ (i.e., the associated bivariate copula function) is given in the
following theorem.

######  Theorem 3.

Let $U_{\nu}$ the Clayton random field in Equation (3) with underlying
correlation $\rho(\bm{h})$. Then the cdf of the bivariate random vector
$\bm{U}_{\nu;ij}=(U_{\nu}(\bm{s}_{i}),U_{\nu}(\bm{s}_{j}))^{\top}$ is given by

$\displaystyle F_{\bm{U}_{\nu;ij}}(\bm{t}_{ij})$ $\displaystyle=$
$\displaystyle
4t_{i}t_{j}\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(\nu/2\right)_{k}(t_{i}t_{j})^{2k/\nu}(\rho^{2}(\bm{h}))^{k+m}}{k!(2k+\nu)^{2}B^{2}\left(k+\nu/2,m+1\right)}{}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{i}\right)$
(13)
$\displaystyle\times\quad{}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{j}\right).$

From the copula point of view, if the dependence parameter
$\rho(\bm{h})\in[0,1)$ is equal to $0$, then it can be shown that independence
is achieved irrespective of $\nu$, that is
$F_{\bm{U}_{\nu;ij}}(\bm{t}_{ij})=F_{U_{i}}(t_{i})F_{U_{j}}(t_{j})=t_{i}t_{j}$.

We now compare the type of bivariate dependence induced by the Gaussian copula
random field $U^{G}(\bm{s})=\Phi\\{Z(\bm{s})\\}$ and the chi-squared copula
random field $U^{CS}(\bm{s})=|1-2\Phi\\{Z(\bm{s})\\}|$ [61] with the one
induced by the proposed Clayton random field
$U_{\nu}(\bm{s})=\\{Y_{\nu,2}(\bm{s})\\}^{\nu/2}$. Although copula theory uses
transformations to uniform marginals, it is better to consider transformation
to Gaussian marginals for identifying the type of dependence ([37], p. 9). In
this case, the plot of the bivariate densities associated with the proposed
model, the Gaussian copula and the chi-squared copula can be obtained by
computing

$f_{\SS_{ij}}(\bm{s}_{ij})=f_{A_{ij}}\\{F_{S}(\bm{s}_{i}),F_{S}(\bm{s}_{j})\\}f_{S}(\bm{s}_{i})f_{S}(\bm{s}_{j}),$
(14)

for $A_{ij}=\bm{U}_{\nu;ij},\bm{U}^{G}_{ij},\bm{U}^{C}_{ij}$ setting
$F_{S}=\Phi$ and $f_{S}=\phi$.

Fig. 1 depicts the contour plots of the three models increasing the underlying
correlation $\rho(\bm{h})\in\\{0.1,0.5,0.9\\}$ (from left to right). The red,
black and green colors are associated with the proposed model, the Gaussian
copula and chi-squared copula respectively. The cases $\nu\in\\{1,2,5\\}$
(first, second, and third row, respectively) are also considered for the
proposed model.

It turns out that when $\nu=2$, both the Gaussian and Clayton cases are
reflection symmetric irrespective of the correlation. However, the Clayton
model shows a non-elliptical dependence. This can be clearly appreciated when
the correlation is stronger while when the correlation is low the two
distributions tend to be similar, as expected. In fact, when the correlation
is zero we obtain independence in both cases.

If $\nu=1$ and $\nu=5$, then it can be appreciated that both the Clayton and
chi-squared models are not reflection symmetric, but our proposed methodology
is more flexible because it is able to model asymmetry to lower tail ($\nu=1$)
and asymmetry to upper tail ($\nu=5$), while the chi-squared copula is able to
model only asymmetry to lower tail; see [37], Sec. 2.14. Hence, our model is
more flexible when modeling reflection (a)symmetries and, as a consequence can
be applied to a wider range of spatial data. Note that other copula models can
capture both lower and upper tail asymmetry, such as mixture of Clayton and
reflected Clayton copulas [37]. However, when applying them in the spatial
context, they are not identifiable with just one realization of the random
field.

|  |
---|---|---
|  |
|  |
Fig. 1: Contour plots of the bivariate Gaussian density (black color), the
bivariate chi-squared copula density (green color) and of the bivariate
density (14) with $F_{S}=\Phi$ and $f_{S}=\phi$ (red color) with underlying
correlation $\rho(\bm{h})\in\\{0.1,0.5,0.9\\}$ (from left to right) when
$\nu=1$ (first row), $\nu=2$ (second row), $\nu=5$ (third row).

###  3.1 Correlation function and geometric properties

The following Theorem provides the correlation of the Clayton random field
$U_{\nu}$ in terms of the Kampé de Fériet function.

######  Theorem 4.

Let $U_{\nu}$ the Clayton random field in (3) with underlying correlation
$\rho(\bm{h})$. Then

$\rho_{U_{\nu}}(\bm{h})=3\left(\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}F_{2;1;0}^{2;2;1}\left[\begin{array}[]{c}\nu/2+1;\nu;1\\\
\nu+1;\frac{\nu}{2};-\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right]-1\right).$
(15)

Fig. 2 depicts the correlation function $\rho_{U_{\nu}}(\bm{h})$ varying the
asymmetry parameter $\nu\in\\{1,2,5\\}$ when the underlying correlation is the
Generalized Wendland model [10] defined as

${\cal
GW}_{\delta,\mu,b}(\bm{h})=\begin{cases}K\left\\{1-\left(||\bm{h}||/b\right)^{2}\right\\}^{\delta+\mu}{}_{2}F_{1}\left(\mu/2,(\mu+1)/2;\delta+\mu+1;1-\left(||\bm{h}||/b\right)^{2}\right),&0\leq||\bm{h}||\leq
b,\\\ 0,&||\bm{h}||>b,\end{cases}$

with
$K=\Gamma(\delta)\Gamma(2\delta+\mu+1)/\\{\Gamma(2\delta)\Gamma(\nu+\mu+1)2^{\mu+1}\\}$,
$\delta\geq 0$, $\mu\geq 0.5(d+1)+\delta$ and $b>0$ is the compact support.
This kind of correlation model is very flexible since it allows to parametrize
the mean squared differentiability of the underlying spatial field and it
includes the Matérn correlation as a special case model after a suitable
reparametrization [9]. Specifically, we set $\rho(\bm{h})={\cal
GW}_{0,4,0.15}(\bm{h})=\left(1-||\bm{h}||/0.15\right)^{4}_{+}$ and we plot
$\rho_{U_{\nu}}(\bm{h})$ versus the distance $||\bm{h}||$ varying the
parameter $\nu\in\\{1,2,5\\}$. It can be seen that $\rho_{U_{\nu}}(\bm{h})$
slightly increases when increasing $\nu$ for each $||\bm{h}||$.

---
Fig. 2: The correlation function $\rho_{U_{\nu}}(\bm{h})$ in (15) for
$\nu\in\\{1,2,5\\}$ when the underlying correlation model is
$\rho(\bm{h})={\cal GW}_{0,4,0.15}(\bm{h}).$

We now focus on the special case $\nu=2$ that is the reflection symmetric
case. The following corollary provide a closed form expression for the
correlation of the Clayton random field $U_{2}$.

######  Corollary 1.

Let $U_{2}$ the symmetric Clayton random field with underlying correlation
$\rho(\bm{h})$. Then

$\rho_{U_{2}}(\bm{h})=\frac{2[\rho^{2}(\bm{h})\\{3\rho^{2}(\bm{h})-1\\}-\\{\rho^{2}(\bm{h})-1\\}^{2}\ln\\{1-\rho^{2}(\bm{h})\\}]}{\rho^{4}(\bm{h})}-3.$
(16)

The simple form of the correlation function $\rho_{U_{2}}(\bm{h})$ allows to
study the geometrical properties of the Clayton random field. It turns out
that nice properties such as stationarity, mean-square continuity, degrees of
mean-square differentiability, long-range dependence and compact support can
be inherited from the underlying Gaussian random field $Z$ with correlation
$\rho(\bm{h})$. The following theorem resumes all these properties.

######  Theorem 5.

Let $U_{2}$ the symmetric Clayton random field with underlying Gaussian random
field $Z$ with correlation $\rho(\bm{h})$. Then,

  1. (i)

$U_{2}$ is weakly stationary;

  2. (ii)

$U_{2}$ is mean-square continuous if $Z$ is mean-square continuous;

  3. (iii)

$U_{2}$ is $0$-times mean-square differentiable if $Z$ is $0$-times mean-
square differentiable; for $m\geq 1$, if $\rho^{(1)}(\bm{0})=0$, $U_{\nu}$ is
$m$-times mean-square differentiable if $Z$ is $m$-times mean-square
differentiable;

  4. (iv)

$U_{2}$ is a long-range dependent process if $Z$ is a long-range dependent
process;

  5. (v)

$\rho(\bm{h})=0\implies\rho_{U}(\bm{h})=0$.

$\nu=1$ | $\nu=2$ | $\nu=5$
---|---|---
|  |
|  |
Fig. 3:  Six realizations of the Clayton random field $U_{\nu}$. First row:
using an underlying zero-time differentiable Generalized Wendland correlation
model ${\cal GW}_{0,4,0.2}$ and reflection asymmetry parameter equal to
$\nu\in\\{1,2,5\\}$ (from left to right). Second row: using an underlying one-
time differentiable Generalized Wendland ${\cal GW}_{1,4,0.2}$ correlation
model and asymmetry parameter equal to $\nu\in\\{1,2,5\\}$ (from left to
right).

The consequence of Theorem 5 is that, as in the Gaussian case, the geometrical
properties, the long range dependence and the compact support of the Clayton
random field can be modeled using a suitable correlation function. For
instance long-range dependence can be achieved using correlations with
specific features such as the Generalized Cauchy model [29, 45, 6]. In
addition, it can be shown that for the Generalized Wendland model ${\cal
GW}_{\delta,\mu,\beta}$ the condition $\rho^{(1)}(\bm{0})=0$ is satisfied when
$\delta\geq 1$; that is, $U_{2}$ inherits the mean square differentiability of
$Z$ if the underlying correlation model is Generalized Wendland.

Even though Theorem 5 is restricted to the case $\nu=2$ we conjecture that the
results are still valid for $U_{\nu}$, $\nu\neq 2$. Some empirical evidence
can be viewed looking at the sample paths of $U_{\nu}$. Fig. 3, first row,
depicts three realizations of $U_{\nu}$ with $\nu\in\\{1,2,5\\}$ (from left to
right) when the underlying correlation is $\varphi_{0,4,0.2}$ (zero time mean
square differentiable). The second row depicts three realizations of $U_{\nu}$
with $\nu\in\\{1,2,5\\}$ (from left to right) when the underlying correlation
is $\varphi_{1,4,0.2}$ (one time mean square differentiable).

It is important to stress that a not mean-square continuous Clayton random
field can be obtained by introducing a nugget effect, i.e., a discontinuity at
the origin of $\rho(\bm{h})$. This can be achieved by replacing the underlying
correlation function $\rho(\bm{h})$ with
$\rho^{*}(\bm{h})=\rho(\bm{h})(1-\tau^{2})+\tau^{2}\mathbb{1}_{0}(\bm{h})$,
where $0\leq\tau^{2}\leq 1$ represents the nugget effect.

Note that, in addition to the correlation function $\rho_{U_{\nu}}(\bm{h})$,
other measures of dependence can be considered. For instance it can be shown
that the Spearman correlation coefficient is equal to the correlation function
in this case. Similarly, using (13) it is possible to obtain the Kendall rank
correlation coefficient and the Blomqvist medial correlation coefficient
defined as [25, 54]. They are respectively given by

$\tau_{\bm{U}_{\nu}}(\bm{h})=4\left\\{\int_{[0,1]^{2}}F_{\bm{U}_{\nu}}(\mathbf{t}_{ij})dF_{\bm{U}_{\nu}}(\mathbf{t}_{ij})-\dfrac{1}{4}\right\\}$

and

$\beta_{\bm{U}_{\nu}}(\bm{h})=4\left\\{F_{\bm{U}_{\nu}}\left(\dfrac{1}{2},\dfrac{1}{2}\right)-\dfrac{1}{4}\right\\}.$

Finally, two common measures of extremal dependence of a bivariate
distribution are the lower and upper tail dependence coefficients [63, 37]
defined respectively as

$\lambda_{L}=\lim\limits_{t\to
0^{+}}\frac{F_{\bm{U}_{ij}}(t,t)}{t}\;\;\hbox{and}\;\;\lambda_{U}=\lim\limits_{t\to
1}\frac{1-2t+F_{\bm{U}_{ij}}(t,t)}{1-t},$

for a given bivariate copula $F_{\bm{U}_{ij}}(t_{i},t_{j})$. Finding closed
form expressions for $\lambda_{L}$ and $\lambda_{U}$ is challenging from a
mathematical point of view when considering the proposed copula
$F_{\bm{U}_{\nu;ij}}$, since it involves a double infinite sum of a squared
hypergeometric function. Based on some numerical evidence, we conjecture that
the proposed copula is tail-independent, that is, $\lambda_{L}=\lambda_{U}=0$.

##  4 Application to beta spatial regression model

As an application of the proposed Clayton random field, we focus on a random
field with marginal distributions of the (reparametrized) beta type that
allows to perform spatial mean regression for spatial data defined on a
bounded support. A random field $S$ with an arbitrary cdf $F_{S}$ using the
Clayton random field, can be obtained as

$S(\bm{s})=F_{S}^{-1}\\{U_{\nu}(\bm{s})\\},$

with bivariate distribution given by

$f_{\SS_{ij}}(\bm{s}_{ij})=f_{\bm{U}_{\nu;ij}}\\{F_{S}(\bm{s}_{i}),F_{S}(\bm{s}_{j})\\}f_{S}(\bm{s}_{i})f_{S}(\bm{s}_{j})$
(17)

and correlation function given by

$\rho_{S_{\nu}}(\bm{h})=\frac{\int\limits_{0}^{1}\int\limits_{0}^{1}F_{S}^{-1}(u_{i})F_{S}^{-1}(u_{j})f_{\bm{U}_{\nu;ij}}(\bm{u}_{ij})d\bm{u}_{ij}-{\rm
E}\left\\{S(\bm{s}_{i})\right\\}{\rm
E}\left\\{S(\bm{s}_{i})\right\\}}{\sqrt{{\rm
var}\left\\{S(\bm{s}_{i})\right\\}{\rm var}\left\\{S(\bm{s}_{j})\right\\}}}.$
(18)

Note that, if we assume continuous random variable then the generalized
inverse distribution function $F_{S}^{-1}$ is Lipschiz of order 1 [66] and,
using results in [3], this implies that mean square continuity and degrees of
mean square differentiability of $S$ are inherited from the uniform random
field, i.e., using Theorem 5, from the underlying Gaussian random field at
least for the case $\nu=2$.

As a special case, we focus on a random field
$\\{B_{\xi,\delta}(\bm{s}),\bm{s}\in A\\}$ with beta marginal distribution. It
can be defined using the Clayton random field, as follows:

$B_{\xi,\delta}(\bm{s})=F^{-1}_{B_{\xi,\delta}}\\{U_{\nu}(\bm{s})\\}.$ (19)

It is important to stress that the parameters of the beta random field (19),
unlike the beta random field defined in Section 2, are not restricted to be
positive integer values.

By construction $B_{\xi,\delta}(\bm{s})\sim\mathcal{B}(\xi,\delta)$ and the
marginal density and cdf are given respectively by

$f_{B_{\xi,\delta}}(y)=\frac{y^{\xi-1}(1-y)^{\delta-1}}{B(\xi,\delta)},\qquad
F_{B_{\xi,\delta}}(y)=I(y,\xi,\delta),\qquad 0<y<1,$

where $\xi>0,\delta>0$ are two shape parameters and $I(y,\xi,\delta)$ is the
regularized incomplete beta function defined as

$I(y,\xi,\delta)=\int_{0}^{y}t^{\xi-1}(1-t)^{\delta-1}dt=\frac{x^{\xi}}{\xi}{}_{2}F_{1}(\xi,1-\delta;\xi+1;y).$

The mean and variance are given respectively by ${\rm
E}\left\\{B_{\xi,\delta}(\bm{s})\right\\}=\xi/(\xi+\delta)$ and ${\rm
var}\left\\{B_{\xi,\delta}(\bm{s})\right\\}=2\xi\delta/\\{(\xi+\delta)^{2}(\xi+\delta+2)\\}$.

Here we consider the parametrization proposed in [23] that is the random field
$\\{B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta},\bm{s}\in A\\}$ with
$0<\mu(\bm{s})<1$. Using this parametrization, it can be shown that the
(spatially varying) mean of $B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}$ is
given by $\mu(\bm{s})$ that can be specified through a regression model with a
logistic link, i.e.,

$E\\{B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}\\}=\mu(\bm{s})=1/(1+e^{-X(\bm{s})^{\top}\bm{\beta}}),$

with $\bm{\beta}\subseteq{\rm I\\!R}^{k}$ a vector of regression parameters.
Finally the bivariate distribution and the correlation function can be
obtained using (17) and (18), respectively.

###  4.1 Nearest neighbors weighted pairwise composite likelihood estimation

Let us assume a random field with arbitrary marginals
$S=\\{S(\bm{s}),\bm{s}\in A\\}$ defined on $A\in\mathbb{R}^{d}$, with ${\rm
E}\left\\{S(\bm{s})\right\\}=\mu(\bm{s})$ and ${\rm
var}\\{S(\bm{s})\\}=\sigma^{2}$, with a parametric model for the underlying
correlation that is $\rho(\bm{h})=\rho(\bm{h},\bm{\vartheta})$, where
$\bm{\vartheta}$ is the vector of correlation parameters. Then the vector of
dependence parameters is given by $\bm{\theta}=(\nu,\bm{\vartheta})^{\top}$,
where $\nu$ is the reflection symmetry parameter.

For any set of distinct points $(\bm{s}_{1},\ldots,\bm{s}_{n})^{\top}$,
$\bm{s}_{i}\in A$, $n\in\mathbb{N}$, we denote by
$\SS_{ij}=(S(\bm{s}_{i}),S(\bm{s}_{j}))^{\top}$, $i\neq j$ and with
$f_{\SS_{ij}}$ we denote the associated probability density functions.

Following [46], the log-composite likelihood is an objective function defined
as a sum of $K$ sub-log-likelihoods, viz.

$CL(\bm{\theta})=\sum_{k=1}^{K}\ell(\bm{\theta};B_{k})w_{k},$

where $\bm{\theta}$ is the vector of unknown parameters, $B_{k}$ is a marginal
or conditional set of $\SS$ and $\ell(\bm{\theta};B_{k})$ is a log-likelihood
calculated by considering only the random variables in $B_{k}$ and $w_{k}$ are
suitable weights that do not depend on $\bm{\theta}$. The maximum CL estimate
is given by
$\hat{\bm{\theta}}=\textrm{argmax}_{\bm{\theta}}\,CL(\bm{\theta})$.

The weighted pairwise composite likelihood estimation method [11] is obtained
by setting $B_{k}=\SS_{ij}$ and in this case we obtain the pairwise log-
likelihood function
$\ell_{ij}(\bm{\theta})=\ln\\{f_{\SS_{ij}}(\bm{\theta})\\}$. The corresponding
weighted pairwise composite log-likelihoods function is given by

$wpl(\bm{\theta})=\sum_{i=1}^{n}\sum_{j\neq
i}^{n}\ell_{ij}(\bm{\theta})w_{ij}$

and $\hat{\bm{\theta}}=\textrm{argmax}_{\bm{\theta}}\,wpl(\bm{\theta})$ is the
associated estimator.

In general, a loss of statistical efficiency is expected for both cases with
respect to the maximum likelihood (ML) estimation and the role of the weights
$w_{ij}$ is to minimize this loss. Using theory of optimal estimating
equations [36], it can be easily seen [12] that the optimal weights require
the computation of the inverse of an $n(n-1)\times n(n-1)$ matrix which is
computationally even harder than what is required for ML estimation. Some
approximations of the optimal weights have been proposed in literature as, for
instance, in [44, 55]. However, the computation of these sorts of weights can
be computationally demanding for large $n$.

To avoid this computational problem, in this work, we adopt a weight function
based on nearest neighbors as proposed in [16]. Specifically, let
$N_{m}(\bm{s}_{l})$ the set of the neighbors of order $m\in\\{1,\ldots,n\\}$
of the point $\bm{s}_{l}\in A$. We make use of the following asymmetric weight
function for $i,j\in\\{1,\ldots,n\\}$ and $i\neq j$

$w_{ij}(m)=\begin{cases}1,&\bm{s}_{i}\in N_{m}(\bm{s}_{j}),\\\
0,&\text{otherwise}.\end{cases}$ (20)

Following the same arguments as in [11], it can be shown that, under
increasing domain asymptotic framework,
$\hat{\bm{\theta}}=\textrm{argmax}_{\bm{\theta}}\,wpl(\bm{\theta})$ is
consistent and asymptotically normal with asymptotic covariance matrix given
by the inverse of the Godambe information matrix defined as

$G(\bm{\theta})=H(\bm{\theta})J(\bm{\theta})^{-1}H(\bm{\theta})^{\intercal},$

where

$H(\bm{\theta})=-{\rm E}\left\\{\nabla^{2}wpl(\bm{\theta})\right\\},\qquad
J(\bm{\theta})={\rm E}\left\\{\nabla wpl(\bm{\theta})\nabla
wpl(\bm{\theta})\,^{\intercal}\right\\}.$

Standard error estimation can be obtained by considering the square root
diagonal elements of $G^{-1}(\widehat{\bm{\theta}})$. Moreover, model
selection can be performed by considering the information criterion, defined
as

$\mbox{PLIC}=-2wpl(\hat{\bm{\theta}})+2\mathrm{tr}(H(\hat{\bm{\theta}})G^{-1}(\hat{\bm{\theta}})),$

which is the composite likelihood version of the Akaike information criterion
(AIC) [69]. Note that the computation of standard errors and PLIC require
evaluation of the matrices $H(\hat{\bm{\theta}})$ and $J(\hat{\bm{\theta}})$.
However, the evaluation of $J(\hat{\bm{\theta}})$ is computationally
unfeasible for large datasets and in this case subsampling techniques can be
used to estimate $J(\bm{\theta})$ as in [12, 35]. A straightforward and more
robust alternative that we adopt in Section 5 is the parametric bootstrap
estimation of $G^{-1}_{a}(\bm{\theta})$.

###  4.2 A numerical example

We simulate, using Cholesky decomposition, $1000$ realizations of
$B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}$ observed at $N=1000$ spatial
location sites uniformly distributed in the unit square. In particular we
consider a reparametrized beta random field
$B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}$ obtained using the Clayton
random field $U_{\nu}$ setting $\nu\in\\{1,2,6\\}$. As underlying correlation
model we consider

${\cal GW}_{0,4,b}(\bm{h})=\left(1-||\bm{h}||/b\right)^{4}_{+}.$ (21)

The use of a compactly supported correlation function is a good strategy from
a computational viewpoint because in this case, when using the pairwise
composite likelihood method of estimation, the bivariate distribution
factorizes as the product of two marginal distributions when the underlying
correlation is zero. We consider two mean regression parameters, i.e., we
assume a spatially varying mean as

$\mu(\bm{s})=\dfrac{1}{1+\exp\left[-\\{\beta_{0}+\beta_{1}u_{1}(\bm{s})\\}\right]},$

with $\beta_{0}=0.2$ and $\beta_{1}=-0.2$, where $u_{1}(\bm{s})$ is a
$\mathcal{U}(0,1)$ random variable. Finally, we set the shape parameters as
$\delta\in\\{0.5,1.5,2.5\\}$ and the compact support parameter $b=0.2$. The
estimation is performed using the weighted pairwise composite likelihood
method using the weight function (20) setting $m=2$.

It is important to stress that, in principle, the parameter $\nu$ can be
estimated because the bivariate distributions in (12) and (17) are defined for
each $\nu>0$. However, recall that the existence of the random field $U_{\nu}$
is guaranteed only for $\nu\in\\{1,2,\ldots\\}$. As a consequence, in this
numerical example and in the real data applications the reflection asymmetry
parameter $\nu$ is assumed fixed. An alternative strategy is a two-step
estimation. In the first step, all the parameters including the asymmetry
parameter $\nu$ are estimated. In the second step $\nu$ is fixed equal to the
rounded value of the estimate obtained at first step and then the marginal and
correlation dependence parameters are estimated.

Table 1 depicts Bias and MSE for each parameter of the nine scenarios. Some
general patterns can be observed. For instance the MSE of $\hat{b}$ is not
affected by the choice of $\delta$ and $\nu$ and the MSE of $\hat{\delta}$ and
of the estimated regression parameters increases when increasing $\delta$,
irrespective of $\nu$. Finally, as a general comment, the estimates are
overall unbiased and symmetric distributed. As an example, Fig. 4 shows the
centered boxplots of the weighted pairwise likelihood estimates for the mean
regression parameters $\beta_{0}$, $\beta_{1}$, the shape parameter when
$\delta=1.5$ and the spatial compact support parameter $b$ for different
$\nu\in\\{1,2,6\\}$.

Table 1: Bias and MSE of each parameter for different setting of $\nu$ and
$\delta$ when estimating the RF $B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}$
with marginal beta distribution using weighted pairwise composite likelihood.
The underlying correlation function is ${\cal
GW}_{0,4,b}(\bm{h})=\left(1-||\bm{h}||/b\right)^{4}_{+}$ with $b=0.2$ and the
mean is specified as
$\mu(\bm{s})=(1+\exp\left[-\\{\beta_{0}+\beta_{1}u_{1}(\bm{s})\\}\right])^{-1}$
with $\beta_{0}=0.2$ and $\beta_{1}=-0.2$.

$\nu$ | $1$ | $2$ | $6$
---|---|---|---
$\delta$ | 0.5 | 1.5 | 2.5 | 0.5 | 1.5 | 2.5 | 0.5 | 1.5 | 2.5
| Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE | Bias | MSE
$\hat{\beta}_{0}$ | 0.0091 | 0.0140 | 0.0091 | 0.0103 | 0.0083 | 0.0080 | -0.0050 | 0.0138 | -0.0031 | 0.0104 | -0.0028 | 0.0080 | -0.0057 | 0.0169 | -0.0030 | 0.0129 | -0.0030 | 0.0099
$\hat{\beta}_{1}$ | -0.0070 | 0.0310 | -0.0066 | 0.0220 | -0.0058 | 0.0167 | 0.0093 | 0.0287 | 0.0051 | 0.0223 | 0.0047 | 0.0171 | -0.0025 | 0.0321 | -0.0074 | 0.0235 | -0.0067 | 0.0179
$\hat{b}$ | -0.0245 | 0.0014 | -0.0247 | 0.0014 | -0.0247 | 0.0014 | -0.0017 | 0.0007 | -0.0025 | 0.0007 | -0.0026 | 0.0007 | -0.0155 | 0.0009 | -0.0156 | 0.0009 | -0.0156 | 0.0009
$\hat{\delta}$ | 0.0034 | 0.0006 | 0.0087 | 0.0068 | 0.0152 | 0.0209 | 0.0039 | 0.0006 | 0.0139 | 0.0070 | 0.0243 | 0.0217 | 0.0014 | 0.0006 | 0.0036 | 0.0069 | 0.0058 | 0.0214

---
Fig. 4: From left to right: centered boxplots of
$\hat{\beta}_{0},\hat{\beta}_{1}$ (first row), $\hat{b}$ and $\hat{\delta}$
(second row) for $\nu\in\\{1,2,6\\}$, respectively, when $\delta=1.5$ when
estimating the RF $B_{\mu(\bm{s})\delta,(1-\mu(\bm{s}))\delta}$ with marginal
beta distribution using weighted pairwise composite likelihood. The underlying
correlation function is ${\cal
GW}_{0,4,b}(\bm{h})=\left(1-||\bm{h}||/b\right)^{4}_{+}$ with $b=0.2$ and the
mean is specified as
$\mu(\bm{s})=(1+\exp\left[-\\{\beta_{0}+\beta_{1}u_{1}(\bm{s})\\}\right])^{-1}$
with $\beta_{0}=0.2$ and $\beta_{1}=-0.2$.

##  5 Application to Vegetation Indexes

In the study of ecological phenomena and dynamics, the use of remote sensing
measurements and associated transformations is a common practice. One
particularly stands out from the rest, the Normalized Difference Vegetation
Index (NDVI); as reviewed in [57], the use of the NDVI gained traction in
recent years. This index is derived from the red near-infrared reflectance
ratio defined as

$NDVI=\dfrac{NIR-RED}{NIR+RED},$

where $NIR$ and $RED$ are the amounts of near-infrared and red light,
respectively, reflected by the vegetation and measured by satellites. This
index ranges from $-1$ to $+1$, where negatives values represent an absence of
vegetation.

To illustrate our proposed methodology, we make use of the NDVI extracted from
Google Earth Engine [30] using data provided from The National Oceanic and
Atmospheric Administration (NOAA) Climate Data Record (CDR) from May 1st,
2018. We analyze our data using R and the GeoModels package [13] in which our
methodology was implemented. In this example, we focus our study on 1000
randomly sampled NDVI data between longitude $-92$ and $-97$, and latitude
between $33$ and $37$, corresponding roughly to the east side of Texas, US.
Fig. 5 show a coloured map after a suitable projection, the normalized
histogram, and the empirical semi-variogram of NDVI data.

In our analysis, we consider two random fields sharing the same marginal
distribution. In particular we consider the (reparametrized) beta random field
(assuming a constant mean) obtained from the Clayton random field and from the
Gaussian copula random field that is we consider:

$B^{C}_{\mu\delta,(1-\mu)\delta}(\bm{s})=F^{-1}_{B_{\mu\delta,(1-\mu)\delta}}\\{U_{\nu}(\bm{s})\\},\qquad
B^{G}_{\mu\delta,(1-\mu)\delta}(\bm{s})=F^{-1}_{B_{\mu\delta,(1-\mu)\delta}}\left[\Phi\\{Z(\bm{s})\\}\right],$

where $\mu=1/(1+e^{-\beta})$ and $\delta>0$. Note that NDVI data are bounded
supported from $-1$ to $+1$. However, the beta random field which is defined
on $(0,1)$ can be easily extended to an arbitrary bounded support
$(a_{1},a_{2})$ through the transformation
$(B_{\mu\delta,(1-\mu)\delta}-a_{1})/(a_{2}-a_{1})$ with $a_{2}>a_{1}$.

We first examine the suitability of the beta marginal model for the NDVI data.
In order to check the reliability of our marginal distribution, we compare the
histogram with the fitted beta distribution, where the marginal parameters
$\beta$ and $\delta$ have been estimated through maximum likelihood assuming
independence ($\hat{\beta}=0.8721,\hat{\delta}=355.52$). From Fig. 5b, it can
be appreciated that our marginal model seems to be a viable hypothesis. Note
that, even in the NDVI data potentially ranges from $-1$ to $1$, our data are
basically concentrated between $0.1$ and $0.7$ approximatively.

In addition, the semivariogram depicted in Fig. 5c shows a clear spatial
dependence and a negligible nugget effect in the data.

|  |
---|---|---
(a) | (b) | (c)
Fig. 5: From left to right: a) coloured map of the NDVI data; b) normalized
histogram of the NDVI data compared with the fitted density of the beta
distribution; c) empirical semivariogram of the NDVI data.

To identify the presence of reflection asymmetry, we apply the normal score
transformation to the NDVI data and we compute the scatterplot. Fig. 6 depicts
the spatial scatterplot for four different neighborhood orders, which clearly
show the presence of reflection asymmetry in the data; in particular,
asymmetry to upper tail can be detected.

|
---|---
|
Fig. 6: Spatial scatterplots of NDVI data for four increasing neighborhood
orders (from left to right) after normal score transformation.

We estimate both beta random fields $B^{C}_{\mu\delta,(1-\mu)\delta}$ and
$B^{G}_{\mu\delta,(1-\mu)\delta}$ with the weighted pairwise composite
likelihood estimation method using, as in the simulation study, the weight
function (20) setting $m=2$. We assume that the underlying isotropic
correlation function is of the generalized Wendland type as defined in (21).

For the proposed $B^{C}_{\mu\delta,(1-\mu)\delta}$ model, the (a)symmetry
parameter is fixed equal to $\nu\in\\{2,4,6\\}$. We consider these specific
values based on the empirical evidence of upper tail asymmetry (see Fig. 6)
and taking into account the discussion on Fig. 1 in Section 3.

Table 2 depicts the composite likelihood estimates with associated standard
errors and, in addition, the composite likelihood AIC (PLIC) as defined in
(4.1) for both types of models. They are computed using parametric bootstrap
as explained in Section 4.2. It turns out that the estimation of the marginal
parameters $\beta$ and $\delta$ are quite similar, as expected. However, the
compact support parameter estimate, i.e., the spatial dependence, is clearly
greater for the proposed model; see also Fig. 7. More importantly, the PLIC
criterion selects the $B^{C}_{\mu\delta,(1-\mu)\delta}$ model with $\nu=4$.
Note that, also in the reflection symmetric case ($\nu=2$) the PLIC criterion
selects the $B^{C}_{\mu\delta,(1-\mu)\delta}$ model with respect to the
$B^{G}_{\mu\delta,(1-\mu)\delta}$ model.

Table 2: Estimates of the beta random fields $B^{C}_{\mu\delta,(1-\mu)\delta}$
(the proposed model) for $\nu\in\\{2,4,6\\}$ and
$B^{G}_{\mu\delta,(1-\mu)\delta}$ (the Gaussian copula model) with
$\mu=1/(1+e^{-\beta})$ when the underlying correlation model is ${\cal
GW}_{0,4,b}(\bm{h})=\left(1-||\bm{h}||/b\right)^{4}_{+}$. In parenthesis, for
each parameter, standard error estimation is reported. Last line shows the
composite likelihood AIC information criterion (PLIC).

Parameter | $B^{G}_{\mu\delta,(1-\mu)\delta}$ | $B^{C}_{\mu\delta,(1-\mu)\delta}$
---|---|---
$\nu=2$ | $\nu=4$ | $\nu=6$
$\hat{\beta}$ |  | 0.8733
---
(0.0127)
| 0.8727
---
(0.0238)
| 0.8715
---
(0.0245)
| 0.8713
---
(0.0245)
$\hat{b}$ |  | 112.4101
---
(12.1615)
| 473.6463
---
(95.6832)
| 433.3901
---
(81.6638)
| 415.3676
---
(77.1154)
$\hat{\delta}$ |  | 355.9033
---
(35.7105)
| 331.6079
---
(41.5803)
| 330.5262
---
(46.4087)
| 329.7499
---
(48.4774)
PLIC | -14349 | -14372 | -14380 | -14368

Finally, in Fig. 7 we provide a graphical comparison between the empirical and
fitted semivariogram of the $B^{G}_{\mu\delta,(1-\mu)\delta}$ model (red line)
and the $B^{C}_{\mu\delta,(1-\mu)\delta}$ model (black line) when $\nu=4$. In
can be seen that the proposed model is able to better capture the spatial
dependence when compared with the Gaussian copula case.

---
Fig. 7: Empirical and estimated semivariogram of NDVI of
$B^{C}_{\mu\delta,(1-\mu)\delta}$ (the proposed model, black line) when
$\nu=4$ and $B^{G}_{\mu\delta,(1-\mu)\delta}$ (the Gaussian copula, red line).

##  6 Concluding remarks

We have introduced a new kind of spatial random field with uniform marginal
distribution as a viable alternative to the spatial Gaussian copula. The main
feature of our model is that, unlike the Gaussian copula case, the type of
dependence can be reflection symmetric or not. In particular the dependence
structure of the proposed random field is indexed by a correlation function
with an additional parameter that characterizes the type of reflection
(a)symmetry and, as a consequence, it can be useful for modeling Gaussian and
non-Gaussian spatial data exhibiting reflection asymmetries. To the best of
our knowledge, this is the first copula that can be used for spatial data that
allows both left and right bivariate asymmetries

As an application of the proposed model we considered a random field with beta
marginals and we applied it to the analysis of Normalized Difference
Vegetation Indexes data. In particular we showed that the beta random field
obtained using the proposed Clayton random field performs better than our
benchmark (the Gaussian copula) in terms of Composite Likelihood Akaike
Information Criterion, providing a practical example of the feasibility of our
results.

For simplicity of presentation, we restrict the treatment to the spatial
Euclidean setting. However, the results presented in this paper can be applied
to more general spaces. The key is the specification of a well-defined
correlation function in a specific space for the underlying Gaussian random
field such as correlation models defined on the space-time setting or on the
sphere [28, 59].

A possible limitation of the proposed model is the lack of amenable
expressions of the associated multivariate distributions, particularly due to
the natural occurrence of involved special functions, such as the Kampé de
Fériet function. This prevents an inference approach based on the full
likelihood and the computation of the optimal predictor. In the first case,
our simulation study shows that an inferential approach based on weighted
pairwise composite likelihood using the bivariate density given in (12) could
be an effective solution for estimating the unknown parameters. In the second
case, even if the prediction has not been investigated in this paper, it can
be performed using the results presented in this work coupled with the methods
proposed for instance in [39].

Finally, the general class proposed in (11) potentially opens the door to new
copula models for spatial data with complex dependencies. In this paper we
have focused in a specific instance of the class. However other types of
complex copula models for spatial data can be developed starting from (11) and
this topic will be subject to future research.

## Acknowledgments

Moreno Bevilacqua acknowledges financial support from grant FONDECYT 1200068
and ANID/PIA/ANILLOS ACT210096 and ANID project Data Observatory Foundation
DO210001 from the Chilean government and project MATH-AMSUD 22-MATH-06
(AMSUD220041). The work of Christian Caamaño-Carrillo was partially supported
by grant FONDECYT 11220066 from the Chilean government and DIUBB 2120538 IF/R
from the University of Bío-Bío. The work of Eloy Alvarado was partially
supported by CONICYT-PFCHA/DOCTORADO-BECAS-CHILE/2018-21180953.

## Appendix

### Proof of Theorem 1

For simplicity of notation we use $\bm{W}=H_{\nu}(\bm{s})$ and
$\bm{R}=N_{\alpha}(\bm{s})$. Under the transformation
$g_{i}=y_{i}r_{i}(1-y_{i})^{-1}$ and $g_{j}=y_{j}r_{j}(1-y_{j})^{-1}$ in (1)
with Jacobian
$J((g_{i},g_{j})\to(y_{i},y_{j}))=r_{i}r_{j}\\{(1-y_{i})(1-y_{j})\\}^{-2}$ and
using series expansion of the hypergeometric function ${}_{0}F_{1}$, we have

$\displaystyle f_{\bm{Y}_{\nu,\alpha;ij}}(\bm{y}_{ij})$
$\displaystyle=\int\limits_{\mathbb{R}_{+}^{2}}f_{\bm{W}_{ij}|\bm{R}_{ij}}(\bm{w}_{ij}|\bm{r}_{ij})f_{\bm{R}_{ij}}(\bm{r}_{ij})Jd\bm{r}_{ij}=\frac{(y_{i}y_{j})^{\nu/2-1}\\{(1-y_{i})(1-y_{j})\\}^{-\nu/2-1}}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)(1-\rho^{2}(\bm{h}))^{(\nu+\alpha)/2}}\int\limits_{\mathbb{R}_{+}^{2}}(r_{i}r_{j})^{(\nu+\alpha)/2-1}e^{-\frac{1}{1-\rho^{2}(\bm{h})}\left(\frac{r_{i}y_{i}}{1-y_{i}}+\frac{r_{j}y_{j}}{1-y_{j}}\right)}e^{-\frac{r_{i}+r_{j}}{1-\rho^{2}(\bm{h})}}$
$\displaystyle\quad\times\quad{}_{0}F_{1}\left(\nu/2;\frac{\rho^{2}(\bm{h})y_{i}y_{j}r_{i}r_{j}}{(1-y_{i})(1-y_{j})\\{1-\rho^{2}(\bm{h})\\}^{2}}\right){}_{0}F_{1}\left(\alpha/2;\frac{\rho^{2}(\bm{h})r_{i}r_{j}}{\\{1-\rho^{2}(\bm{h})\\}^{2}}\right)d\bm{r}_{ij}$
$\displaystyle=\frac{(y_{i}y_{j})^{\nu/2-1}\\{(1-y_{i})(1-y_{j})\\}^{-\nu/2-1}}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{(\nu+\alpha)/2}}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{I(k,m)}{k!m!\left(\nu/2\right)_{k}\left(\alpha/2\right)_{m}}\left\\{\frac{\rho^{2}(\bm{h})y_{i}y_{j}}{(1-y_{i})(1-y_{j})\\{1-\rho^{2}(\bm{h})\\}^{2}}\right\\}^{k}\left\\{\frac{\rho^{2}(\bm{h})}{\\{1-\rho^{2}(\bm{h})\\}^{2}}\right\\}^{m},$
(22)

where, using Fubini’s Theorem and (3.381.4) of [31], we obtain

$\displaystyle I(k,m)$
$\displaystyle=\int\limits_{\mathbb{R}_{+}}r_{i}^{(\nu+\alpha)/2+k+m-1}e^{-\left\\{\frac{1}{(1-\rho^{2}(\bm{h}))(1-y_{i})}\right\\}r_{i}}dr_{i}\int\limits_{\mathbb{R}_{+}}r_{j}^{(\nu+\alpha)/2+k+m-1}e^{-\left\\{\frac{1}{(1-\rho^{2}(\bm{h}))(1-y_{j})}\right\\}r_{j}}dr_{j}$
$\displaystyle=\Gamma^{2}\left((\nu+\alpha)/2+k+m\right)\left\\{(1-\rho^{2}(\bm{h}))(1-y_{i})\right\\}^{(\nu+\alpha)/2+k+m}\left\\{(1-\rho^{2}(\bm{h}))(1-y_{j})\right\\}^{(\nu+\alpha)/2+k+m}$
(23)

and combining (Proof of Theorem 1) and (Proof of Theorem 1), we obtain

$\displaystyle f_{\bm{Y}_{\nu,\alpha;ij}}(\bm{y}_{ij})=$
$\displaystyle\frac{(y_{i}y_{j})^{\nu/2-1}\\{(1-y_{i})(1-y_{j})\\}^{\alpha/2-1}\Gamma^{2}\left((\nu+\alpha)/2\right)}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{(\nu+\alpha)/2}}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left((\nu+\alpha)/2\right)_{k+m}^{2}}{k!m!\left(\nu/2\right)_{k}\left(\alpha/2\right)_{m}}\left\\{\rho^{2}(\bm{h})y_{i}y_{j}\right\\}^{k}\left\\{\rho^{2}(\bm{h})(1-y_{i})(1-y_{j})\right\\}^{m}$
$\displaystyle=$
$\displaystyle\frac{(y_{i}y_{j})^{\nu/2-1}\\{(1-y_{i})(1-y_{j})\\}^{\alpha/2-1}\Gamma^{2}\left((\nu+\alpha)/2\right)}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{-(\nu+\alpha)/2}}F_{4}\left((\nu+\alpha)/2,(\nu+\alpha)/2;\nu/2,\alpha/2;\rho^{2}(\bm{h})y_{i}y_{j},\rho^{2}(\bm{h})(1-y_{i})(1-y_{j})\right),$

where $F_{4}(a,b;c,c^{\prime};w,z)$ is the fourth Appell hypergeometric
function of two variables defined in (5).

### Proof of Theorem 2

We first cosider the $(a,a)$-th product moment ${\rm
E}\left\\{Y_{\nu,\alpha}^{a}(\bm{s}_{i})Y_{\nu,\alpha}^{a}(\bm{s}_{j})\right\\}$
of the beta random field $Y_{\nu,\alpha}$ with underlying correlation
$\rho(\bm{h})$. Using $c=(\nu+\alpha)/2$, then by definition

$\displaystyle{\rm
E}\left\\{Y_{\nu,\alpha}^{a}(\bm{s}_{i})Y_{\nu,\alpha}^{a}(\bm{s}_{j})\right\\}$
$\displaystyle=\frac{\Gamma^{2}\left(c\right)}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{-(\nu+\alpha)/2}}\int\limits_{[0,1]^{2}}y_{i}^{\nu/2+a-1}y_{j}^{\nu/2+a-1}\\{(1-y_{i})(1-y_{j})\\}^{\alpha/2-1}$
$\displaystyle\quad\times\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(c\right)_{k+m}^{2}}{k!m!\left(\nu/2\right)_{k}\left(\alpha/2\right)_{m}}\left\\{\rho^{2}(\bm{h})y_{i}y_{j}\right\\}^{k}\left\\{\rho^{2}(\bm{h})(1-y_{i})(1-y_{j})\right\\}^{m}d\bm{y}_{ij}$
$\displaystyle=\frac{\Gamma^{2}\left(c\right)}{\Gamma^{2}\left(\nu/2\right)\Gamma^{2}\left(\alpha/2\right)\\{1-\rho^{2}(\bm{h})\\}^{-c}}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(c\right)_{k+m}^{2}}{k!m!\left(\nu/2\right)_{k}\left(\alpha/2\right)_{m}}I(k,m)\rho^{2}(\bm{h})^{k+m},$
(24)

where, using Fubini’s Theorem,

$I(k,m)=\int\limits_{0}^{1}y_{i}^{\nu/2+a+k-1}(1-y_{i})^{\alpha/2+m-1}dy_{i}\int\limits_{0}^{1}y_{j}^{\nu/2+a+k-1}(1-y_{j})^{\alpha/2+m-1}dy_{j}.$

Then using (3.251) of [31], we obtain

$I(k,m)=\dfrac{\Gamma\left(\nu/2+a+k\right)\Gamma^{2}\left(\alpha/2+m\right)\Gamma\left(\nu/2+a+k\right)}{\Gamma\left(c+a+m+k\right)\Gamma\left(c+a+m+k\right)}$
(25)

and combining (25) with (Proof of Theorem 2), we obtain

${\rm
E}\left\\{Y_{\nu,\alpha}^{a}(\bm{s}_{i})Y_{\nu,\alpha}^{a}(\bm{s}_{j})\right\\}=\frac{\Gamma^{2}\left(c\right)\Gamma\left(\nu/2+a\right)^{2}\\{1-\rho^{2}(\bm{h})\\}^{c}}{\Gamma^{2}\left(\nu/2\right)\Gamma\left(c+a\right)^{2}}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(c\right)_{k+m}^{2}\left(\nu/2+a\right)^{2}_{k}\left(\alpha/2\right)_{m}^{2}\\{\rho^{2}(\bm{h})\\}^{k+m}}{k!m!\left(\nu/2\right)_{k}\left(\alpha/2\right)_{m}\left(c+a\right)^{2}_{k+m}}.$

The latter equation can be reduced to

${\rm
E}\left\\{Y_{\nu,\alpha}^{a}(\bm{s}_{i})Y_{\nu,\alpha}^{a}(\bm{s}_{j})\right\\}=\frac{\Gamma^{2}\left(c\right)\Gamma\left(\nu/2+a\right)^{2}\\{1-\rho^{2}(\bm{h})\\}^{c}}{\Gamma^{2}\left(\nu/2\right)\Gamma\left(c+a\right)^{2}}F_{2;1;0}^{2;2;1}\left[\begin{array}[]{c}c;\frac{\nu}{2}+a;\frac{\alpha}{2}\\\
c+a;\frac{\nu}{2};-\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right],$
(26)

where $F_{E;G;H}^{A;B;C}$ is the Kampé de Fériet function defined in (7).

Using (26) with $a=1$ and the fact that ${\rm
E}\left\\{Y_{\nu,\alpha}(\bm{s})\right\\}=\nu/(\nu+\alpha)$ and ${\rm
var}\left\\{Y_{\nu,\alpha}(\bm{s})\right\\}=2\nu\alpha/\\{(\nu+\alpha)^{2}(\nu+\alpha+2)\\}$,
the correlation can be written as

$\rho_{Y_{\nu,\alpha}}(\bm{h})=\frac{\nu(c+1)}{\alpha}\left[\\{1-\rho^{2}(\bm{h})\\}^{c}A-1\right],$

where

$A=F_{2;1;0}^{2;2;1}\left[\begin{array}[]{c}c;\frac{\nu}{2}+1;\frac{\alpha}{2}\\\
c+1;\frac{\nu}{2};-\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right]=F_{2;0;1}^{2;1;2}\left[\begin{array}[]{c}c;\frac{\alpha}{2};\frac{\nu}{2}+1\\\
c+1;-;\frac{\nu}{2}\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right].$

### Proof of Theorem 3

If $U_{\nu}$ is a uniform random field with underlying correlation function
$\rho(\bm{h})$, then the bivariate pdf is given by

$f_{\bm{U}_{\nu;ij}}(\bm{u}_{ij})=\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}F_{4}\left(\nu/2+1,\nu/2+1;\nu/2,1;\rho^{2}(\bm{h})(u_{i}u_{j})^{2/\nu},\rho^{2}(\bm{h})(1-u_{i}^{2/\nu})(1-u_{j}^{2/\nu})\right).$

This implies that the bivariate cdf is given by

$\displaystyle F_{\bm{U}_{\nu;ij}}(\bm{t}_{ij})$
$\displaystyle=\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}\int\limits_{0}^{t_{j}}\int\limits_{0}^{t_{i}}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(\nu/2+1\right)_{k+m}^{2}}{k!m!\left(\nu/2\right)_{k}\left(1\right)_{m}}\left\\{\rho^{2}(\bm{h})(u_{i}u_{j})^{2/\nu}\right\\}^{k}\left\\{\rho^{2}(\bm{h})(1-u_{i}^{2/\nu})(1-u_{j}^{2/\nu})\right\\}^{m}d\bm{t}_{ij}$
$\displaystyle=\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(\nu/2+1\right)_{k+m}^{2}}{k!m!\left(\nu/2\right)_{k}\left(1\right)_{m}}I(k,m)\\{\rho^{2}(\bm{h})\\}^{k+m},$
(27)

where using Fubini’s Theorem with the change of variable $x=u^{2/\nu}$ and the
definition $(8.391)$ of the incomplete beta function in [31], we obtain

$\displaystyle I(k,m)$
$\displaystyle=\int\limits_{0}^{t_{i}}u_{i}^{2k/\nu}(1-u_{i}^{2/\nu})^{m}du_{i}\int\limits_{0}^{t_{j}}u_{j}^{{}^{2k/\nu}}(1-u_{j}^{2/\nu})^{m}du_{j}$
$\displaystyle=\frac{\nu^{2}(t_{i}t_{j})^{2k/\nu+1}}{(2k+\nu)^{2}}{}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{i}\right){}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{j}\right)$
(28)

and combining (Proof of Theorem 3) and (Proof of Theorem 3), we obtain

$\displaystyle F_{\bm{U}_{\nu;ij}}(\bm{t}_{ij})$
$\displaystyle=\nu^{2}t_{i}t_{j}\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}\sum\limits_{k=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{\left(\nu/2+1\right)_{k+m}^{2}(t_{i}t_{j})^{2k/\nu}\\{\rho^{2}(\bm{h})\\}^{k+m}}{k!(2k+\nu)^{2}\left(\nu/2\right)_{k}(1)^{2}_{m}}$
$\displaystyle\quad\times\quad{}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{i}\right){}_{2}F_{1}\left(k+\nu/2,-m;k+\nu/2+1;t^{2/\nu}_{j}\right).$

Finally, simplifying the Pochhammer values we obtain (13).

### Proof of Theorem 4

Since ${\rm E}\left\\{U(\bm{s}_{i})U(\bm{s}_{j})\right\\}={\rm
E}\left\\{Y_{\nu,\alpha}^{a}(\bm{s}_{i})Y_{\nu,\alpha}^{a}(\bm{s}_{j})\right\\}$
and setting $\alpha=2$ and $a=\nu/2$ in (26) we obtain:

${\rm
E}\left\\{U(\bm{s}_{i})U(\bm{s}_{j})\right\\}=\frac{\\{1-\rho^{2}(\bm{h})\\}^{\nu/2+1}}{4}F_{2;1;0}^{2;2;1}\left[\begin{array}[]{c}\nu/2+1;\nu;1\\\
\nu+1;\frac{\nu}{2};-\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right].$

Then since ${\rm E}\left\\{U(\bm{s})\right\\}=0.5$ and ${\rm
var}\left\\{U(\bm{s})\right\\}=1/12$, we obtain (15):

### Proof of corollary 1

Replacing $\nu=2$ in equation (15) we get

$\displaystyle\rho_{U}(\bm{h})$
$\displaystyle=3\left[\\{1-\rho^{2}(\bm{h})\\}^{2}F_{2;1;0}^{2;2;1}\left[\begin{array}[]{c}2;2;1\\\
3;1;-\end{array}\middle|\rho^{2}(\bm{h}),\rho^{2}(\bm{h})\right]-1\right]=3\left[\\{1-\rho^{2}(\bm{h})\\}^{2}\sum_{k=0}^{\infty}\sum_{m=0}^{\infty}\dfrac{(2)_{k+m}^{2}(2)_{k}^{2}(1)_{m}}{k!m!(3)_{k+m}^{2}(1)_{k}}\\{\rho^{2}(\bm{h})\\}^{k+m}-1\right]$
$\displaystyle=3\left[\\{1-\rho^{2}(\bm{h})\\}^{2}\left\\{\sum_{k=0}^{\infty}\dfrac{4\rho^{2k}(\bm{h})}{(2+k)^{2}}\sum_{m=1}^{k+1}m^{2}\right\\}-1\right]=3\left[\\{1-\rho^{2}(\bm{h})\\}^{2}\left\\{\dfrac{2}{3}\left(\sum_{k=0}^{\infty}2k\rho^{2k}(\bm{h})+\sum_{k=0}^{\infty}\dfrac{\rho^{2k}(\bm{h})}{k+2}+\sum_{k=0}^{\infty}\rho^{2k}(\bm{h})\right)\right\\}-1\right]$
$\displaystyle=3\left[\\{1-\rho^{2}(\bm{h})\\}^{2}\left\\{\dfrac{-2\\{\rho^{2}(\bm{h})-3\rho^{4}(\bm{h})+(\rho^{2}(\bm{h})-1)^{2}\ln(1-\rho^{2}(\bm{h}))\\}}{3\rho^{4}(\bm{h})(\rho^{2}(\bm{h})-1)^{2}}\right\\}-1\right]$
$\displaystyle=\dfrac{2[\rho^{2}(\bm{h})\\{3\rho^{2}(\bm{h})-1\\}-\\{\rho^{2}(\bm{h})-1\\}^{2}\ln\\{1-\rho^{2}(\bm{h})\\}]}{\rho^{4}(\bm{h})}-3.$

For the double summation we make use of (1.4) in [17], the Lerch’s trascendent
defined as $\Phi(z,s,a)=\sum_{n=0}^{\infty}(a+n)^{-s}z^{n}$ and its link to
the polylogarithm function given by $\Phi(z,s,2)=z^{-2}(Li_{s}(z)-z).$ Since
$Li_{1}(\rho^{2}(\bm{h}))=-\ln(1-\rho^{2}(\bm{h}))$, in our case, the term

$\sum\limits_{k=0}^{\infty}\dfrac{\rho^{2k}(\bm{h})}{k+2}=\Phi(\rho^{2}(\bm{h}),1,2)$

simplifies to

$\dfrac{-\rho^{2}(\bm{h})-\ln\\{1-\rho^{2}(\bm{h})\\}}{\rho^{4}(\bm{h})}.$

### Proof of Theorem 5

If the underlying Gaussian random field $Z$ is a weakly stationary process
with correlation $\rho(\bm{h})$ then from (8) it is straightforward to see
that $U$ is also weakly stationary.Following [65], the mean-square continuity
and $m$-times mean-square differentiability of $U$ are equivalent to the
continuity and $2m$-times differentiability of $\rho_{U}(\bm{h})$ at
$\bm{h}=0$. From (16), using simple limits properties it can be easily seen
that $\lim_{\bm{h}\rightarrow 0}\rho_{U}(\bm{h})=1$ if $\rho(\bm{h})=1$.
Hence, $U$ is mean-square continuous if $Z$ is mean-square continuous.
For the mean square differentiability, let the underlying Gaussian random
field $Z$ $m$-times differentiable. Then it can be shown that
$\rho_{U}^{(2m)}(\bm{h})\big{|}_{\bm{h}\bm{=0}}<\infty$, hence, $U$ is
$m$-times times mean square differentiable. For instance, assume that the
random field $Z$ is $1$-times differentiable that is
$\rho^{(i)}(\bm{h})\big{|}_{\bm{h}\bm{=0}}<\infty$ for $i\in\\{1,2\\}$. The
first derivative is given by

$\displaystyle\rho_{U}^{(1)}(\bm{h})$ $\displaystyle=$
$\displaystyle\frac{-4\rho^{(1)}(\bm{h})}{\rho^{5}(\bm{h})}\left[-2\ln\\{1-\rho^{2}(\bm{h})\\}+2\rho^{2}(\bm{h})(\ln\\{1-\rho^{2}(\bm{h})\\}-1)+\rho^{4}(\bm{h})\right].$

Computing the limit at $\bm{h}\rightarrow 0$ and using L’Hôspital’s rule
conveniently, we obtain
$\rho_{U}^{(1)}(\bm{h})\big{|}_{\bm{h}\bm{=0}}=4\rho^{(1)}(\bm{h})$. The
second derivative is given by

$\displaystyle\rho_{U}^{(2)}(\bm{h})$ $\displaystyle=$
$\displaystyle\frac{4}{\rho^{6}(\bm{h})}\Bigg{[}(\rho^{(1)}(\bm{h}))^{2}\left[-10\ln\\{1-\rho^{2}(\bm{h})\\}+2\rho^{2}(\bm{h})(3\ln\\{1-\rho^{2}(\bm{h})\\}-5)+\rho^{4}(\bm{h})\right]$
$\displaystyle-\rho(\bm{h})\rho^{(2)}(\bm{h})\left[-2\ln\\{1-\rho^{2}(\bm{h})\\}+2\rho^{2}(\bm{h})\left[\ln\\{1-\rho^{2}(\bm{h})\\}-1\right]+\rho^{4}(\bm{h})\right]\Bigg{]}.$

Taking the limit at $\bm{h}\rightarrow 0$ with L’Hôspital’s rule and assuming
$\rho^{\prime}(\bm{0})=0$, it can be shown that
$\rho_{U}^{(2)}(\bm{h})\big{|}_{\bm{h}\bm{=0}}<\infty$ and hence the uniform
random field is $1-$times differentiable. Similarly for $m>1$, the
$(2m)$-derivatives of $\rho_{U}(\bm{h})$ evaluated at $\bm{h}=\bm{0}$ takes
finite values if $\rho^{\prime}(\bm{0})=0$.
Following [45], the process $U$ is long-range dependent if the correlation of
$U$ is such that
$\int_{\mathbb{R}_{+}^{n}}|\rho_{U}(\bm{h})|d^{n}\bm{h}=\infty$. Using the
expression of $\rho_{U}(\bm{h})$ inside integral and integrating term by term
it is sufficient to show that one of the integral diverges. Applying the
inequality

$\int_{\mathbb{R}_{+}^{n}}\rho(\bm{h})d^{n}\bm{h}\leq\int_{\mathbb{R}_{+}^{n}}\frac{1}{\rho^{2}(\bm{h})}d^{n}\bm{h}$

for $\rho(\bm{h})\in[0,1]$ to one of the term shows that
$\int_{\mathbb{R}_{+}^{n}}|\rho_{U}(\bm{h})|d^{n}\bm{h}=\infty$ if and only if
$\int_{\mathbb{R}_{+}^{n}}|\rho(\bm{h})|d^{n}\bm{h}=\infty$. Hence, the random
field $U$ has long-range dependence if $Z$ has long-range dependence.

Finally, to prove that $\rho(\bm{h})=0\Rightarrow\rho_{U}(\bm{h})=0$, it is
sufficient to apply L’Hôspital’s rule two times and show that:

$\displaystyle\lim_{\rho(\bm{h})\rightarrow 0}\rho_{U}(\bm{h})=0.$

## References

  * Anderes et al. [2020] E. Anderes, J. Møller, J. G. Rasmussen, Isotropic covariance functions on graphs and their edges, The Annals of Statistics 48 (2020) 2478 – 2503.
  * Banerjee et al. [2004] S. Banerjee, B. P. Carlin, A. E. Gelfand, Hierarchical Modeling and Analysis for Spatial Data, Chapman & Hall/CRC Press, Boca Raton: FL, 2004.
  * Banerjee and Gelfand [2003] S. Banerjee, A. Gelfand, On smoothness properties of spatial processes, Journal of Multivariate Analysis 84 (2003) 85–100.
  * Bapat [1989] R. B. Bapat, Infinite divisibility of multivariate gamma distributions and m-matrices, Sankhy$\bar{a}$ A 51 (1989) 73–78.
  * Bárdossy [2006] A. Bárdossy, Copula-based geostatistical models for groundwater quality parameters, Water Resources Research 42 (2006).
  * Berg et al. [2008] C. Berg, J. Mateu, E. Porcu, The dagum family of isotropic correlation functions, Bernoulli 14 (2008) 1134–1149.
  * Bevilacqua et al. [2021] M. Bevilacqua, C. Caamaño-Carrillo, R. B. Arellano-Valle, V. Morales-Oñate, Non-gaussian geostatistical modeling using (skew) t processes, Scandinavian Journal of Statistics 48 (2021) 212–245.
  * Bevilacqua et al. [2020] M. Bevilacqua, C. Caamaño-Carrillo, C. Gaetan, On modeling positive continuous data with spatiotemporal dependence, Environmetrics 31 (2020) e2632.
  * Bevilacqua et al. [2022] M. Bevilacqua, C. Caamaño-Carrillo, E. Porcu, Unifying compactly supported and matérn covariance functions in spatial statistics, Journal of Multivariate Analysis 189 (2022) 104949.
  * Bevilacqua et al. [2019] M. Bevilacqua, T. Faouzi, R. Furrer, E. Porcu, Estimation and prediction using generalized Wendland covariance functions under fixed domain asymptotics, The Annals of Statistics 47 (2019) 828 – 856.
  * Bevilacqua and Gaetan [2015] M. Bevilacqua, C. Gaetan, Comparing composite likelihood methods based on pairs for spatial gaussian random fields, Statistics and Computing 25 (2015) 877–892.
  * Bevilacqua et al. [2012] M. Bevilacqua, C. Gaetan, J. Mateu, E. Porcu, Estimating space and space-time covariance functions for large data sets: A weighted composite likelihood approach, Journal of the American Statistical Association 107 (2012) 268–280.
  * Bevilacqua et al. [2023] M. Bevilacqua, V. Morales-Oñate, C. Caamaño-Carrillo, GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis, 2023\. R package version 1.0.7.
  * Blasi et al. [2022] F. Blasi, C. Caamaño-Carrillo, M. Bevilacqua, R. Furrer, A selective view of climatological data and likelihood estimation, Spatial Statistics 50 (2022) 100596. Special Issue: The Impact of Spatial Statistics.
  * Brychkov and Saad [2017] Y. A. Brychkov, N. Saad, On some formulas for the appell function $f_{2}(a,b,b^{\prime};c,c^{\prime};w;z)$, Journal Integral Transforms and Special Functions 25 (2017) 1465 –1483.
  * Caamaño-Carrillo et al. [2023] C. Caamaño-Carrillo, M. Bevilacqua, C. López, V. Morales-Oñate, Nearest neighbours weighted composite likelihood based on pairs for (non-)Gaussian massive spatial data with an application to Tukey-hh random fields estimation, 2023. Research Square.
  * Choi [2003] J.-S. Choi, Notes on formal manipulations of double series, Communications of the Korean Mathematical Society 18 (2003) 781–789.
  * Cressie and Wikle [2011] N. Cressie, C. Wikle, Statistics for Spatio-Temporal Data., Wiley Series in Probability and Statistics. Wiley, 2011\.
  * De Oliveira [2006] V. De Oliveira, On optimal point and block prediction in log-gaussian random fields, Scandinavian Journal of Statistics 33 (2006) 523–540.
  * Diggle et al. [1998] P. Diggle, J. Tawn, R. Moyeed, Model-based geostatistics, Journal of the Royal Statistical Society: Series C (Applied Statistics) 47 (1998) 299–350.
  * Eisenbaum and Kaspi [2006] N. Eisenbaum, H. Kaspi, A characterization of the infinitely divisible squared gaussian processes, The Annals of Probability (2006) 728–742.
  * Erhardt et al. [2015] T. M. Erhardt, C. Czado, U. Schepsmeier, Spatial composite likelihood inference using local c-vines, Journal of Multivariate Analysis 138 (2015) 74–88. High-Dimensional Dependence and Copulas.
  * Ferrari and Cribari-Neto [2004] S. Ferrari, F. Cribari-Neto, Beta regression for modelling rates and proportions, Journal of Applied Statistics 31 (2004) 799–815.
  * Gelfand and Schliep [2016] A. E. Gelfand, E. M. Schliep, Spatial statistics and gaussian processes: A beautiful marriage, Spatial Statistics 18 (2016) 86–104. Spatial Statistics Avignon: Emerging Patterns.
  * Genest and MacKay [1986] C. Genest, R. J. MacKay, Copules archimédiennes et familles de lois bidimensionnelles dont les marges sont données, Canadian Journal of Statistics 14 (1986) 145–159.
  * Genton and Zhang [2012] M. G. Genton, H. Zhang, Identifiability problems in some non-Gaussian spatial random fields, Chilean Journal of Statistics 3 (2012).
  * Gneiting [2002] T. Gneiting, Stationary covariance functions for space-time data, Journal of the American Statistical Association 97 (2002) 590–600.
  * Gneiting [2013] T. Gneiting, Strictly and non-strictly positive definite functions on spheres, Bernoulli 19 (2013) 1327–1349.
  * Gneiting and Schlather [2004] T. Gneiting, M. Schlather, Stochastic models that separate fractal dimension and the hurst effect, SIAM Rev. 46 (2004) 269–282.
  * Gorelick et al. [2017] N. Gorelick, M. Hancher, M. Dixon, S. Ilyushchenko, D. Thau, R. Moore, Google earth engine: Planetary-scale geospatial analysis for everyone, Remote Sensing of Environment (2017).
  * Gradshteyn and Ryzhik [2007] I. Gradshteyn, I. Ryzhik, Table of Integrals, Series, and Products, Academic Press, New York, 7 edition, 2007.
  * Gräler [2014] B. Gräler, Modelling skewed spatial random fields through the spatial vine copula, Spatial Statistics 10 (2014) 87–102.
  * Griffiths [1970] R. C. Griffiths, Infinitely divisible multivariate gamma distributions, Sankhy$\bar{a}$ Ser. A 32 (1970) 393–404.
  * Guolo and Varin [2014] A. Guolo, C. Varin, Beta regression for time series analysis of bounded data, with application to Canada Google® flu trends, The Annals of Applied Statistics 8 (2014) 74–88.
  * Heagerty and Lele [1998] P. Heagerty, S. Lele, A composite likelihood approach to binary spatial data, Journal of the American Statistical Association 93 (1998) 1099 –1111.
  * Heyde [1997] C. Heyde, Quasi-Likelihood and Its Application: A General Approach to Optimal Parameter Estimation, Springer, New York, 1997\.
  * Joe [2014] H. Joe, Dependence modeling with copulas, Chapman and Hall/CRC, Boca Raton, FL, 2014\.
  * Kazianka and Pilz [2010a] H. Kazianka, J. Pilz, Copula-based geostatistical modeling of continuous and discrete data including covariates, Stochastic Environmental Research and Risk Assessment 24 (2010a) 661–673.
  * Kazianka and Pilz [2010b] H. Kazianka, J. Pilz, Spatial interpolation using copula-based geostatistical models, Springer Netherlands, Dordrecht, 2010b, pp. 307–319.
  * Kibble [1941] W. Kibble, A two-variate gamma type distribution, Sankhyā: The Indian Journal of Statistics (1941) 137–150.
  * Krishnaiah and Rao [1961] P. R. Krishnaiah, M. M. Rao, Remarks on a multivariate gamma distribution, The American Mathematical Monthly 68(4) (1961) 342–346.
  * Krishnamoorthy and Parthasarathy [1951] A. S. Krishnamoorthy, M. Parthasarathy, A multivariate gamma-type distribution, Ann. Math. Statist. 22 (1951) 549–557.
  * Krupskii et al. [2018] P. Krupskii, R. Huser, M. G. Genton, Factor copula models for replicated spatial data, Journal of the American Statistical Association 113 (2018) 467–479.
  * Li and Sang [2018] F. Li, H. Sang, On approximating optimal weighted composite likelihood method for spatial models, Stat 7 (2018) e194.
  * Lim and Teo [2009] S. Lim, L. Teo, Gaussian fields and Gaussian sheets with Generalized Cauchy covariance structure, Stochastic Processes and Their Applications 119 (2009) 1325–1356.
  * Lindsay [1988] B. Lindsay, Composite likelihood methods, Contemporary Mathematics 80 (1988) 221–239.
  * Mai and Scherer [2014] J.-F. Mai, M. Scherer, Financial engineering with copulas explained, Springer, 2014.
  * Malov [2001] S. V. Malov, On Finite-Dimensional Archimedean, Birkhäuser Boston, Boston, MA, 2001.
  * Marshall and Olkin [1988] A. W. Marshall, I. Olkin, Families of multivariate distributions, Journal of the American Statistical Association 83 (1988) 834–841.
  * Masarotto and Varin [2012] G. Masarotto, C. Varin, Gaussian copula marginal regression, Electronic Journal of Statistics 6 (2012) 1517–1549.
  * McNeil and Nešlehová [2009] A. J. McNeil, J. Nešlehová, Multivariate Archimedean copulas, $d$-monotone functions and $\ell_{1}$-norm symmetric distributions, Annals of Statistics 37 (2009) 3059 – 3097.
  * Miller and Samko [2001] K. S. Miller, S. G. Samko, Completely monotonic functions, Integral Transforms and Special Functions 12 (2001) 389–402.
  * Morales-Navarrete et al. [2022] D. Morales-Navarrete, M. Bevilacqua, C. Caamaño-Carrillo, L. M. Castro, Modeling point referenced spatial count data: A poisson process approach, Journal of the American Statistical Association 0 (2022) 1–14.
  * Nelsen [2006] R. B. Nelsen, An introduction to copulas, Springer Science & Business Media, 2006\.
  * Pace et al. [2019] L. Pace, A. Salvan, N. Sartori, Efficient composite likelihood for a scalar parameter of interest, Stat 8 (2019) e222.
  * Palacios and Steel [2006] M. B. Palacios, M. F. J. Steel, Non-Gaussian Bayesian geostatistical modeling, Journal of the American Statistical Association 101 (2006) 604–618.
  * Pettorelli et al. [2005] N. Pettorelli, J. O. Vik, A. Mysterud, J.-M. Gaillard, C. J. Tucker, N. C. Stenseth, Using the satellite-derived NDVI to assess ecological responses to environmental change, Trends in ecology & evolution 20 (2005) 503–510.
  * Plemmons [1977] R. Plemmons, M-matrix characterizations. i – nonsingular m-matrices, Linear Algebra and its Applications 18 (1977) 175 –188.
  * Porcu et al. [2016] E. Porcu, M. Bevilacqua, M. G. Genton, Spatio-temporal covariance and cross-covariance functions of the great circle distance on a sphere, Journal of the American Statistical Association 111 (2016) 888–898.
  * Quessy and Durocher [2019] J.-F. Quessy, M. Durocher, The class of copulas arising from squared distributions: Properties and inference, Econometrics and Statistics 12 (2019) 148–166.
  * Quessy et al. [2016] J.-F. Quessy, L.-P. Rivest, M.-H. Toupin, On the family of multivariate chi-square copulas, Journal of Multivariate Analysis 152 (2016) 40–60.
  * Royen [2004] T. Royen, Multivariate Gamma distributions II, in: Encyclopedia of Statistical Sciences, New York: John Wiley & Sons, 2004, pp. 419–425.
  * Sibuya [1960] M. Sibuya, Bivariate extreme statistics, Annals of the Institute of Statistical Mathematics 11 (1960) 195 – 210.
  * Srivastava and Karlsson [1985] H. M. Srivastava, P. W. Karlsson, Multiple Gaussian Hypergeometric Series, Ellis Horwood Ltd, 1985.
  * Stein [1999] M. Stein, Interpolation of Spatial Data. Some Theory of Kriging, Springer-Verlag, New York, 1999.
  * Sudakov [2008] V. Sudakov, Lipschitz continuity of quantile functions on spaces of random variables, Journal of Mathematical Sciences 152 (2008) 941 – 943.
  * Suroso and Bárdossy [2018] S. Suroso, A. Bárdossy, Investigation of asymmetric spatial dependence of precipitation using empirical bivariate copulas, Journal of Hydrology 565 (2018) 685–697.
  * Varin [2008] C. Varin, On composite marginal likelihoods, Advances in Statistical Analysis 92 (2008) 1–28.
  * Varin and Vidoni [2005] C. Varin, P. Vidoni, A note on composite likelihood inference and model selection, Biometrika 52 (2005) 519–528.
  * Vere-Jones [1997] D. Vere-Jones, Alpha-permanents and their applications to multivariate gamma, negative binomial and ordinary binomial distributions, New Zealand J. Math 26 (1997) 125–149.
  * Williamson [1956] R. E. Williamson, Multiply monotone functions and their laplace transforms, Duke Math. J. 23 (1956) 189 – 207.
  * Xua and Genton [2017] G. Xua, M. G. Genton, Tukey g-and-h random fields, Journal of the American Statistical Association 112 (2017) 1236 –1249.
  * Zhang and El-Shaarawi [2010] H. Zhang, A. El-Shaarawi, On spatial skew-Gaussian processes and applications, Environmetrics 21(1) (2010) 33–47.