covatest : An R Package for Selecting a Class of Space-Time Covariance Functions

Although a very rich list of classes of space-time covariance functions exists, speciﬁc tools for selecting the appropriate class for a given data set are needed. Thus, the main topic of this paper is to present the new R package, covatest , which can be used for testing some characteristics of a covariance function, such as symmetry, separability and type of non-separability, as well as for testing the adequacy of some classes of space-time covariance models. These last aspects can be relevant for choosing a suitable class of covariance models. The proposed results have been applied to an environmental case study.

Until now, a space-time class of covariance models is selected in a subjective way, without considering the main characteristics exhibited by the empirical covariance surface. In some contributions of the last decade the authors analyzed some well-known properties of a covariance function, which can be suitably used in the choice of an appropriate class of models for a given variable under study (De Iaco, Posa, and Myers 2013;Cappello, De Iaco, and Posa 2018;Sherman 2007, 2008a).
Due to the increasing demand of spatial and spatio-temporal data analysis, there exist various packages which perform univariate and multivariate geostatistical analysis in two and three dimensions, such as Geo-EAS (Englund 1991), GEOSAN (Dutter 1996), the specific routine for assessing the fit of a variogram model (Barry 1996), the Geostatistical Software Library GSLib (Deutsch and Journel 1998) and its version for Windows WinGSLIB (LLC Statios 2001), SGeMS (Remy, Boucher, and Wu 2009), as well as the commercial software ISATIS (Geovariances 2015) and the cgeostat Fortran software for geostatistical analysis of complex-valued random fields (De Iaco 2017). Moreover, the R programming language (R Core Team 2020) has various packages devoted to Geostatistics, such as the gstat package for geostatistical modeling, prediction and simulation (Pebesma and Wesseling 1998;Pebesma 2004;Pebesma and Bivand 2005;Gräler, Pebesma, and Heuvelink 2016), the packages sp (Pebesma and Bivand 2005) and spacetime (Pebesma 2012;Bivand, Pebesma, and Gomez-Rubio 2013) for dealing with spatial and spatio-temporal data, geoR (Ribeiro and Diggle 2001;Diggle and Ribeiro 2007), which contains functions for model-based geostatistics, as well as the geospt (Melo, Santacruz, and Melo 2015), for geostatistical analysis and design of optimal spatial sampling networks. In addition, the R package RandomFields (Schlather et al. 2015(Schlather et al. , 2019 supports the modeling of multivariate random fields, as well as provides kriging, conditional simulation, covariance functions and maximum likelihood function fitting for a wide range of spatio-temporal covariance models. Padoan and Bevilacqua (2015) built the R package CompRandFld for the analysis of spatial and spatio-temporal Gaussian and binary data, and spatial extremes through composite likelihood inferential methods. The RGeostats (MINES ParisTech/ARMINES 2019) performs geostatistical operations simultaneously on several variables in a space of any dimension. The Geostatistical Team of Centre de Géosciences of Mines ParisTech also developed R2I which permits to transfer data between RGeostats and the software ISATIS. Ad-hoc routines for graphical interface in the ecological modeling software Bio7 (Austenfeld and Beyschlag 2012) or GeoXp for exploratory spatial data analysis (Laurent, Ruiz-Gazen, and Thomas-Agnan 2012) are also available. Moreover, a review of geostatistical softwares is given in Goovaerts (2010) and a list of some R packages on spatial statistics has been provided in Pebesma, Bivand, and Ribeiro (2015). A more complete list of packages on the Comprehensive R Archive Network (CRAN) related to spatial and spatio-temporal topics is available on CRAN Task Views "Spatial" (Bivand 2020) and "SpatioTemporal" (Pebesma 2020). Other contributions concern specialized routines and packages which perform spatio-temporal geostatistical analysis (De Cesare, Myers, and Posa 2002;De Iaco, Myers, Palma, and Posa 2010;De Iaco and Posa 2012;Pebesma 2012;Gabriel, Rowlingson, and Diggle 2013).
However, there is no routine or package that can help to cope with the problem of selecting of a suitable class of space-time covariance functions for a given data set. At most, existing packages provide criteria on how well some admissible models fit the empirical covariance function (e.g., mean difference between empirical and theoretical covariance surfaces, crossvalidation), but they give no support in picking out the class of covariance functions. This gap has been filled through the R package covatest, which provides tools for testing some features of space-time covariance functions (such as symmetry, separability and type of nonseparability) and some specific well-known classes of space-time covariance models. In other terms, the R package covatest provides some statistical measures on key properties of spatiotemporal covariance functions (De Iaco, Posa, Cappello, and Maggio 2019b). Moreover, some calculations provided by this package depend on other existing libraries, such as spacetime, gstat or stats, as well as lubridate (Grolemund and Wickham 2011) and zoo (Zeileis and Grothendieck 2005) for treating dates.
It will be shown that the above mentioned features support practitioners to properly choose a suitable class of space-time covariance functions for a given data set. In particular, the functions of the proposed package are based on the procedure given by Li et al. (2007) and Cappello et al. (2018), which are used for testing symmetry and separability, as well as for testing the type of non-separability (positive and negative) and some well known classes of spatio-temporal covariance models. This paper is organized as follows: in Section 2, a review on the testing procedure of some characteristics of a covariance function for stationary space-time random fields (STRF) is given. In Section 3 the statistical tests for some widely used classes of covariance models (such as the product-sum, the integrated-product and the Gneiting classes) are introduced. In Section 4 all the functions of covatest are briefly explained and the testing procedure is discussed. A more detailed overview of the functionality is given in the reference manual of covatest (De Iaco, Cappello, and Posa 2020), which is available from CRAN at https:// CRAN.R-project.org/package=covatest. Finally, in Section 5, the use of the implemented functions is shown through a case study regarding a well-known environmental database. It is worth pointing out that the package presented in this paper has a practical use in choosing a class of space-time covariance functions (not a specific model) for a given data set as well as the proposed examples can be easily extended to a wide range of applications.

Testing some features of space-time covariance functions
Before introducing the covatest package, it is convenient recalling the theoretical framework of the statistical tests used to check for the symmetry and separability assumptions, as well as to check for positive and negative non-separability and for some classes of space-time covariance models. Some details and applications of the testing procedure, even in the multivariate context, and a thorough analysis of size and power of these tests can be found in De Iaco, Palma, and Posa (2019a) and Cappello et al. (2018).
Let Z(s, t) be a random variable at the location s, in space, and the point t, in time, and let {Z(s, t) : s ∈ R d , t ∈ R} be a strictly stationary spatio-temporal random field (in the sense that the multivariate distributions are translation invariant), with zero expected value and finite covariance function which depends only on the spatial and temporal lags, denoted by h and u, respectively. The assumption of zero mean expected value reduces notation complexity with no loss of generality (Adler 1981), indeed by assuming that the mean function can be adequately modeled, one can reasonably focus on the covariance structure of the residual random field Z (Li, Genton, and Sherman 2008b).

Testing symmetry
It is well known that a stationary space-time covariance function C is fully symmetric if As previously specified, Li et al. (2007) provided a technique for testing the symmetry assumption, based on the asymptotic joint normality of the space-time covariance estimator of Z over D n = S × T n , with a fixed space S ⊂ R d and regularly spaced times T n = {1, ..., n}. Thus, given the set Λ of space-time lags, with cardinality equal to m, the following result is considered ∈ Λ} is the vector of space-time covariances at the lags in Λ.
• |T n | is the cardinality of T n .
The null hypothesis applied to check for full symmetry of a spatio-temporal covariance function is formulated as follows which could be written in compact form as H 0 : AG = 0 where A is a contrast matrix of dimension (q × m), whose row rank is q.
An example regarding the elements involved in the symmetry testing procedure is illustrated below.
Example 2.1 Given the set of space-time lags Λ = {(h 1 , u 1 ), (h 1 , −u 1 ), (h 2 , u 2 ), (h 2 , −u 2 )}, with m = 4, then the vector G of the covariances involved in the test and the contrast matrix A are defined as follows Under the hypothesis of symmetry, AG = 0.
Thus the test statistic used to check for the symmetry assumption is: which converges in distribution, under H 0 , to a χ 2 q with q degrees of freedom. At a significance level α, if the value of the test statistic is greater than the critical value, then the null hypothesis H 0 is rejected, meaning that there is a statistical evidence against the assumption of fully symmetry.

Testing separability
As widely discussed in the literature, a stationary space-time covariance function C is separable (Sherman 2011 or a spatio-temporal correlation function ρ is separable if In the contribution of Li et al. (2007) a nonparametric test for assessing the separability has been also proposed. In particular, if f = (f 1 , . . . , f r ) is a vector of real-valued functions that are differentiable at G, then using the asymptotic joint normality of G n and the multivariate delta theorem (Mardia, Kent, and Bibby 1979), the following asymptotic distribution is derived where the generic element of the matrix The hypothesis applied to check for the separability of a spatio-temporal covariance function is formulated as follows which could be written in compact form as H 0 : Af (G) = 0, where A is a contrast matrix of dimension (q × r), whose row rank is q and f (G) is a vector of r ratios (as the ones compared in the null hypothesis given in (8)).
An example regarding the elements involved in the separability testing procedure is illustrated below.
Example 2.2 Given the set of space-time lags Λ = {(0, 0), (h 1 , u 1 ), (h 1 , −u 1 ), (h 1 , 0), (0, u 1 )}, with m = 5, then the vector G of covariances involved in the test, the vector f (G) of the suitable ratios between the elements of G and the contrast matrix A are defined as follows Under the hypothesis of separability, Af (G) = 0.
The test statistic for checking separability is which converges in distribution, under H 0 , to a χ 2 q with q degrees of freedom. As explained in Li et al. (2007), the degrees of freedom depend on the number of space-time lags considered and thus on the number of contrasts. Size and power analysis of this separability test was provided by the same authors.
At a significance level α, if the value of the test statistic (9) is greater than the critical value, then the null hypothesis H 0 is rejected, meaning that there is a statistical evidence against the assumption of separability.

Testing for positive and negative non-separability
If the separability assumption is rejected, a non-separable covariance model is required and the investigation on the type of non-separability (i.e., computation of the non-separability index and the related testing procedure) is necessary. In De Iaco and Posa (2013) the authors proposed a generalization of the definitions of positive and negative non-separability given in Rodrigues and Diggle (2010).
Let ρ be the spatio-temporal correlation function, corresponding to the stationary covariance C. The non-separability index ratio is defined as follows where ρ(h, u) > 0, ρ(h, 0) > 0 and ρ(0, u) > 0. The spatio-temporal covariance function is uniformly positive or negative non-separable, if r(h, u) > 1 or r(h, u) < 1, respectively, for any spatio-temporal lag (h, u). Furthermore, if r(h, u) > 1, for some lags (h, u), the covariance function is pointwise positive non-separable at (h, u); alternatively, it is pointwise negative non-separable at (h, u), if r(h, u) < 1 for some lags (h, u). Then, a covariance function which is uniformly positive (negative) non-separable, is also pointwise positive (negative) nonseparable, but the converse is not true. Note that the type of non-separability can be also studied for a covariance function model, in this case the type of non-separability may depend not only on the lag (h, u) but also on the values of the parameters that characterize the model. Given the set of space-time lags Λ = {(0, 0), (h i , u i ), (h i , 0), (0, u i ), i = 1, . . . , l}, with cardinality equals to m = 3l + 1, with l ∈ N, l ≥ 1, the test statistic for the type of nonseparability is based on the following function where: • G = {C(h, u) : (h, u) ∈ Λ} is the vector of space-time covariances, whose cardinality is m.
• |T n | is the cardinality of T n .
• Σ = lim n→∞ |T n |cov( G n , G n ) is the asymptotic variance of the vector G n .
The null hypotheses on negative or positive non-separability and the corresponding alternative hypotheses are defined respectively as follows Note that a generic element of Af (G) is of the following type: At the significance level α, if H (-) 0 is set against H (+) 1 , as in (12), the test is conducted on the right tail (z α , +∞) of the standard normal distribution, or alternatively if H (+) 0 is set against H (-) 1 , as in (13), the test is conducted on the left tail (−∞, −z α ) of the standard normal distribution.
Failing to reject H (-) 0 formulated against H (+) 1 , given in (12), means that there is no statistical evidence that the space-time covariance function is positive non-separable for the specified lags; on the other hand, failing to reject H (+) 0 formulated against H (-) 1 , given in (13), means that there is no statistical evidence that the space-time covariance function is negative nonseparable for the specified lags. Note that for l = 1 the pointwise positive or negative nonseparability can be tested.
A simple way to fix the direction of the one-tailed test on the type of non-separability and the set of lags to be tested is to analyze as a first step the box-plots of the sample non-separability ratios, classified by spatial and temporal lags . Indeed, it is advisable to use a right tailed test for a given space-time lag, if the box-plots of the sample non-separability ratios support the assumption of negative non-separability (sample ratios less than one in average); analogously, for the left tailed test. According to the contribution of Cappello et al. (2018); Li et al. (2007), even the test on the type of non-separability can regard specific lags; thus, the corresponding decision of rejection/non rejection is referred to the conjectured type of non-separability for those lags (characterized by non-separability).
Moreover, if l > 1 only homogeneous lags, i.e., the spatio-temporal lags characterized by the same type of non-separability, should be used to perform the test.
An example regarding the elements involved in the test statistic for the type of non-separability is given below.

Example 2.3
Given the set of space-time lags Λ = {(0, 0), (h 1 , u 1 ), (h 1 , 0), (0, u 1 )}, with m = 4, then the vector G of covariances involved in the test, the vector f (G) of the suitable ratios between the elements of G and the contrast matrix A are defined as follows

Symmetry Separability Type non-separability Covariance models Positive
Negative Table 1: Properties of well-known spatio-temporal classes of covariance models. Remarks: • The asymptotic variance Σ is provided by applying a blocking technique (Section 2.3 in Li et al. (2007)).
• For the sake of simplicity it is assumed in the following that the expectation of Z is known and equal to zero; however if this assumption is removed, it is enough to denote with C * n (h, u) and G * n the mean-corrected estimators of C(h, u) and G, since G * n and G n have the same asymptotic properties (see Lemma A.6 in Li et al. 2008a andSherman 2011, p. 135).
In Table 1

Testing some classes of space-time covariances
As pointed out in the introduction, wide families of space-time covariance functions have been defined in the last years. They can be classified according to the type of non-separability (De Iaco and Posa 2013), in: • Uniformly non-separable models: for example, the Gneiting class of space-time covariance models (Gneiting 2002) as well as the class of space-time stationary covariance functions generated by positive power mixtures Ma (2002, Equation 2.3) are uniformly positive non-separable. Otherwise the product-sum class (De Iaco, Myers, and Posa 2001) is uniformly negative non-separable.
• Models with different types of non-separability: wide classes of non-separable stationary covariance functions, which can be uniformly positive or negative non-separable for some values of the parameters, or pointwise positive or negative non-separable, for some choices of spatio-temporal lags and parameters, belong to this family. For example, the classes generated by a linear combination of product-sum covariance models, the integrated product-sum covariance models, an example given by Cressie (1993), the Rodrigues and Diggle (2010) models, the more general positive power mixture model, proposed by Ma (2002, Equation 2.1), and the metric model (Dimitrakopoulos and Luo 1994), could be either uniformly positive and negative non-separable, or pointwise non-separable.
In Cappello et al. (2018) a non parametric testing procedure has been proposed to support the choice of a specific class of spatio-temporal covariance models. The tests on the classes of covariance models are based on some mathematical characteristics of the same classes which have been used to appropriately define the fundamental elements of the hypothesis testing (null hypothesis, contrast matrix and test statistics).
For this testing procedure the null hypothesis on a specified class of models can be formulated as H 0 : Af (G) = 0 and the test statistic is the one given in (9), which converges in distribution, under H 0 , to a χ 2 with degrees of freedom equal to the row rank of the contrast matrix A.
The testing procedure has been defined for some well-known and widely used classes of covariance models (the class of product-sum models, the class of integrated product models, the class of Gneiting models). For more details, see Cappello et al. (2018).

Testing the class of product-sum models
The class of product-sum models is defined as: where C s is a spatial covariance in As proved in the paper of De Iaco and Posa (2013), the class of product-sum models can only describe stationary space-time covariance functions that are uniformly negative non-separable.
Moreover, for any set of spatial lags h i , i = 1, 2, 3, and temporal lag u, the increments if and only if the class of models in (14) is considered.
Analogously, for any spatial lag h and any set of temporal lags u i , i = 1, 2, 3, Then, for any set of spatial lags h 1 , h 2 , and h 3 , the specific null hypothesis for the class of product-sum models is Analogously, for any set of temporal lags u 1 , u 2 , and u 3 , the specific null hypothesis for the class of product-sum models is An example regarding the elements involved in testing procedure on the class of product-sum models is given below.
} be the set of space-time lags, with m = 12, and G be the vector of covariances at any lag of Λ, then the elements of f (G) for the product-sum model and the contrast matrix A are defined as follows Under the null hypotheses (17) and (18), Af (G) = 0.

Testing the class of integrated product models
The following class of integrated product covariance models (De Iaco, Myers, and Posa 2002) with σ 2 > 0, smoothness parameters α, γ ∈]0, 1] and range parameters a, b > 0, can just describe stationary space-time covariance functions that are uniformly positive non-separable (De Iaco and Posa 2013).
For any set of spatial lags h i , i = 1, 2, 3, and temporal lag u, the increments 1 if and only if the class of models in (19) is considered.
Analogously, for any spatial lag h and any set of temporal lags u i , i = 1, 2, 3, Then, for any set of spatial lags h 1 , h 2 , and h 3 , the specific null hypothesis for the class of integrated product models is Analogously, for the set of temporal lags u 1 , u 2 , and u 3 An example regarding the elements involved in the testing procedure on the class of integrated models is given below.

