Computational advances for spatio-temporal multivariate environmental models

In multivariate Geostatistics, the linear coregionalization model (LCM) has been widely used over the last decades, in order to describe the spatial dependence which characterizes two or more variables of interest. However, in spatio-temporal multiple modeling, the identification of the main elements of a space–time linear coregionalization model (ST-LCM), as well as of the latent structures underlying the analyzed phenomenon, represents a tough task. In this paper, some computational advances which support the selection of an ST-LCM are described, gathering all the necessary steps which allow the analyst to easily and properly detect the basic space–time components for the phenomenon under study. The implemented algorithm is applied on space–time air quality data measured in Scotland in 2017.


Introduction
In environmental studies, the analysis of the relationships among two or more variables, measured in space and time, is of wide interest. For prediction purposes, the first stage to be tackled in Geostatistics concerns the identification of a suitable multivariate model which describes the spatio-temporal covariance matrix computed for the observed variables. As thoroughly reviewed by Genton and Kleiber (2015), one can choose among various models of spatial multiple correlation; however, the LCM is still very common in practice, since it is computationally flexible and convenient (Babak and Deutsch 2009;Emery 2010;Gneiting et al. 2010;Bevilacqua et al. 2015).
Its adequacy can be also tested by applying the statistics proposed by Li et al. (2008).
B Sandra De Iaco sandra.deiaco@unisalento.it Nevertheless, in the space-time context the definition of the main elements of the ST-LCM, namely the basic structures characterizing the variables of interest, represents yet a tough task, due to the lack of software to analyze multivariate spatio-temporal data. A first attempt to propose a procedure to support the identification of an ST-LCM for a multivariate data set was given in De Iaco et al. (2012), where the authors assumed to model the basic components through the product-sum class. Although the restriction on this specific class of covariance models to be used for describing the basic components was abandoned in some further applications, only the contribution in De  pointed out that the basic components can be modelled by using different families of covariance models and the choice of the family for one basic component does not influence the choice for the other components, since the selection of each class can be reasonably relied on the main features of the sample basic covariance surfaces. However, the procedure of choosing suitable LCMs for spatio-temporal data was presented without discussing the related computational issues, as also clarified by the same authors. Thus, in this paper, the proposed algorithm focuses, in a very compact and analytic way, on the implementation of the multivariate modeling procedure of an ST-LCM, where no specific class of covariance models is assumed a priori for the basic components. Moreover, differently from the common approach which consists in using the same class of models for all the components of the ST-LCM as in De Iaco et al. (2003), the new algorithm allows the users to obtain a flexible ST-LCM, whose components can have different properties in terms of symmetry/asymmetry, separability/non-separability and, in case of non-separability, positive/negative nonseparability . To our knowledge, these advances have not been tackled in previous works and represent a remarkable step forward of the current state of computational progress in spatio-temporal multivariate analysis. The steps of the presented algorithm can be easily carried out by implementing some ad-hoc R functions (R Core Team 2015), as done by the authors.
The computational advances are applied to model the spatio-temporal multivariate covariance function of three different pollutants, which are still critical in Scotland, i.e. nitrogen dioxide (NO 2 ), nitric oxide (NO) and fine particles (PM 10 ). The data consist of hourly measurements recorded in 2017 at 18 monitoring stations around Glasgow. However, the same algorithm can be useful even for other environmental applications when the identification of a spatio-temporal multivariate model, able to capture the common latent structures underling the analyzed variables, is required.
The paper outline is as follows. Section 2 provides some basic notions on multivariate spatio-temporal Geostatistics; Sect. 3 offers a detailed description of the algorithm for the ST-LCM selection, then, the application on the data set concerning the above-mentioned air pollutants is presented through the various steps of the proposed algorithm and some interesting practical aspects are highlighted (Sect. 4). Then, the adequacy of the chosen ST-LCM for the studied data set is evaluated and successively the model is used for prediction purposes. Finally, in Sect. 5, the proposed procedure is compared with respect to the exiting fitting method based on the use of the same class of covariance models for all the basic components.