Example 3.2 Let
} be the set of space-time lags, with m = 6, and G be the vector of covariances at any lag of Λ, then the elements of f (G) for the integrated product model and the contrast matrix A are defined as follows Under the null hypotheses (22) and (23), Af (G) = 0.

Testing a class of the Gneiting models
Consider the following family among the more general class of the Gneiting models (Gneiting 2002): with σ 2 , a, b ∈ R + ; γ, α, β ∈]0, 1] and d is the spatial dimension. As proved in De Iaco and , this class of models is characterized by uniform positive non-separability.
If the spatial and temporal marginals present a linear behavior near the origin (α = 0.5, γ = 0.5), for any set of spatial lags h 1 , h 2 , and h 3 in a bidimensional spatial domain (d = 2), such that ||h 1 || − ||h 2 || = ||h 2 || − ||h 3 ||, the specific null hypothesis for the Gneiting model is: or alternatively Analogously, for any set of temporal lags u 1 , u 2 , and u 3 , such that An example regarding the elements involved in the testing procedure on the class of the Gneiting models in (24) } be the set of spacetime lags, with m = 6, and G be the vector of covariances at any lag of Λ, then the elements of f (G) for the Gneiting model in (24) and the contrast matrix A are defined as follows Under the null hypotheses (25) and (27), Af (G) = 0.
Note that, since f (G) depends on the parameter β, the test can be computed by considering different admissible values of this parameter, whose domain is in the interval [0, 1]. However, if the hypothesis of separability is rejected, the value β = 0 (corresponding to the separable model) has to be excluded from the parameter space since it would not be consistent with the rejection of the separability hypothesis.
In the R programming language some packages are devoted to the spatio-temporal analysis.
In Table 2 a list of the most important R packages, their functionality and some implemented covariance models are provided.