A synthetic multivariate framework
Let X be a multiple second-order stationary random function (MSRF), defined on a spatio-temporal domain D × T ⊆ R d × R, which consists of m components X 1 (s, t), . . . , X m (s, t), with (s, t) ∈ D × T and m ≥ 2. Its first and second order moments are respectively: where μ i is the expected value of corresponds to the direct and cross-covariance functions of the MSRF X. In particular, for i = j, C ii is the direct covariance function of the X i , while for i = j, C i j is the cross-covariance function between X i and X j . As already clarified in the introduction, the ST-LCM is still extensively used to model the matrix C in various scientific fields. This is done through the following linear combination of basic covariance functions: where c l (δ s , δ t ) are covariance functions and B l = [b l i j ], l = 1, . . . , L, with L ≤ m, are (m × m) positive definite matrices and are called coregionalization matrices. Note that c l can be interpreted as the covariance functions of the latent variables which characterize the random field X.

Algorithms for an ST-LCM selection
The main algorithm, shown below, recalls the following steps: I. detection of the basic components determined through simultaneous diagonalization of the covariance matrices; II. identification of suitable models for the basic components, by taking into account the results concerning hypothesis testing of their main features; III. computation of coregionalization matrices and check of their admissibility through the singular value decomposition. The above-mentioned steps can be executed in the R environment.

Algorithm 1: Identify basic components
The first algorithm has the objective to determine the basic components of the MSRF X and to identify the corresponding covariance functions and their scales of variability.
Given the spatio-temporal data set x of size (N × m) related to the N observations of m variables over the spatio-temporal domain, let C(δ s , δ t ) k = [ C i j (δ s , δ t ) k ] be Main algorithm 1: Load the multivariate data set of m ≥ 2 spatio-temporal variables 2: Recall the Algorithm 1: "Identify basic components" 3: For each basic component l (l ≤ L), test the empirical characteristics of the sample basic covariance surface and define the appropriate model by recalling the Algorithm 2: "Test main features and select basic models" 4: Recall the Algorithm 3: "Compute coregionalization matrices" 5: Define the ST-LCM by using equation (2) the symmetric (m × m) matrix of direct and cross-covariance functions estimated at each lag k, with k = 1, . . . , K . Then, the matrices C(δ s , δ t ) k are simultaneously diagonalized through a (m × m) orthonormal matrix , i.e.
where D(δ s , δ t ) k are diagonal, or nearly diagonal, (m × m) matrices. Among the different methods proposed in the literature to perform simultaneous diagonalization of several symmetric matrices, the one used in the present work is based on an extension of the Jacobi technique where the criterion of simultaneous diagonalization is iteratively minimized under plane rotations applied to the matrices to be diagonalized (Cardoso and Souloumiac 1996).
Given the elements d i j,k of the matrices D(δ s , δ t ) k , with i, j = 1, . . . , m, k = 1, . . . , K , for the K spatio-temporal lags, the goodness of the diagonalization procedure is assessed by evaluating the following index: Then, the matrix is used to compute the empirical independent components y l , characterizing the data matrix x as follows: so that the m basic covariance functions c l are estimated, at K user-chosen space-time lags, from the m latent components y l and denoted withĉ l , l = 1, . . . , m. Note that y l indicates the column vector, of size N , of the (N × m) matrix y.
A joint visual inspection of the plots of the covariance surfacesĉ l and their spatial and temporal marginals, is useful to detect the different scales of variability, which correspond to the different lags where the surfaces of the basic covariances decay. Thus, the analyst is able to define the L ≤ m distinct scales of spatio-temporal variability and to select the L ≤ m basic components of the ST-LCM.
The instructions of Algorithm 1 can be specified as follows.
Algorithm 1: Identify basic components 1: Load the multivariate data set of m ≥ 2 spatio-temporal variables and the coordinates of the spatial locations 2: Load the set of K spatio-temporal lags 3: Load the matrix (K × m 2 ) of the spatio-temporal direct and cross-covariance estimates 4: Perform the simultaneous diagonalization of the real-valued (m × m) matrices of the sample direct and cross-covariances estimated at the K user-defined spatio-temporal lags 5: Determine the diagonality index as in (4) 6: Obtain a realization of the m independent components y l from equation (5) 7: Determine the space-time covariance surfaces of the m independent components 8: Plot the space-time covariance surfaces and determine (even graphically) the corresponding scales of variability. Then compare the scales of variability and define the L ≤ m basic components with distinct scales of variability

Algorithm 2: Test main features and select basic models
The second algorithm supports the steps for assessing the main features of the basic components; it ends with the selection of the appropriate classes of covariance models as well as with their fitting. First of all the full symmetry hypothesis is tested. In particular, the null hypothesis is formulated as with (δ s , δ t ) ∈ Q l , (δ s , −δ t ) ∈ Q l , l = 1, . . . , L, where Q l denotes the set of the space-time lags fixed by the analyst at each scale l, where its cardinality is p l ( p l ≥ 2). Alternatively, the hypothesis in (6) is also expressed as H 0 : A l G l = 0, where A l is usually called contrast matrix (with rank(A l ) = r l ) and G l is the vector of the covariances c l at the lags (δ s , δ t ) ∈ Q l . By recalling the assumptions in Li et al. (2007), the hypothesis (6) can be tested by applying the following statistic where G ln is the vector of the covariance estimators based on a sample extracted from the domain D × T n , with D ⊆ R d and T n = {1, . . . , n}, |T n | is the cardinality of the set T n and l = lim n→∞ |T n |Cov( G ln , G ln ) and χ 2 denotes the chi square distribution.
If the full symmetry test does not end with the rejection of H 0 , then the separability condition can be tested.
Analogously to the previous test, it is necessary to specify the null hypothesis of separability, i.e.
with (δ s , δ t ) ∈ Q l , l = 1, . . . , L, and the test statistic is where r l = rank(A l ), B l is a matrix whose elements are Harville 2001). Subsequently, if separability is rejected, it is convenient to investigate on the type of non-separability. To begin with, the empirical non-separability ratios are computed, i.e.
Depending on the values that this index takes on, the analyst can draw some tips regarding the type of non-separability. Indeed, • if the index r l (δ s , δ t ) is greater (smaller) than 1 for all the given lags, the uniform positive (negative) non-separability is reasonably suggested; • otherwise non uniform non-separability can be supposed.
During this step, the box diagrams of the indexes in (10) can support the formalization of the hypotheses to be tested, i.e.
or alternatively where 1 = [1] r and A l is the associated contrast matrix (with rank(A l ) = r l ), l = 1, . . . , L. These hypotheses can be tested by using the test statistic T S 3l , l = 1, . . . , L, which is obtained through the standardization of the following statistic (Cappello et al. 2018) The instructions of Algorithm 2 are listed in the proper frame.
Algorithm 2: Test main features and select basic models 1: Load the l-th basic component 2: Test the validity of full symmetry of the l-th component 3: If the hypothesis of full symmetry is not rejected, the separability assumption must be verified 4: If the hypothesis of separability is unlikely from the empirical evidence and it is rejected, compute the sample non-separability ratios, plot the values classified by spatial and temporal lags, then test the type of non-separability 5: Choose the class of covariance models according to the tested features and fit the selected model 6: Compute indexes of fitting goodness