Testing procedures
The main purpose of this paper is to provide some useful tools to support the R users in the choice of an appropriate class of space-time covariance functions for a given data set.
In this section the functions of the covatest package, which are relevant in the testing procedure, are described and the related outputs are analyzed.  arbitrary names for the outputs are proposed (these names can be changed by the users). It is evident that the outputs of some functions are arguments of other functions of the package.
The input data file of the implemented tests contains spatio-temporal data stored in: • A data.frame in which the data are in form of long formats (each record reflects a single time and space combination). This data base has to contain at least the (x, y) coordinates, the temporal codes/dates and the values of the variable under study.
• An object of class STFDF/STSDF (for more details see the package spacetime).
Alternatively, if the data for each spatio-temporal point are given by row in a text file, the utility function read.STdata(), available in the package, should be used in order to convert it into a data.frame or STFDF.
For each test, the R users have to define preliminary the couples of spatial points and the corresponding temporal lags to be analyzed, through the function couples().
It is worth pointing out that regarding the specific spatial points to be selected, the choice can be justified on the basis of different reasons. In general, the analyst can inspect the geometry of the sample points and start with few pairs of spatial points spread out over the domain, which are representative of the space-time correlation of the data set.
In some empirical cases, the selection of the pairs of spatial points might be based on intrinsic characteristics of the phenomenon under study with a preferential choice for those points with few missing values. For example, they can be chosen by taking into account the pairs of points with the smallest or, alternatively, with the largest ratio between the east-west component and the north-south component of the spatial lag, as well as the pairs of points with the shortest distance h . In particular, for wind speed data, the pairs of spatial points might be selected along the prevalent wind direction over the study area, as suggested for the test of separability by Li et al. (2007).
It is also common to combine various spatial lags or couples of spatial points with different temporal lags (usually from 2 to 6 lags) where the correlation is still strong, thus the choice does not fall on lags where the empirical covariance function decays. Moreover, for the test on the type of non-separability, it is advisable to avoid the spatial/temporal lags for which the sample non-separability ratios are close to one (therefore, before performing the test on the type of non-separability, it is advisable to execute the function sepindex() in order to compute the sample non-separability ratios).
Note that the structure of the spatio-temporal lags, which can be fixed through the function couples(), depends on the type of test; in particular: • For the test on symmetry, separability and type of non-separability, both positive and negative temporal lags are automatically considered for each spatial couple of spatial points. However, if the symmetry hypothesis is not rejected, only positive temporal lags might be considered for the test on separability and the type of non-separability, hence the setzero method must be used to neglect the negative temporal lags and thus to set them equal to zero.
• For the test on classes of space-time covariance functions, at least 3, or multiple of 3, of spatial and positive temporal lags are considered. For the test on the Gneiting class of space-time covariance functions the spatial points must be used to create at least 3, or multiple of 3, spatial couples such that each triplet of spatial lags h 1 , h 2 , h 3 satisfies the condition h 1 − h 2 = h 2 − h 3 , moreover 3, or multiple of 3, temporal lags must be fixed, such that each temporal triplet of lags u 1 , u 2 , u 3 satisfies the condition The spatial points and the corresponding temporal lags to be considered in the testing procedure can be fixed as follows couples(sel.staz = sel.staz, sp.couples.in = sp.couples.in, t.couples.in = t.couples.in, typetest = typetest, typecode = numeric()) where the argument • sel.staz is a vector of ID codes corresponding to the spatial points to be analyzed.
• sp.couples.in is a two-column matrix, each row corresponds to the couples of different spatial points, chosen among the ones fixed in sel.staz argument.
• t.couples.in is a vector of only positive (negative) temporal lags, different in absolute value.
• typetest denotes the test to be performed ("sym" for the symmetry test, "sep" for the separability test, "tnSep" for the type of non-separability test, "productSum", "intProduct" and "gneiting" for the model tests).
• typecode specifies the codification (numeric or character) of the spatial points in the data set under study.
Moreover, if some spatio-temporal lags are not required, they could be set equal to zero through the specific setzero method (see the manual of covatest for more details). The function couples() creates, as output, an object of the new class called couples. This object consists of five slots: • couples.st, matrix containing in the first two columns the couples of spatial points (denoted with order numbers, associated with the sequence of ID codes given in input) and in the other columns the temporal lags associated with each spatial couple.
• sel.staz, vector of the ID codes of the selected spatial points.
• sp.couples, matrix containing the couples of order numbers associated with the spatial points (the order is given by the way in which the spatial points are arranged in sel.staz) and the corresponding couples of the ID codes to be analyzed.
• tl.couples, vector of the temporal lags considered.
• typetest, the code of the test to be performed.
As discussed in Sections 2 and 3, the estimation of the matrix Σ, which is usually unknown, is required for each test. The sub-sampling technique, which consists of defining blocks by using a moving window along time (Li et al. 2007), has been implemented in the package for the estimation of the covariance matrix.
The matrix Σ can be estimated through the functions blocks() and covablocks().
The function blocks() allows the users to extract all possible overlapped blocks of data of the same length from the time series associated with the spatial points previously selected (see the slot stpairs@sel.staz): blocks(lb = lb, ls = ls, lt = lt, matdata = matdata, pardata1 = pardata1, pardata2 = pardata2, stpairs = stpairs) where the arguments are the length of each block (lb), the number of overlapped data between two consecutive blocks (ls), the spatio-temporal STFDF/STSDF or data.frame (matdata), and the spatio-temporal lags to be analyzed (i.e., stpairs, the output of couples(), see Table 3). Moreover, if the spatio-temporal dataset is given as STFDF/STSDF • pardata1 denotes the number of variables.
• pardata2 represents the slot in which the values of the variable to be analyzed are stored.
Otherwise if the spatio-temporal dataset is given as data.frame • pardata1 denotes the column in which the spatial IDs are stored.
• pardata2 represents the column in which the values of the variable are stored.
The number of blocks which can be extracted from the time series of each spatial point is computed as: where l t is the length of the time series of the spatial point. Note that a stop occurs if in the time series of each selected spatial point more than 75% of consecutive data are missing, or if in the blocks more than 80% of consecutive data are missing. Moreover, this function returns a warning message if more than 15% of the data in the last extracted block, of each selected spatial point, are missing. This warning message is very useful for practitioners since the estimation of the covariance matrix based on a large number of missing values is not reliable. If the warning message occurs, it is recommended to run again the function blocks() with different values of lb and ls in order to overcome the problem related to missing values.
The function blocks() creates, as output, an object of the new class called blocks. This consists of three slots: • array.block, array of dimension (lb × n b × n sp ), where lb is the length of each block, n b is number of blocks that can be extracted from the time series related to each spatial point and n sp is the number of spatial points defined in stpairs@sel.staz (slot sel.staz of the output stpairs). Each table of this array contains the blocks extracted for each spatial location.
• mat.block, matrix of dimension (lb × t n b ), where t n b = n b × n sp is the overall number of blocks that can be extracted from all time series related to the selected spatial points.
• sel.staz, vector of the ID codes of the selected spatial points.
The function covablocks() allows the user to compute the spatio-temporal covariances (for the specified lags) for each block of data and to estimate the matrix Σ: covablocks(stblock = stblock, stpairs = stpairs, typetest = typetest) where the arguments are the extracted blocks (i.e., stblock, the output of blocks(), see Table 3), the spatio-temporal lags to be analyzed (i.e., stpairs, the output of couples(), see Table 3) and the type of test to perform (typetest).
The function covablocks() creates, as output, an object of the new class called covablocks. This object consists of four slots: • mat.cova, matrix of sample spatio-temporal covariances for each block, computed for the spatial and temporal lags given in stpairs. Depending on the test the empirical variance, the sample spatial and temporal marginal covariances for each block of data are also computed.
• mat.cova.h, matrix of sample spatial marginal covariances for the spatial lags (spatial pairs) specified in stpairs.
• mat.cova.u, matrix of sample temporal marginal covariances for the temporal lags specified in stpairs.
• mat.cova.cova, matrix of sample covariances between space-time covariances computed for the blocks.
Depending on the type of the test, some of the above outputs (such as mat.cova.h and mat.cova.u) are not required, thus they are not computed.
Given the spatial and temporal lags specified in stpairs, the function covastat() determines the sample spatio-temporal covariance and other statistics useful for computing symmetry, separability and type of non-separability tests, whereas the function covastatM() determines the sample spatio-temporal covariance and other statistics useful for testing some classes of space-time covariance models. These functions can be called as follows covastat(matdata = data, pardata1 = pardata1, pardata2 = pardata2, stpairs = stpairs, typetest = typetest) covastatM(matdata = data, pardata1 = pardata1, pardata2 = pardata2, stpairs = stpairs, typetest = typetest, beta.data = NULL) where the arguments are the STFDF/STSDF or the spatio-temporal data.frame (matdata), two parameters regarding the dataset (pardata1 and pardata2, already specified in the description of the function blocks()), the spatio-temporal lags to be analyzed (stpairs, i.e., the output of class couples, see Table 3), the test to be performed (typetest).
The functions covastat() and covastatM() create, as output, an object of the new class called covastat and covastatM, respectively. These objects consist of several slots: • G, matrix containing the spatio-temporal covariances for the specified lags. For all tests, except for the symmetry test (typetest = "sym"), the sample variance and the sample spatial and temporal marginal covariances are also computed.
• cova.h, matrix containing the sample spatial marginal covariances for the spatial lags (i.e., spatial pairs) specified in stpairs.
• cova.u, matrix containing the sample temporal marginal covariances for the temporal lags specified in stpairs.
• f.G, array containing the computation of specific functions of the elements of G.
• B, matrix containing the computation of the derivatives of each element of f.G with respect to each element of G.
• A, the contrast matrix.
• beta.data, vector containing the different values of the parameter β, used for typetest = "gneiting" (slot available only for objects of covastatM class).
Depending on the type of the test, some of the above outputs (such as cova.h, cova.u, f.G, B and beta.data) are not required, thus they are not computed. Note that if the typetest, fixed as input of covablocks, covastat or covastatM, is different from the one used to define stpairs, a warning message informs the user about this discrepancy; however this message can be ignored if, by decision of the user, the same stpairs, previously defined for a fixed type of test, is used for other types of test.
Finally, the function covaprop() allows the users to test some properties (symmetry, separability, type of non-separability) of spatio-temporal covariance functions and some classes of space-time covariance models: covaprop(cblock = cblock, cstat = cstat, nonseptype = NULL, sign.level = 0.05) Note that (a) the objects of classes covablocks and covastat/covastatM (i.e., cblocks and cstat respectively, see Table 3) are used as input arguments of covaprop; (b) the specific test to be performed, through the function covaprop(), is implicitly retrieved from the inputs.
The function covaprop() creates, as output, an object of the new class called covaprop. This object consists of four slots: • test.statistic, • p.value, • df, the degrees of freedom (available only if the test statistic converges in distribution to a χ 2 ), • typetest, the code of the test performed.
Moreover, on the user's console a message appears: this helps to decide to either reject the null hypothesis in favor of the alternative or not reject it, at a 5% level of significance.
It is worth noting that the argument nonseptype in covaprop() is required only when the type of non-separability has to be tested (which implies that the arguments cblock and cstat in covaprop() have been defined for the same kind of test). This argument can assume values equal to 0 (for testing the null hypothesis that the non-separability is negative) or equal to 1 (for testing the null hypothesis that the non-separability is positive). The computation of the empirical non-separability index ratio through the function sepindex() and the inspection of the box-plots of the sample non-separability ratios, obtained through the boxplot() method, will help practitioner in setting nonseptype argument.
• nt and ns denote the number of temporal and spatial lags, respectively, in vario_st.
• globalSill is the value of the variance (estimated through the function var() of the package stats or through the limiting value of the sample variogram).
Note that the function returns a message of the percentage of non admissible values for the non-separability index (i.e., values less than zero), moreover the function sepindex() creates, as output, an object of the new class called sepindex. This object consists of four slots: • sep.index.ratio, the empirical non-separability index ratio, • cov.tm, the purely temporal sample covariances at the specified lags, • cov.sp, the purely spatial sample covariances at the specified lags, • cov.st, the spatio-temporal sample covariances at the specified lags.
Moreover, the boxplot() method, available for object of class sepindex, is useful to create the box-plots of the sample non-separability ratios, classified for spatial lags and temporal lags and to detect the type of non-separability.
It is important to point out that all the functions in covatest package return an alert if one or more arguments are not consistent with the admissible values of the parameters.