Algorithm 3: Compute coregionalization matrices
The third algorithm concerns the computation of the elements of the (m × m) coregionalization matrices B l , l = 1, . . . , L.
The entries of B l represent the contribution to be associated to each basic component at the distinct scales of variability. To this aim, it is important to fix the value reached by C i j (for each i, j = 1, . . . , m) at each selected scale of variability. This can be determined, for each i, j = 1, . . . , m, by averaging the sample covariances C i j that fall in a neighborhood of the l-th scale of variability.
Then, the elements b l i j of matrices where C i j (δ s , δ t ) 0 = C i j (0, 0). Note that the difference at the numerator corresponds to the net contribution of C i j at each scale of variability. However, the matrices B l are admissible if they are positive definite. Consequently, it is necessary to check their validity or equivalently, to verify if their eigenvalues are non negative. This can be carried out through the spectral decomposition of the matrices B l , i.e.
where V l is the eigenvectors matrix and l is the diagonal matrix of the eigenvalues. In case of negative eigenvalues, it is suggested to modify them by forcing the negative eigenvalues to be equal to zeros and to produce the following admissible matrix B l : where the new diagonal matrix of the eigenvalues l is obtained from l by substituting the negative eigenvalues with zeros. The instructions of Algorithm 3 are specified in the proper frame. Note that there is not a specific order in executing the last two algorithms (2 and 3), since the corresponding results do not depend from each other.

An application on air pollutants
The present case study aims to describe the various steps of the algorithm previously proposed, in order to identify an apt ST-LCM for the data of interest. The analyzed variables concern air pollution which still represents one of the major health risks, notwithstanding the introduction of national and supranational regulations for the environmental protection and control.
Moreover, as air pollution is the result of the interaction of various contaminants, which change in space and in time, a proper and punctual air quality assessment requires a deep multivariate analysis of the joint spatio-temporal behavior of the variables of interest.
In the proposed case study, the evolution in space and time of NO 2 , NO and PM 10 hourly concentrations measured in 2017 at 18 environmental monitoring stations in Scotland is analyzed in order to detect possible latent independent components which are common to the variables and describe their direct and cross-correlation in spacetime.

The sampling data: main features
In Scotland, the concentration levels of the most hazardous air pollutants, such as particulate matter (i.e. PM 10 , PM 2.5 ) and nitrogen oxides (i.e. NOx, NO 2 ) have gradually reduced over the past decades, nevertheless the limit values fixed by the European Union for these pollutants are often exceeded at some Scottish locations. In 2017, the Lancet Commission on pollution and health, declared "Glasgow as one of the most polluted parts of the UK".
The analyzed data, which are openly available at the following web site http://www. scottishairquality.scot/data/, come from the air quality monitoring sites operated by Defra, the Scottish Government and Local Authorities (Defra 2020). They refer to hourly measurements of NO 2 , NO and PM 10 , expressed in µg/m 3 . The data have been collected in 2017 in Scotland, at 18 monitoring stations close to Glasgow (Fig. 1), most of which are named "urban traffic stations" due to their location at city centers and close to roadways where the main emission sources are vehicles and home heating. On the other hand, the "urban background station" is located in an urban residential area and reflects mainly the city-wide background conditions; while the "rural background station" is situated in countryside where the principal emission source consists of the regional long-range transport.
As known, several studies have shown that long-term exposures of these air pollutants caused negative effects on human health, such as respiratory as well as acute and chronic cardiovascular outcomes (Ayyagari et al. 2007;Hoek et al. 2002). For these reasons the interest of the scientific community is mainly in studying and modeling the joint behaviour of these hazardous pollutants, which exhibit the same evolution in time, with a daily periodicity characterized by very low values early in the morning and high values during the work hours (from 6 a.m. till 6 p.m.), as shown in Fig. 2. Moreover, it is evident that the lowest mean values for the three pollutants occurred The 24-h periodicity, which immediately appears in Fig. 2 for the variables under study, has been removed, station by station, by using the moving average technique applied to homogeneous subsets of the temporal span. Note that the program REMOVEMULT (Brockwell and Davis 1987;De Iaco et al. 2010), utilized to treat the systematic component in time, performs also a linear time interpolation in the presence of a maximum of 5 consecutive missing values. However, it is worth pointing out that, in the analyzed data set, the number of missing values is very low, indeed for all the stations the percentage of missing values is on average at most 4.3%, with a maximum of 14.2% for one of the stations. Moreover, there are few sequences of 5 or more consecutive missing values (on average, 5 sequences for all the stations).
The effects of various factors (such as traffic volumes, seasonality, working/nonworking days) which may affect the air pollution magnitude and cycle have been essentially removed by filtering out, station by station, the systematic component. In addition, it is worth pointing out that, on the basis of a descriptive analysis of the data, the altitude of the locations has no significant effects on the daily behavior of the pollutants, this is mainly due to the orography of the area under study which is also not very extended. Thus the NO 2 , NO and PM 10 residuals have been examined in order to detect the most appropriate ST-LCM, by following the procedure previously described.

Basic components' identification
As described in Sect. 3, the first step towards the identification of the most suitable ST-LCM, concerns the computation of the sample covariances, both direct and cross, for a proper selection of lags fixed by the researcher on the basis of the spatio-temporal sampling.
For the present case study, 8 lags in space and 37 lags in time, in total 296 lags in space-time, have been chosen, thus the spatio-temporal direct and cross-covariance surfaces have been computed at the selected lags for NO 2 , NO and PM 10 residuals (Fig. 3).
Successively, the 296 symmetric (3 × 3) matrices of the estimated covariance functions have been nearly diagonalized through the package Jade (Miettinen et al. 2017) developed in the R environment (R Core Team 2015). The diagonalization accuracy has been quantified by the relative index (4) determined to measure the closeness of the nearly diagonalized matrices to diagonal matrices. In this case, index (4) has shown very low values (approximately 8 out of 10 values are less than 0.01) confirming thus the diagonalization's goodness.
Then, by the product between the matrix of the NO 2 , NO and PM 10 residuals and the following orthonormal matrix which has nearly diagonalized the 296 sample covariance matrices, three sample independent latent components have been computed, as indicated in (5). Afterwards, the spatio-temporal covariance functions have been estimated for the obtained independent components and, through the visual exploration of their surfaces, two basic components (L = 2) have been retained since they have exhibited different scales of variability in space-time: (a) a short distance variability scale around 10 km and 12 h, which represents the short distance component (SDC), and (b) a long distance variability scale around 13 km and 24 h which represents the long distance component (LDC). The remaining component has been discarded considering that its spatio-temporal variability is quite similar to one of the chosen components.
The covariance surfaces corresponding to the two basic components are shown in Fig. 4 together with the respective marginals in space and time, whose behavior is linear at the origin and concave for small spatial and temporal lags.

Basic components' modeling
As second step of the proposed selection procedure, the main features (full symmetry, non-separability and type of non-separability) of the chosen basic components, need to be tested. In particular, the aforementioned tests have been performed through the covatest package (Cappello et al. 2020). The results of these tests immediately indicate the most apt classes of covariance models for the retained components.

Main characteristics of the selected basic components: tests' results
The full symmetry test (7) has been carried out for each basic component on the set of lags where the correlation in space and time was significant; in particular spatial lags equal to 2.2 and 4.5 km, and temporal lags equal to 1 and −1 h, have been fixed.
At 5% significance level, it results that the null hypothesis of full symmetry cannot be rejected for both components, since T S 11 = 3.561 ( p value = 0.169), for the SDC and T S 12 = 1.113 ( p value = 0.573), for the LDC.
Then the separability condition has been checked after selecting the following spatio-temporal lags: On the basis of the above results it is reasonable to assume that the two selected components are fully symmetric and non-separable; hence the type of non-separability has been checked. For this purpose, the non-separability ratios (10) have been computed and classified per spatial and temporal lags (Fig. 5). This graphical tool represents a very simple way to have information on the type of non-separability which will be successively tested. In this case, the non-separability ratios (10) computed for the SDC are uniformly smaller than 1, thus a uniform negative non-separability can be reasonably assumed; while for the LDC, the sample ratios (10) are greater than 1 for all spatio-temporal lags, hence a positive non-separability can be assumed.
At the significance level α = 0.05, the test in (13) is conducted, for the hypotheses (11), on the right tail of the standard normal distribution for the SDC, while the test  (13) is applied, for the hypotheses in (12), on the left tail of the standard normal distribution for the LDC. In particular, after choosing different spatial lags (namely, lags equal to 2.2 km, 4.5 km, 6.5 km and 8.9 km) and temporal lags (namely, from 1 to 6 h), the null hypotheses previously fixed are not rejected at 5% significance level.