A case study
In this section the package covatest useful for choosing a class of space-time covariance functions has been applied to the well known data set airBase. This database, provided by the European Environmental Agency (http://acm.eionet.europa.eu/databases/airbase), is available with the package spacetime of the R environment.
In particular, the P M 10 daily data, measured from the 1st of January 2005 to the 31st of December 2006, at 13 rural background stations ( Figure 2) located in central Germany area have been retained since they are characterized by a very small percentage of missing data (less than 3%).

R> round(length(which(is.na(rr_13@data[1])))/9490, digits = 3)
[1] 0.023 The proposed testing procedure starts from checking for symmetry and separability (test statistics given in (4) and (9)). Then, if the hypothesis of separability is rejected, the type of non-separability (test statistic given in (11)) has to be checked. Finally, according to the previous results, one can proceed with testing the hypothesis on the suitable class of space-time covariance models.

Symmetry test
As pointed out in Section 4, each test requires, as a first step, to set the spatial and temporal lags through the function couples().
For this case study, the couples of stations to be used for this test have been formed on the basis of the criteria of the minimum distance. Then, the couples of stations with the smallest Euclidean distances have been selected through the function spDists() of the sp package and some pairs of spatial points have been chosen to perform the test (Table 4).
In particular, for testing the symmetry, 12 spatial points, hence 6 pairs of stations ( Figure 2) and 2 positive and 2 negative temporal lags, equal to ±1 and ±2, have been considered. The use of the function couples() and the appropriate settings are illustrated below.
Note that t.couples.in argument contains only positive temporal lags since the corresponding negative temporal lags will be set by default.
The information contained in couples.sym can be summarized by using the summary method, as shown below.

Number of temporal lags = 4 Number of spatial points involved = 12 Number of spatial couples = 6
After the definition of the spatio-temporal lags, the blocks, created through the technique of moving window along time, have to be used to estimate the covariance matrix Σ.
In couples.sym (object of class couples) 6 spatial couples (nlag sp = 6) and 4 temporal lags (nlag t = 4) have been considered, then the covariance for 24 lags (nlag sp × nlag t = nlag) has to be computed in order to perform the symmetry test. In particular, on the basis of this last aspect and taking into account the formula in (28), 24 blocks have been obtained by fixing the length of each block (lb) equal to 40 and the overlapped data between two consecutive blocks (ls) equal to 10: R> block.sym <-blocks(lb = 40, ls = 10, matdata = rr_13, pardata1 = 1, + pardata2 = 1, stpairs = couples.sym) A text message informs the user on the number of blocks obtained. The overall text output of block.sym is provided in the Appendix A.
The object blocks.sym of class blocks consists of three slots: • @sel.staz, vector of length equal to 12 (n sp ) of the ID codes of the selected spatial points.
• @array.block, array of dimensions (lb × n b × n sp ), where lb = 40, n b = 24 and n sp = 12. Each table of the array is associated with each spatial point in couples.sym. In each table there are as many rows as the length of blocks (in this case 40) and there are as many columns as the number of blocks (in this case 24).
As a further step, the sample spatio-temporal covariances, for the specified spatial and temporal lags given in couples.sym, and the contrast matrix have been determined by calling the function covastat() as follows: R> covast.sym <-covastat(matdata = rr_13, pardata1 = 1, pardata2 = 1, + stpairs = couples.sym, typetest = "sym") covast.sym@G is a matrix of dimensions (nlag × 1), which contains the 24 sample spatiotemporal covariances for the specified lags: R> covast.sym@G Furthermore, the contrast matrix of dimensions (npair ×nlag), where npair = 12 is the number of testing pairs (pairs of covariances on which the test is performed), could be inspected via covast.sym@A. Note that 12 testing pairs have been defined in couples.sym, consequently a contrast matrix of row rank 12 has been obtained. Each testing pair is of the typeĈ(h, u) andĈ(h, −u), as can be deduced from the null hypothesis in (3) This function returns the values of the test statistic and the p value equal to 2.184 and 0.999, respectively. Therefore there is no evidence of rejecting the null hypothesis of full symmetry.
Hereinafter the procedure to perform the separability test will be discussed. A complete demo for this test can be run through demo("symmetry", package = "covatest").

Separability test
For the test on separability the same spatial couples and temporal lags used for the symmetry test have been considered. Note that for the use of the function couples() the same procedure described in Section 5.1 has been followed.
However, negative temporal lags can be neglected since the hypothesis of symmetry has not been rejected in the previous step. Hence, all negative temporal lags have been set equal to zero through the specific setzero method, as shown below.
R> couples.sep <-setzero(x = couples.sep, zero = TRUE, value = 0) In this case, the zero argument has been fixed equal to TRUE since all negative temporal lags have to be set equal to zero (see the manual of covatest for more details). By taking into account that 6 spatial couples (nlag sp ) and 2 temporal lags (nlag t ) have been considered in couples.sep, approximately 12 (6 × 2) blocks have to be generated in order to perform the separability test. In particular 13 blocks (n b ) have been used by fixing lb and ls equal to 80 and 27, respectively: R> block.sep <-blocks(lb = 80, ls = 27, matdata = rr_13, pardata1 = 1, + pardata2 = 1, stpairs = couples.sep) A text message informs the user on the number of blocks obtained. The overall text output of block.sep is provided in Appendix A.
• @mat.cova.h, matrix of sample spatial covariances of dimensions (n b × nlag sp ), where n b = 13 and nlag sp = 6; in particular the 6 columns represent the number of spatial lags defined in couples.sep previously.
• @mat.cova.u, matrix of sample temporal covariances of dimensions (n b × nlag t ), where n b = 13 and nlag t = 2; in particular the 2 columns represent the number of temporal lags defined in couples.sep previously.
On the other hand, covast.sep@f.G is a matrix of dimensions (nfg ×1), with nfg = 14, which contains the ratios between the space-time covariances and the corresponding spatial marginal covariances (first 12 elements) and the ratios between the temporal marginal covariances and the sample variance (last 2 elements): R> covast.sep@f.G covast.sep@B is a matrix of dimensions (nlag * × nfg), with nlag * = 21 and nfg = 14, which contains the computation of the derivatives of each element of f (G) with respect to each element of G. While in the slot covast.sep@A is stored the contrast matrix, of dimensions (npair × nfg), with npair = 12 and nfg = 14. Note that 12 testing pairs (as many as the number of spatio-temporal lags specified in couples.sep) have been generated, consequently a contrast matrix of row rank 12 has been defined. Each testing pair is of the , as can be deduced from the null hypothesis in (8).
Finally, the function covaprop() is applied for testing the null hypothesis on the separability of the space-time covariance:  This function returns the values of the test statistic and the p value equal to 229.478 and 2.55e−42, respectively. Since the hypothesis of separability is rejected, hence a non-separable class of models is required for the data set under study. A complete demo for this test can be run through demo("separability", package = "covatest").
Hereinafter the procedure to detect and test the type of non-separability will be discussed.

Type of non-separability test
The type of non-separability can be first empirically analyzed through the function sepindex(), which preliminary requires the computation of the global sill and the spatio-temporal variogram. These two objects can be easily obtained through the functions var (stats package) and variogramST (gstat package), respectively, by using the commands below: R> C00_13 <-var(rr_13[,,"PM10"]@data[[1]], na.rm = TRUE) R> vv_13 <-variogramST(PM10~1, rr_13, width=60, cutoff = 220, tlags=0:15) The computed variogram is also available in the package covatest; thus it can be loaded and plotted by running the following command lines: R> data("vv_13") R> library("lattice") R> plot(vv_13, wireframe = T, zlab = NULL, xlab = list("distance (km)", + rot = 30), ylab = list("time lag (days)", rot = -35), + scales = list(arrows = F, z = list(distance = 5))) The object vv_13, shown in Figure 3, is characterized by 16 temporal lags (time lags from 0 to 15) and 4 spatial lags (i.e., 0, 30, 90 and 150). These data are essential for using the function sepindex(), which requires to specify the arguments ns and nt, corresponding to the number of spatial and temporal lags fixed to estimate the spatio-temporal variogram. Thus the settings for the function sepindex() are shown below: R> nonsep.index <-sepindex(vario_st = vv_13, nt = 16, ns = 4, + globalSill = C00_13) Warning message: 12.5% of non admissible values for the non-separability index (<0) have been removed Note that non valid values of the index have been found and removed, as pointed out in the warning message. Then, the new boxplot method, available for objects of class sepindex, can be used to produce the graphs of the non-separability index ratio, classified by spatial and temporal lags.
R> boxplot(nonsep.index, ylab="Non-separability ratio") Looking at the box-plots in Figure 4 it is evident that the sample ratios between the empirical space-time covariances and the product of the corresponding sample spatial and temporal marginals, are always less than one. Then a uniform negative non-separability can be detected and a right tailed test (12) can be applied.
Even for the test on the type of non-separability (typetest = "tnSep") it is necessary to specify first the spatial and temporal lags by applying the function couples(). For this test the spatial pairs used are the same as the ones used for the symmetry and separability tests and the temporal lags selected are equal to 3, 4 and 5 (only positive lags). It is worth nothing that the temporal lags 1 and 2 have not been selected since the non-separability index is very close to one (separability condition). Then, the following functions have been called: R> t.couples.in.tns <-c(3, 4, 5) R> couples.tns <-couples(sel.staz = sel.staz.sym, sp.couples.in = + sp.couples.in.sym, t.couples.in = t.couples.in.tns, + typetest = "tnSep", typecode = character()) This is a preview of the couples of the spatial points and the temporal lags to be analyzed.
As for the separability test, all negative lags have been excluded. Thus the setzero method has been used as follows: R> couples.tns <-setzero(x = couples.tns, zero = TRUE, value = 0) Details on the structure of couples.tns can be obtained by using the summary method: R> summary(couples.tns)

Number of temporal lags = 3 Number of spatial points involved = 12 Number of spatial couples = 6
For this test, the number of blocks can be chosen independently from the number of spatiotemporal lags fixed in couples.tns. In particular, 19 blocks (n b ) have been retained by fixing the length of each block (lb) equal to 60 and the overlapped data between two consecutive blocks (ls) equal to 23: R> block.tns <-blocks(lb = 60, ls = 23, matdata = rr_13, pardata1 = 1, + pardata2 = 1, stpairs = couples.tns) A text message informs the user on the number of blocks obtained. The overall text output of block.sym is provided in Appendix A.
The sample spatio-temporal covariances for the specified spatial and temporal lags given in couples.tns and other statistics useful for computing the test have been determined by calling the function covastat(): R> covast.tns <-covastat(matdata = rr_13, pardata1 = 1, pardata2 = 1, + stpairs = couples.tns, typetest = "tnSep") As clarified in Section 5.2, covast.tns@G is a matrix of dimensions (nlag * × 1), with nlag * = 28; covast.tns@f.G is a matrix of dimensions (nfg × 1), with nfg = 21, which contains the ratios between the space-time covariances and the corresponding spatial marginal covariances (first 18 elements) and the ratios between the temporal marginal covariances and the sample variance (last 3 elements) are provided. Moreover, covast.tns@B is a matrix of dimensions (nlag * × nfg) with nlag * = 28 and nfg = 21; while in the slot covast.tns@A is stored the contrast matrix, of dimensions (npair × nfg), with npair = 18 and nfg = 21. Note that 18 testing pairs (as the ones defined for the separability test) have been generated in covast.tns, thus a contrast matrix of row rank 18 has been defined.
Finally, the partial results (the objects covabl.tns and covast.tns of classes covablocks and covastat, respectively) have been used as input arguments in covaprop() in order to test the null hypothesis that the empirical space-time covariance function is negative nonseparable (nonseptype = 0): This function returns the values of the test statistic and the p value equal to −0.626 and 0.734, respectively. Hence the null hypothesis of negative non-separability is not rejected, at 5% significance level. A complete demo for computing the non-separability index and performing this test can be obtained through demo("typenonseparability", package = "covatest").

Testing models
On the basis of the previous results, the class of product-sum models (i.e., uniformly negative non-separable covariance model) can be suitably selected in order to describe the spacetime correlation exhibited by the data. For the sake of exposition the test on the class of covariance models will be applied on the product-sum class, but also on the Gneiting and the integrated product classes. These last two classes are not consistent with respect to the type of non-separability empirically detected, thus it is expected that the hypothesis on the appropriateness of these classes will be rejected. To perform the tests on these classes of models (typetest = "productSum" or "intProduct" or "gneiting"), 9 pairs of stations, hence 3 spatial triplets (each triplet forms a group of 3 pairs of spatial points), and 3 temporal lags (equal to 1, 2 and 3), thus 1 temporal triplet, have been considered. Note that, although there are no constrains in setting the spatial and temporal triplets for the tests on the productsum model and integrated product model, the spatial and temporal triplets have been fixed according to the standards required by the test on Gneiting model (see Section 3 for details). In this way, the same object of class couples can be used to perform all the 3 tests on the selected classes of models.
The use of the function couples() and the appropriate settings for the test on the productsum class of models (typetest = "productSum") are illustrated below.
A contrast of type (18), (23) and (27) can be defined if there is a temporal triplet (3 non zero temporal lags), for each spatial couple (i.e., for a fixed h). Analogously, a contrast of type (17), (22) and (25) can be defined if there is a spatial triplet, for each temporal lag (i.e., for a fixed u). Thus, 6 contrasts can be obtained for each spatial triplet and each temporal triplet. In order to avoid redundant contrasts (see the manual of covatest for more details), 3 temporal lags (one for each spatial triplet) can be arbitrarily set equal to zero through the setzero method, hence 4 contrasts, instead of 6, can be conveniently defined: R> zero.index <-matrix(data=c(3, 7, 6, 7, 9, 7), ncol=2, byrow = TRUE) R> couples.mod <-setzero(x = couples.mod, zero = FALSE, index = zero.index, + value = 0) On the basis of the spatial couples and the temporal lags selected, 12 comparisons have been set and 14 blocks have been retained by fixing the length of each block (lb) equal to 60 and the overlapped data between two consecutive blocks (ls) equal to 10: R> block.mod <-blocks(lb = 60, ls = 10, matdata = rr_13, pardata1 = 1, + pardata2 = 1, stpairs = couples.mod) A text message informs the user on the number of blocks obtained. The overall text output of block.mod is provided in Appendix A.
The object block.mod has been used as input value for the estimation of the space-time covariance in covablocks().
R> covabl.ps <-covablocks(stblocks = block.mod, stpairs = couples.mod, + typetest = "productSum") Similarly to the separability and type of non-separability tests, the object covabl.ps (object of class couples) consists of four slots. The sample spatio-temporal covariances for the specified spatial and temporal lags given in couples.mod and other statistics useful for computing the test have been determined by calling the function covastatM: R> covast.ps <-covastatM(matdata = rr_13, pardata1 = 1, pardata2 = 1, + stpairs = couples.mod, typetest = "productSum", beta.data = NULL) R> test.ps <-covaprop(cblock = covabl.ps, cstat = covast.ps, + nonseptype = NULL, sign.level = 0.05) Finally, the partial results (the objects covabl.ps and covast.ps, of classes covablocks and covastatM, respectively) have been used as input values in covaprop() in order to test the adequacy of the product-sum class of models: This function returns the values of the test statistic and the p value equal to 7.214 and 0.843, respectively. Then, it is clear that the null hypothesis of the adequacy of the product-sum class is not rejected, at 5% significance level.
In this case, the user has to run again covablocks() and covastatM() and set the appropriate code of typetest.
The command lines for testing the integrated product class are given below: R> covabl.ip <-covablocks(stblocks = block.mod, stpairs = couples.mod, + typetest = "intProduct") Warning message: the argument typetest is different from the one defined in stpairs R> covast.ip <-covastatM(matdata = rr_13, pardata1 = 1, pardata2 = 1, + stpairs = couples.mod, typetest = "intProduct", beta.data = NULL) Warning message: the argument typetest is different from the one defined in stpairs Note that the above warning messages inform the user that the input argument stpairs used in covabl.ip and covast.ip has been previously defined for testing another model (i.e. the product-sum model). This is intentional, since at the beginning of this section it has been pointed out that the same pairs of stations and temporal lags are considered for testing the selected class of models. For the integrated product class (typetest = "intProduct"), the test statistic and the p value have been equal to 53.614 and 3.202e − 07, respectively. Then it is clear that the null hypothesis of the adequacy of the integrated product model is rejected, at 5% significance level.
Afterwards, the test on the Gneiting class of models (typetest = "gneiting") has been performed. Note that in this test, it is necessary to specify in covastatM() the value of the parameter β (beta.data argument).
R> covabl.gn <-covablocks(stblocks = block.mod, stpairs = couples.mod, + typetest = "gneiting") Warning message: the argument typetest is different from the one defined in stpairs R> covast.gn <-covastatM(matdata = rr_13, pardata1 = 1, pardata2 = 1, + stpairs = couples.mod, typetest = "gneiting", beta.data = 1) Warning message: the argument typetest is different from the one defined in stpairs Similarly to the previous case, the above warning messages inform the user that the input argument stpairs in covabl.gn and covast.gn has been defined for testing the product-sum model whereas in this part of the case study is used for testing another type of model, i.e. the Gneiting model. For the Gneiting class (typetest = "gneiting"), the test statistic and the p value have been equal to 414.175 and 3.760e − 81, respectively. It is clear that this class of models is not able to describe the spatio-temporal correlation exhibited by the data, since the null hypothesis is rejected. A complete demo for these tests can be run through demo("typemodel", package = "covatest").

Conclusions
The complexity of spatial-temporal processes does not often allow an adequate description of these phenomena through well-known physical laws; in these cases, it would be useful to have a guide for selecting a suitable class of space-time covariance models for a given data set. At most, existing packages provide criteria on how well some admissible models fit the empirical covariance function (e.g., mean difference between empirical and theoretical covariance surfaces, cross-validation), but they give no support in picking out the class of covariance functions. Thus, the main topic of this paper has been to propose the R package covatest in order to provide some statistical measures on key properties of spatio-temporal covariance functions and consequently to choose an appropriate class of space-time covariance functions, on the basis of some characteristics, such as symmetry, separability and type of non-separability. Moreover, tests on some well-known classes of space-time covariance models (easy extendable to some other classes) have also been provided. These last aspects have never been analyzed in the literature and support the modeling decision.

A. Supplementary output
In this section the whole output the functions couples() and blocks(), previously shown in Section 5, are provided.