Modeling the selected basic components
The results concerning the test on the main characteristics of the two basic components drive the analyst towards the identification of the most appropriate covariance models . In this case study, there is a statistical evidence to model the surface of the SDC with a fully symmetric, negative non-separable space-time covariance function, while the LDC's surface can be described by a fully symmetric, positive non-separable space-time covariance function. In particular, the product-sum covariance function, which is a fully symmetric and negative non-separable spacetime covariance model (De Iaco et al. 2001;, can be reasonably used for SDC (l = 1) with k 11 > 0, k 21 ≥ 0, k 31 ≥ 0, C s 1 the spatial exponential covariance model in R d , (d = 2) and C t 1 the temporal exponential covariance model in R, with practical ranges a s 1 and a t 1 , respectively. On the other hand, the Gneiting class of model (Gneiting 2002), which is characterized by full symmetry and positive non-separability (De Iaco and Posa 2013), has where a s 2 and a t 2 are the practical ranges in space and time, respectively, while N and C are positive values representing the contributions to the variance and β ∈ [0, 1]. The parameters' estimates of the above mentioned covariance models are shown in Table 1: they have been obtained by carrying out the non-linear regression method developed in SPSS on the basis of the Newton-Lagrange approach (Alt and Malanowski 1993).
Hence, the fitted ST-LCM is the following: where the coregionalization matrices B 1 and B 2 are estimated as described hereafter.

Estimation of the coregionalization matrices
As clarified in Sect. 3.3, the estimation of the coregionalization matrices is the object of the third algorithm. In this case, the computation of the entries of B l , l = 1, 2, as indicated in (14), requires the values c 1 (0, 0) and c 2 (0, 0) of the two basic components, as well as the information coming from the sample covariance of the variables under study, i.e.
• c i j (0, 0), i, j = 1, 2, 3, • c i j (δ s , δ t ) 1 and c i j (δ s , δ t ) 2 , detected at the first and the second scales of spatiotemporal variability, i.e. at 10 km and 12 h, and at 13 km and 24 h, respectively, where i = 1 refers to NO 2 residuals, i = 2 to NO residuals and i = 3 to PM 10 residuals. By looking at Figs. 3 and 4, the above-mentioned values are easily identified and consequently the elements of the coregionalization matrices are immediately deter-mined. In particular, for the elements of B 1 it results: whose eigenvalues are non negative ensuring, in this way, the positive definiteness of the B l , l = 1, 2, and consequentially the admissibility of the fitted ST-LCM, which is now completely specified: (20)

Adequacy of the fitted ST-LCM
A way to evaluate the adequacy of the fitted ST-LCM consists in exploring the results of the cross-validation procedure which produces leave-one-out predictions on the basis of the fitted model and the data set under study. At each point of the domain, the observed value is temporarily eliminated and estimated by using the fitted ST-LCM and the measurements available in a neighborhood properly defined. In the spatiotemporal multivariate case, cokriging estimations can be computed for cross-validation purposes through the routine "COK2ST" implemented in the GSLib package (De Iaco et al. 2010).
In this case study, the previously mentioned program has been executed in order to obtain leave-one-out predictions of the residuals for each variable in turn (NO 2 , NO  (20). Moreover, the suitability of the selected ST-LCM has been evaluated even for a decreasing number of observations; in other terms, the cross-validation procedure has been performed by considering, alternatively, three different input databases, composed of 100%, 80% or 60% of the data (residuals). Note that the cut data sets have been obtained from the original data set through a random selection. The Spearman correlation coefficients calculated between the residuals of the analyzed variables and the cross-validation estimates have been reported in Table 2, together with the relative mean absolute estimation errors (RMAE).
Note that the linear correlation coefficients are very high for all the analyzed variables (at 1% significance level), even when the reduced data sets are used. Indeed, the reduction of Spearman coefficients is at most of 3.4% for the cut data sets. The RMAEs are approximately constant for all the alternative databases and show that the relative discrepancies (between the observed and the estimated values) are on average less than 20%. Thus, these results can be used to confirm the adequacy of the fitted ST-LCM.

A comparative analysis
In this section, the proposed procedure has been compared with respect to the exiting fitting method, which consists in using the same class of covariance models for all the basic components, without considering the main features of the basic components that the statistical tests can reveal. To this aim, two other alternative ST-LCMs have been fitted, as specified below: • the ST-LCM based on the product-sum covariance model, which has been fitted either for the SDC and the LDC, i.e.
where P S c 1 (δ s , δ t ) is as in (16) and with k 12 = 17.015, k 22 = 18.572, k 32 = 40, C s 2 is the spatial exponential covariance model in R d and C t 2 the temporal exponential covariance model in R, with practical ranges a s 2 = 13 km and a t 2 = 24 h, respectively; • the ST-LCM, obtained by using the Gneiting class of model for both the SDC and the LDC, i.e.
Then the adequacy of the ST-LCM in (20), defined by using the proposed algorithm, has been compared to the one of the models (21) and (22), developed through the usual fitting procedure. First of all, the comparative analysis has been conducted by computing the discrepancies between the empirical and the theoretical surfaces associated to the direct and cross-covariance functions. To this aim, the relative mean difference ( ) has been determined and analyzed as the number of lags increases, starting from small spatial and temporal distances and including step by step larger distances. This index represents an adequate synthesis of the errors between the fitted model and empirical covariance surface and is defined as follows: . . . , m, (23) where c i j (δ s , δ t ) k are the values of the direct (i = j) and cross (i = j) sample covariance, while c i j (δ s , δ t ) k are the values of the fitted covariance model at the K fixed lags. Table 3 shows the values of index (23) for each fitted ST-LCM, computed for an increasing number K of lags, as specified below: • K = 52, when the first 4 and 13 lags in space and time are analyzed; • K = 150, when the first 6 and 25 lags in space and time are considered; • K = 296, when all lags fixed in space and time for the structural analysis are examined.
Thus, from the comparison between the empirical and the theoretical covariance surfaces, the , obtained on the basis of the proposed algorithm, has performed better than the other two models, with an improvement (in terms of reduction of the index ) of 8% on average. It is worth pointing out that: • the improvement, with respect to the model (21), rises up to 13.3% for K = 296; thus the discrepancy increases when large spatial and temporal distances are included in the computation of the index , since in model (21), the LDC was forced to be modelled by using the product-sum class;  • the improvement, with respect to the model (22), rises up to 14.1% for K = 52; thus the discrepancy increases mainly for small spatial and temporal distances, since in the model (22), the SDC was forced to be modelled by using the Gneiting class.
In addition, the superiority of the ST-LCM in (20), coming out from the proposed algorithm, has been also evaluated on the basis of the cross-validation results. In particular, the leave-one-out estimates of each variable under study (NO 2 , NO and PM 10 ) have been also computed through cokriging, by using the models in (21) and (22). As done in Sect. 4.5, the cross-validation technique has been performed by considering, alternatively, three different input databases, composed of 100%, 80% or 60% of the original residual data. Then the RMAEs have been determined and compared with respect to the ones of the . From the results given in Table 4, it is worth pointing out that the RMAEs associated to the  are significantly lower than the ones of the other models. Indeed, the percentage decrease varies from 16 to 28%: in detail, it is in the range 15.56-16.20% when the complete data set is used, or in the range 20.42-21.07% when the 80% of the data are considered, or in the range 25.39-27.65% for the smaller data set.

Summary
This paper aimed to present some computational advances regarding the selection of an ST-LCM to be used for modeling the joint spatio-temporal evolution of environmental variables. It is also worth pointing out that the flexibility and good reliability of the developed procedure support an easy reproduction of all stages for data sets concerning either the same information but updated in space and time, or different variables. Moreover, the algorithms can be performed with the support of ad-hoc R functions (R Core Team 2015), thus practitioners can use this procedure for other applications. For this reason, as a future development, the effort should be concentrated on the construction of an R package.