The cgeostat Software for Analyzing Complex-Valued Random Fields

Given a vectorial data set in two dimensions, a representation on a complex domain is often convenient. This representation is rarely considered in geostatistics, although interesting applications can be found in environmental sciences and meteorology (e.g., for wind ﬁelds). In such a case, some computational diﬃculties are related to the lack of software for estimating and modeling a complex covariance function, for predicting complex variables as well as for representing the output results. In this paper, the new Fortran software cgeostat for geostatistical analysis of complex-valued random ﬁelds is presented and an application is demonstrated.


Introduction
Complex-valued random fields are often appropriate to model vectorial data with two components: This kind of data can be found in many areas of research, such as in hydrology and meteorology (i.e., in wind or stream studies). The formalism of complex-valued random fields can reflect the specific characteristics of the components of a vectorial random field: These components are associated with a physical phenomenon described by homogeneous quantities, characterized by the same unit of measure, such as wind velocity, force, electric or magnetic field, which are available at the same spatial points over the domain (De Iaco, Maggio, Palma, and Posa 2012). In this case, the corresponding realization of a vectorial field in two dimensions is an expression of a single entity, i.e., a complex number, where the decomposition in modulus and angle is natural and has a physical interpretation. Thus, a compact and unified formalism has to be suitably adopted for this kind of data. In other terms, the use of a complex formalism and consequently of complex kriging for prediction purposes is not an arbitrary choice, but it is strictly motivated by the nature of the phenomenon. Nevertheless, geostatistical applications are usually concerned with scalar variables or only with the magnitude of vectorial variables, such as the norm of the wind vector (Haslett and Raftery 1989;Cressie and Huang 1999).
The specific formalism of complex-valued random fields for vectorial data with two components, as well as various formulations of the predictor for bivariate random fields with a representation on a complex domain, were proposed by Lajaunie and Béjaoui (1991) and Wackernagel (2003). In particular, Lajaunie and Béjaoui (1991) applied the complex kriging and spent some effort to describe valid covariance functions. They also compared the complex kriging to simple kriging or direct cokriging of the two components of the vectorial variable. Some other contributions concern admissible complex covariance models, such as the one given by De Iaco, Palma, and Posa (2003), where the complex covariance models are generated by starting from a real covariance model and its spectral representation, and by introducing a shifting factor in order to translate the spectral density function.
Even if a lot can still be done from a theoretical point of view, computational advances in complex geostatistical analysis are also in great demand. In particular, for a wider use of the geostatistical theory of complex-valued random fields, some computational aspects on complex prediction should be developed.
In the literature, 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)
In this paper, after introducing the framework of complex-valued random fields, a discussion on complex kriging and some details on complex covariance models are presented (Section 2). Then, the new software cgeostat, to be applied for spatial analysis of phenomena with a representation on a complex domain, is proposed (Section 3). This is a software package which has been written in Fortran by using the Geostatistical Software Library GSLib, given in the version called GSLIB77. Note that these GSLib source codes can be freely downloaded from the website http://www.statios.com/Quick/gslib.html and their availability allows the basic algorithms to be used for more advanced customized programs. The cgeostat software package presented in this paper includes the program ckriging for predicting complex vari-ables and some other programs for: representing complex data (clocmap), computing and plotting components of a sample complex covariance function (ccovaplt), fitting and computing components of a complex covariance model (cfitting and ccovamodel). The program ckriging provides various options for prediction, such as complex kriging for a grid of points or blocks, cross-validation and jackknife; moreover the user can choose between simple kriging, if the expected value is known, and ordinary kriging, if the expected value is unknown. In Section 4, the new software is applied to study a complex-valued random variable and to provide predictions of the same variable at some unsampled locations on the basis of the available data. Note that cross-validation and jackknife options are also considered and some descriptive and inferential statistics regarding prediction performance are discussed.

Theoretical framework on complex-valued random fields
Let U and V be two real-valued random fields, which represent the components of the complexvalued random field, defined as follows where i = √ −1 and d is the dimension of the spatial domain (d ∈ N + , d ≤ 3). The random field W is stationary of second order, if its expected value, denoted with m, exists and is constant for any u over the domain and the covariance function C depends on the lag separation vector h, for any two points u and u + h, that is where the complex conjugation is denoted with ·. The above complex covariance function can be explicitly written as, where the real part is defined by the sum of the covariances of U and V , i.e., C U (h) and C V (h), respectively, while the imaginary part is given by the difference between the related cross-covariances C VU (h) and C UV (h), Note that the real part of (2) is an even function and, because of the convexity property (Chilés and Delfiner 1999), it is a covariance function, whereas the imaginary part of (2) is an odd function and is equal to zero in presence of isotropy and cannot be considered as a covariance function.
Moreover, for a complex-valued random field W , which is stationary of second order, the positive definiteness condition of C is given by: for any n ∈ N + and any choice of points u i ∈ R d , and a i ∈ C, i = 1, . . . , n.

Some details on complex covariance models
For prediction purposes, it is necessary to model C Re (defined as the sum of the covariances C U and C V ) and C Im (defined as the difference of the two cross-covariances C VU and C UV ) in a consistent way such that C is admissible, i.e., that it is a complex covariance function. In particular, De Iaco et al. (2003) showed that is a family of complex covariance functions, whereC(h) is a real covariance function and c ∈ R d , c = 0. In other words, given a real covariance functionC(h) and the corresponding even spectral density function f (s), complex covariance models can be obtained by shifting the function f (s), i.e., given a vector c = 0, the function f (s−c) is a not-even spectral density function ∀c ∈ R d , c = 0, hence, the spectral representation based on this translated function can be used to generate a complex covariance model.
In this paper, the above complex covariance models are used for modeling the correlation structure of the variables under study, since: • the family of complex covariance models in (4) is generated by starting from real covariance functions, thus various examples of complex covariance functions can be easily obtained by using such a model; • the fitting procedure is shown to be flexible enough.
As described in De Iaco, Posa, and Palma (2013), the family of complex covariance functions, given in (4), can be fitted by using an easy step by step procedure. After estimating the covariance functions for U and V and the corresponding cross-covariances for different lags and directions of the spatial domain (step I), the presence of anisotropy with type and direction of maximum continuity can be detected (step II). Then the components of the complex covariance function in (2), i.e., the sum of the sample direct covariances and the difference of the sample cross-covariances for relevant directions, have to be computed (step III). At last, the shifting vector c has to be estimated, for example through least squares regression (step IV), and an appropriate model forC can be easily chosen by inspecting the real part of the sample complex covariance function (step V).

Prediction of complex random fields
Let W be a second order stationary complex-valued random field, defined as in (1). Prediction of W at a location u 0 can be performed through different forms of kriging, which can be considered as an extension of the ones known for random fields on a real domain: The simple complex kriging and the ordinary complex kriging can be used depending on whether the expected value of W is known or not, as described below.
Let the predictor of the complex random variable W (u 0 ) at an unsampled location u 0 be where W represents the vector [W (u 1 ), W (u 2 ), . . . , W (u n )] of the complex random variables W (u α ) at the sampled locations u α and ω corresponds to the vector of the complex weights [ω 1 , ω 2 , . . . , ω n ] . In terms of the real and imaginary parts, predictor (5) can be written as where ω Re and ω Im are the vectors of the real and imaginary parts of the vector ω, while Under the unbiasedness condition, if it is required that the error variance is minimum, then the following complex ordinary kriging system, based on (2n + 2) equations and (2n + 2) parameters, is obtained where , is known or separately estimated over the study area, the residuals U and V of the real and imaginary parts can be obtained as follows Then, the residual complex random field W (u) = W (u) − m(u) can be defined as Hence, complex simple kriging of W , based on 2n equations and 2n parameters, can be performed by solving the following system where C Re = [C Re (u α − u β )] (n×n) and C Im = [C Im (u α − u β )] (n×n) , with α, β = 1, 2, . . . , n, while c Re = [C Re (u 0 − u α )] (n×1) and c Im = [C Im (u 0 − u α )] (n×1) , with α = 1, 2, . . . , n. Note that the above components of the complex covariance function C are related to the variable W .
As will be clarified in the following section, the complex covariance model in (4) and the above systems (6) and (7) have been suitably implemented in the new software cgeostat.

Remarks
Some general aspects on complex formalism for bivariate data and a theoretical comparison of complex kriging with respect to other techniques (such as separate kriging and complex cokriging) were proposed by Wackernagel (2003). Thus, any empirical comparison would reflect the already known theoretical results, as highlighted in De Iaco et al. (2013). As pointed out in the introduction, in addition to these aspects, the difference between complex formalism and other (univariate or multivariate) alternatives is not only formal, but conceptual, since any realization of a vectorial field in two dimensions (such as a wind field, streams, an electric or magnetic field) is expression of a single entity, i.e., a complex-valued function; for this reason, the compact and unified complex formalism is suitably adapted for this kind of data.
From a theoretical and practical point of view, the advantages of using complex formalism, whenever the phenomenon requires this formalism, are described below: • Complex kriging is adequate for predicting vectorial data if the components, characterized by homogeneous quantities, are correlated and the decomposition in modulus and angle is physically meaningful. An interesting discussion on the rotation and the scaling effects produced, to vectorial data, by the complex kriging weights, can also be found in De Iaco et al. (2013).
• The complex covariance completely characterizes the correlation structure of vectorial data in two dimensions; this enables the computation of spatial correlations in an easier and more flexible way compared to bivariate techniques.
• Using complex kriging simplifies the prediction step, since the complex kriging system depends on 2n parameters (where n is the number of available data), while the "complex cokriging" system, which is equivalent to the well-known formalism used for the cokriging predictions of a real-valued random vector or multivariate random field of two components (bivariate cokriging), measured at the same sample locations (isotopic case), depends on 4n parameters. If bivariate cokriging is used, it is not enough to have a model for the sum and the difference of direct and cross covariances, but a valid model for every direct and cross covariance function is also required.

Software for complex random fields analysis (tutorial)
The principal problem when analyzing bivariate data, with a reasonable representation on a complex domain, is the lack of software. Thus, customized GSLib (Deutsch and Journel 1998) Fortran programs for complex data representation (clocmap), estimation and plotting components of a sample complex covariance function (ccovaplt), fitting and computing a complex covariance model (cfitting and ccovamodel), complex kriging prediction (ckriging) have been implemented and are presented in this paper.
Note that the parameters (such as variables to be used, input/output files, options) have to be set by the user through an ad-hoc file, named parameter file (see Appendix A for details).
These Fortran programs, together with the GSLib library (see Appendix B for details on dependencies), and the corresponding executable programs are freely available, moreover a directory with the parameter files, the input and output files used for the example given in this article are also included.

Program clocmap
The program clocmap can be used to plot a location map for a complex-valued variable, measured at some spatial points over a two-dimensional domain. The data values for the two components U and V are plotted by using the usual representation on a complex domain, that is each pair of values is illustrated with an arrow whose modulus and direction depend on U and V . If the program clocmap is applied to represent the complex kriging estimates for U and V , then their error variance (with the option to apply some transformation, such as the square root) can be also illustrated with a circle whose diameter is proportional to the magnitude of the error variance. Moreover, the overlapping of the location map of the complex data used to computed the estimates (called conditioning data) can be also considered. Different colors for the location map of the estimates and the location map of the conditioning data can be chosen. Note that if the data set is referred to a 3D domain, a 2D representation can be obtained by using an option in the program clocmap, where the user can specify the 2D section, for a given coordinate z, to be plotted; in this way, the program can be used to produce different profiles of the data corresponding to various sections.

Parameter, input and output files
In the following, the input and output files of the program clocmap and the associated parameter file are described.
There are two input files: • The first file is mandatory and represents the data file with at least the two components of the complex data (and if necessary a third variable) for each location, in Geo-EAS format.
• The second one is optional and represents the data file with the conditioning data (sample data for the two components, which have been used for structural analysis), in Geo-EAS format; indeed, this is required only when the overlay option is set equal to 1.
These files are ASCII files (or text files), their extension is usually .dat and their names can be arbitrarily fixed in the parameter file.
Regarding the output, the program clocmap creates a postscript file.
An example of a parameter file for the program clocmap is given below.
Parameters for clocmap *********************** START OF PARAMETERS: In particular, apart from the parameters which are similar to the ones that are used in GSLib (Deutsch and Journel 1998), it is worth noting that: • the parameter that specifies the color code might take an integer number from 1 to 8 (1 = red, 2 = yellow, 3 = green, 4 = light blue, 5 = dark blue, 6 = violet, 7 = white, 8 = black); • the scaling parameter might be fixed equal to 0, if no scaling transformation is applied to the components U and V , while it is fixed equal to 1 if the logarithmic transform is applied; • regarding the parameters that specify the size of the symbol according to the modulus scale, if data minimum (dmin) is set equal or greater than data maximum (dmax), minimum and maximum are found from the data. Moreover, if minimum size (smin) is set equal or greater than maximum size (smax), minimum size and maximum size are set equal to the data minimum and maximum; • if the overlay option is set equal to 1, then the parameters regarding color, scale transformation and symbol size for the additional variable (usually the error variance, if available) can be appropriately fixed.

Run clocmap
If the parameter file has been appropriately set, the user can run the program clocmap. After typing the name of the parameter file, the program starts working.
The output is saved as a postscript file, whose name is defined in the parameter file. The ghostscript program would be a useful supplement for visualizing the graphical output.
Besides the error alerts which are common to all programs, specific error messages for the program clocmap may be issued: • to point out if the columns indicated for the coordinates and the components U and V are missing or not valid, • to detect if the number of data to be stored exceeds the maximum limit (which is equal to 10000).

Program ccovaplt
After computing the direct and cross directional sample covariance functions by using the available program gamvm of the GSLib library, it is necessary to determine and plot the components of the sample complex covariance function along different directions. Then the program ccovaplt can be applied.

Parameter, input and output files
In the following, the input and output files of the program ccovaplt and the associated parameter file are described.
First of all, the input files might be: • the direct and cross sample covariance functions, determined for different directions, in the format produced by the output of the program gamvm of the GSLib library; • the components of the complex covariance model, in the format produced by the output of the program ccovamodel (details of this program are given in Section 3.4).
These files are ASCII files (or text files) and their names can be arbitrarily fixed.
The output of this program is a postscript file and its name is also specified in the parameter file.
The parameter file required by the program ccovaplt is given below.
• for each covariance function to plot, the following parameters have to be set: the serial number that identifies the covariance to plot, selected among the covariances recorded (collected) in the indicated "file with covariance data"; the style options for the line: dashed line (integer numbers from 1 to 10, according to the desired length of the dash), point between dashes (1 = yes, 0 = no), continuous line (1 = yes, 0 = no) and color (integer numbers from 1 to 8, according to the color codes which are 1 = red, 2 = yellow, 3 = green, 4 = light blue, 5 = dark blue, 6 = violet, 7 = white, 8 = black).

Run ccovaplt
If the parameter file has been suitably specified, the user can run the program ccovaplt and type the name of the parameter file.
When the program has terminated, the output is saved as a postscript file, whose name is defined in the parameter file.
Moreover, the program ccovaplt produces a specific error message to acknowledge if the number of covariances exceeds the maximum value (equal to 10).

Program ccfitting
After computing and representing the components of the empirical complex covariance function along different directions, it is useful to estimate the translating vector c through nonlinear least squares. The program cfitting can be applied to this end.
Thus, the input files are the direct and cross sample covariance functions, determined for different directions, in the format produced by the output of the program gamvm of the GSLib library.
The output of this program is an ASCII file (or text file) and its name is arbitrarily defined in the parameter file.
The parameter file required by the program cfitting is given below.

Program ccovamodel
Another program is ccovamodel which calculates the components of a complex covariance model at certain lags and writes the result in a format compatible with ccovaplt. The parameter file required by this program is illustrated below.
Parameters for ccovamodel ************************* The output of this program contains alternatively the real component or the imaginary component of the complex covariance model in (4) according to the specification given at the second line of the parameter file, which is set equal to 1 for the real part computation, or equal to 2 for the imaginary part computation. Moreover, apart from some parameters that are similar to the ones used in GSLib (Deutsch and Journel 1998), such as the ones regarding the lags and directions where to calculate the components of a complex covariance model, the following details are also given: • the parameter nst, which corresponds to the number of structures ofC, used in (4); • the parameter associated with the nugget effect; • the parameters c 1 ,c 2 ,c 3 , which represent the components of the translating vector c; • the parameter it, which corresponds to the model type of each structure (Table 1); • the parameter cc, which corresponds to the contribution of each structure; • the parameters ang1,ang2,ang3, which correspond to the spatial angles, and a_hmax, a_hmin,a_vert, which correspond to the ranges or scales of variability of each structure.

Run ccovamodel
After specifying the parameter file, the user can run the program ccovamodel and type the name of the ad-hoc parameter file.
The output is saved in a text file whose extension (usually .var) is defined in the parameter file.
In the program ccovamodel, some specific error messages may be issued: • to acknowledge if the indicated model component does not respect the prefixed code; • to point out if the number of structures is not valid (that is equal or less than zero).

Program ckriging
The program ckriging can perform simple and ordinary complex kriging. According to what has been presented in Section 2, the program ckriging allows the system of equations (7) Type of covariance Integer code model real part imaginary part Spherical 1 5 Exponential 2 6 Gaussian 3 7 Hole effect 4 8 Table 1: Types of covariance models and the corresponding integer codes.
to be used in case of simple complex kriging. Alternatively, in case of ordinary complex kriging, the system of equations (6) can be used. The solutions of the systems (6) and (7) with respect to the parameter vector ω and the Lagrange multiplier M (in case of complex ordinary kriging) require the knowledge of the complex covariance model. The parameters, required by the program ckriging for setting the correlation structure given in (4) and the kriging options, are clarified below.

Parameter, input and output files
In the following, the input and output files of the program ckriging and the associated parameter file are described.
First of all, the input files are: • the data file which contains information on the sample locations and the observed values for the components U and V , in Geo-EAS format; • the file for jackknife estimation in Geo-EAS format, in case the jackknife option is set. The jackknife file contains information on the locations where the estimation is required and the available true values for the components U and V to be used for comparison.
These files are ASCII files (or text files), their names can be arbitrarily fixed and their extension is usually .dat.
The output files are: • the file which contains some descriptive and inferential statistics for the true data values (if available) and the estimates. Moreover, in case of cross-validation and jackknife, a statistical test on the differences between the means of the true values and the estimated ones for the components U and V is conducted and the p values for the null hypotheses that the differences are equal to zero are computed; • the debug file, with the usual GSLib debugging levels; • the file which contains the estimated values for the two components U and V and the estimation variance. Moreover, if the kriging option is set to cross-validation or jackknife, the true values for the components U and V and the estimation errors for the two components are also indicated in the output file, as shown in the output example given below. The parameter file required by the program ckriging, in order to provide cross-validation and spatial predictions of a complex random field, is illustrated below.

COMPLEX
Parameters for ckriging ***************************  Note that some parameters are clearly explained through the brief documentation reported on the right side of the parameter file; on the other hand, some specific parameters, required by the program ckriging, are described in order of appearance in the following list: • the kriging option: koption = 0 for complex kriging for a grid of points or blocks, koption = 1 for cross-validation with the data in datafl, koption = 2 or koption = 3 for jackknife with the data indicated in the jackknife file. If some data locations coincide with some points in the jackknife file, then the option 2 of jackknife allows the coincident data to be removed one at a time and estimated from the remaining neighboring data as in cross-validation, while the option 3 of jackknife allows the coincident data to be included in the searching neighborhood (this option produces exact estimates); • the name of the output file with the statistics regarding the available data and the estimated data; • in the subsequent 3 lines, the definition of the grid in an ad-hoc coordinate system, which is necessary for grid predictions (koption = 0); this is done in terms of the coordinates at the center of the first node/block (xmn, ymn, zmn), the number of grid nodes/blocks (nx, ny, nz) and the size/spacing of the nodes/block (xsiz, ysiz, zsiz). More details are given in the GSLib manual (Deutsch and Journel 1998); • in the last 6 lines, the kriging type and the complex covariance model: the kriging type is set to 0 if simple kriging is performed; it is set to 1 if ordinary kriging is performed; in case of simple kriging, the known mean values for the two components must be defined; the specifications of the complex covariance model in (4) can be provided for the real part or alternatively for the complex part, since given the translating vector c, the complex part can be defined starting from the real part and vice versa; hence, the indicator is set to 1 if the specifications are given for the real part, while it is set to 2 if the specifications are given for the imaginary part; the parameter nst corresponds to the number of structures ofC, used in (4); the parameter c0 corresponds to the nugget effect; the parameters c1, c2 and c3 represent the components of the translating vector c; the parameter it corresponds to the model type of each structure (Table 1); the parameter cc corresponds to the contribution of each structure ofC; the parameters ang1, ang2, ang3 correspond to the spatial angles and a_hmax, a_hmin, a_vert correspond to the ranges or scales of variability of each structure ofC.
In Table 1, the integer codes that can be used to specify the covariance models for theC in (4) are provided. It is worth highlighting that the admissible complex covariance models, generated through Equation 4, have been implemented in the subroutine ccov3 which is called by the program ckriging. Note that the same subroutine ccov3 is also called by the program ccovamodel described in Section 3.4.

Run ckriging
If the parameter file has been suitably filled, then the user can run the program ckriging and type the name of the prepared parameter file. The parameters set in the ad-hoc file are read and simultaneously displayed on the screen, then the input data files are read and saved in a temporary array, and some basic statistics are computed, as shown in the following template: As soon as the program starts the kriging procedure, the note Working on the complex kriging is shown; moreover, a report on the working progress from time to time is also provided on the screen, as given below: Working on the complex kriging currently on estimate 12 currently on estimate 24 currently on estimate 36 .... currently on estimate 456 At the end, some descriptive statistics on the estimates of the two components are shown on the screen and written to the debug file, then a closing message on the screen announces that the program has terminated, as shown in the following output template. The output estimates and some statistics regarding available data and estimated data are saved in two separate text files whose extensions (usually .out and .txt respectively) are defined in the parameter file.
As specified at the beginning, when the program ckriging is run, some error messages might appear. In particular, for this specific program, some alerts may be issued • to point out if the kriging option or the kriging type are not valid; • to detect if the indicated model type of a basic structure does not exist (it is not one of the models in Table 1).

An application on complex data
In the following case study, a vectorial data set with two components is interpreted as a finite realization of a complex-valued random field W , with components U and V (or modulus and direction in a polar system), and is analyzed by using the software cgeostat. A directory, called replication_materials, contains the executable files, the parameter files, the input and output files used for the case study of this article.

Data set
The whole data set of the spatial vectorial values (data19x24.txt) refers to 456 locations regularly distributed over a grid (19 × 24). However, a grid of (10 × 12) is selected and the corresponding vectorial data (data10x12.txt) is used for structural analysis; the remaining data is considered for comparison purposes in jackknife prediction.
The vectorial representation of the available complex data over the (19 × 24) grid is produced by applying the program clocmap with the parameter file clocmap19x24.par given in Section 3.1. On the other hand, the location map for the selected complex data is similarly obtained by applying the program clocmap with the parameter file clocmap10x12.par given below.

Structural analysis and modeling
Given the class of models in (4), a suitable complex covariance function for the data under study is found through • the estimation of the covariance functions for components U and V of the complex random field W , denoted as C U and C V , and the corresponding cross-covariance functions C UV and C VU for different lags and directions of the spatial domain; • the estimation of the components, C Re and C Im , of the complex covariance function given in (2), by computing the sum of the sample direct covariances C U and C V and the difference of the sample cross-covariances C VU and C UV for relevant directions, and fitting of the parameter vector c; • selection of an appropriate anisotropic covariance modelC(·; θ) with parameter vector The sample covariance functions C U and C V and the corresponding cross-covariance functions C UV and C VU for different lags and directions are estimated by using the original GSLib program gamv with the appropriate parameter files. In particular, four parameter files are filled in for directional sample covariance functions of the components U (covaU.par) and V (covaV.par), as well as for directional sample cross-covariance functions of UV (crosscovaUV.par) and VU (crosscovaVU.par).
The graphical representations of the estimated direct and cross covariance functions are obtained by applying the program ccovaplt with ad-hoc parameter files ( Figure 2). Hence, four parameter files are filled in for plotting the directional sample covariance functions of the components U (plt_covaU.par) and V (plt_covaV.par), and the directional sample cross-covariance functions of UV (plt_covaUV.par) and VU (plt_covaVU.par). For example the parameter file plt_covaU.par, which is used to plot the directional sample covariance function of the component U , is given below.
Then the real and imaginary components of the complex covariance are estimated and modeled. In particular, given the class (4) of complex covariance models, the least squares technique is first used for c 1 and c 2 estimation. The program cfitting with the parameter file given in Section 3.2, provides the least squares estimates for c 1 and c 2 . Then the following complex covariance model is fitted: whereC(h), withC(0) = 21.5, is a Gaussian covariance model with geometric anisotropy, whose direction of maximum continuity is 45 • with maximum range equal to 98 (corresponding to an effective range of 294) and minimum range in the perpendicular direction equal to 35 (corresponding to an effective range of 105), that is θ = (21.5, 294, 105), c = (−0.0025962, 0.0018629, 0.0). Note that the geometric characteristics ofC(h) are fixed through visual inspection of the empirical direct and cross covariance functions. In particular, the real and imaginary parts of the complex covariance model are given below: where h = h 2 x + h 2 y . The program ccovamodel is applied to compute the real and imaginary parts of the complex covariance model. The parameter file used to produce the real component is given in Section 3.4, while the parameter file used to determine the imaginary part is illustrated below.
Parameters for ccovamodel ************************* Both parameter files are filled in for computing the covariance model along 4 directions and 20 lags. Note that the model specifications fixed in these parameter files are also used for complex kriging (see Section 4.3). As specified before, the output of this program contains alternatively the real component or the imaginary part of the complex covariance model (4) in accordance with the specification given at the second line of the parameter file, which is set equal to 1 for the real part computation, or equal to 2 for the imaginary part computation.

START
The graphical representations of the sample directional components of the complex covariance function and the corresponding model are obtained by using the program ccovaplt. Four parameter files are appropriately filled in, one for each direction. For example, the parameter file plt_ccova_SE.par is illustrated in Section 3.2 and it is prepared so that the sample directional components (real and imaginary) of the complex covariance function and the corresponding model along the direction South-East are plotted. Note that, for a fixed direction, the sample real component is obtained by summing, lag by lag, the sample direct covariance functions for U and V , while the sample imaginary component is obtained by applying the difference, lag by lag, of the sample cross-covariance function for VU with respect to the one for UV . This reflects the theoretical equations given in (3). The reader can find the other three parameter files in the attached directory with the case study results. Figure 3 shows the real and imaginary parts of the sample complex covariance function and the corresponding models computed for different directions.
In Section 4.3, cross-validation, grid predictions and jackknife predictions are computed by applying the above model characterizations.

Complex validation and prediction
The program ckriging is used for both cross-validation and prediction purposes (by using the grid option or the jackknife options) with the parameter files appropriately filled in.
First of all, the complex covariance model for the corresponding random field is cross-validated on the basis of the available data at the 120 spatial points located over a grid (10× 12). The parameter file ckriging_cross.par used for cross-validation is presented below.
Parameters for ckriging ***************************  0.00 0.00 \ mean value(i), i = 1, 2 (nvar = 2) 1 \ model for real part (= 1), imag part (= 2) 1 0.0 -.0025962 0.001862 0.0 \ nst, c0, c1, c2, c3 3 21.5 45.0 0.0 0.0 \ it, cc, ang1, ang2, ang3 294.00 105.00 0.0 \ a_hmax, a_hmin, a_vert Note that • the complex kriging option is set equal to 1, in order to produce cross-validation; • the parameters related to jackknife predictions (such as the file name and the columns for the coordinates and the components U and V in the jackknife file) do not have to be specified in this case; • the file names which contain the debugging data, the cross-validation results and the corresponding statistics are set equal to complex_cross.dbg, complex_cross.out and statistics_cross.out, respectively; • regarding the grid specifications, they are useful in GSLib algorithms for defining local search neighborhoods apart from the selected kriging option. In this case, the coordinates at the first point are (53, 39, 0), the grid nodes are 10 along the East direction, 13 along the North direction and 0 in the vertical direction (indicating that a 2D grid is defined) and the spacing between the nodes are, respectively, 26, 26 and 1. The size of this search network is usually much larger than the final estimation or simulation grid node spacing (Deutsch and Journel 1998); • regarding block discretization, the value 1 results in point kriging being performed; • the characteristics of the search neighborhood are such that the estimation is performed by using from 3 to 8 sample data which fall in the ellipsoid whose maximum axis, along the direction of maximum continuity (corresponding to 45 • clockwise from the azimuth direction) is 200, and the minimum axis is 100. These values are fixed taking into account the ranges of the covariance model; • regarding the complex covariance model, the parameters are such that the fitted model (8) is considered.
Tables 2 and 3 illustrate some descriptive and inferential statistics for the original data values and the estimates, which are computed and saved in the output file statistics_cross.out. They can be compared and used to assess the estimation procedure and indirectly the fitted covariance function.
Indeed, a statistical test is conducted in order to verify that the differences between the means of the true values and the estimated ones for the components U and V are nil. The p values of the test-statistics are determined (Tables 2 and 3) indicating that the null hypothesis is retained. Hence, the complex covariance function can be considered appropriate for modeling the spatial correlation of the analyzed vectorial components, as also confirmed by computing the mean absolute error (MAE) and the root mean square error (RMSE) between the true values and the estimated ones for the two components (Tables 2 and 3).   Table 3: Some statistics on jackknife estimations with the kriging options 2 and 3.
After cross-validation, the new program ckriging is applied to predict the variable under study through complex ordinary kriging.
First of all, the grid option is considered; in particular, the above discussed covariance model is used for the estimation of the two components U and V over a grid (19×24). The parameter file used for the grid prediction is given below.
Parameters for ckriging *************************** Successively, the entire data set (not used for structural analysis) is also considered for the jackknife options 2 and 3. The two parameter files used for jackknife predictions, which are named ckriging_jack2.par and ckriging_jack3.par, are equal to the previous one except for the kriging option which is set equal to 2 for the former and 3 for the latter (the parameter file is illustrated in Section 3.5).
In Table 3, some descriptive and inferential statistics regarding the jackknife results, with both options, are given. It is worth underlining that the estimation results over the grid (19×24) obtained through the grid option (koption = 0) are equal to the ones obtained with the jackknife option 3 (koption = 3) over the same grid (19×24), except that this last option provides the comparisons with the true values.
The statistical test on the mean difference of the two components confirms that the null hypothesis of zero means can be accepted with low margin of uncertainty, in terms of p values. Moreover, note that both MAE and RMSE highlight the capability of the algorithm to reproduce the vectorial random field. Thus, the jackknife results show that the method can recover the unknown truth when every part of the model is correctly specified. Figure 4 shows the prediction map of the complex-valued random field for the study area, together with the error standard deviation (square root of the error variance), and the conditioning data (available data over the grid (10 × 12)), marked with red arrows. This location map is obtained by using the program clocmap with the parameter file clocmap_jack.par given below.

Summary
In this paper, some dedicated Fortran programs for geostatistical analysis of a finite realization of a complex-valued random field are presented. Complex kriging and some other programs, for representing complex data, computing and plotting components of a sample complex covariance function, computing and plotting components of a complex covariance model, have been implemented and an application on vectorial data with two components has been proposed. In particular, these routines are customized programs generated, on the basis of the GSLib library (specifically the version GSLIB77), for analyzing complex-valued random variables on a spatial domain R d (d ≤ 3).
The new GSLib programs ckriging, clocmap, ccovaplt, cfitting and ccovamodel, provide some computational solutions, which had not been implemented in the GSLib library before.
In particular, the program ckriging faces the problem of estimating a complex-valued random field. Differently from the original kriging program KT3D, the following computational aspects have been considered in the program ckriging: • complex kriging systems, given in (6) and (7), have been implemented; • the program COVA3, which belongs to the GSLib library, has been modified in order to include the admissible complex covariance model, defined in expression (4). This modified program has been renamed ccov3 and it has been included in the main program ckriging.for. In this way, it has not been necessary to modify the GSLib library; • simultaneous estimations of the two variables, which are the two components of the complex-valued random field, and the associated error variance are computed; • another kriging option for jackknife estimation has been included: the option 3 of jackknife allows the coincident data (the ones which coincides with the data locations) to be included in the searching neighborhood. This produces exact estimates associated with the jackknife points which are coincident with that data locations; • at the end, some descriptive statistics on the estimates of the two components are shown on the screen and written in a text file; moreover the mean absolute errors and the mean square errors for both components are computed, shown on the screen and written on the same text file.
Similarly, the other programs provide innovative tools for analyzing and representing complex variables and their correlation structure.
In the future, various complex covariance models and different forms of predictors for complex variables can be considered. It is worth highlighting that the implementation in the R environment and the extension to a spatio-temporal context might be also interesting as well as the idea of understanding if the complex-valued random fields models can adhere to solutions of particular partial differential equations should be investigated in the future.

A. Data files and parameter files
This appendix provides some GSLib standards in terms of the format of the data files and the parameter files required by the executable files.
There is no choice for the user to specify explicitly the input format (usually Geo-EAS format), however the accessibility of the source code allows this to be changed. Missing data are codified with large negative or positive numbers (e.g., −9999, or 9999, or −1.0e21 or 1.0e21).
It is important to highlight that the programs read numerical values and not alphanumeric strings. In the presence of alphanumeric variables, these can be transformed to integers or alternatively the source code can be modified.
All the variables, options and names of input/output files are contained in a parameter file.
In the parameter files, the parameter values can be appropriately stated on the left side, while a brief documentation of the parameters are reported on the right side. The user can add as many lines of comments at the top of the parameter file as desired, but formatted input starts immediately after the string "START". Example of parameter files are included in the paper. Each program requires its own parameter file, whose name is given by the user from a keyboard entry. Indeed, consistently to the GSLib standards, when one of the provided programs of the cgeostat software is run, the following question appears: Which parameter file do you want to use?
Then, the user has to type the name of the parameter file, previously filled, with its absolute path or relative to the current working directory. If no file name is typed, then the program would look for the parameter file which has the same name as the program and the extension .par (for example, the program ckriging would automatically look for ckriging.par). If the parameter file is not found a blank parameter file is created and the program stops.
Otherwise if the parameter file is found, the parameter values are read from the parameter file and are simultaneously displayed on the screen.
A closing message on the screen announces that the program has terminated, as shown in the following output template Stop -Program terminated.
The input and output files are saved in the folder the users are working in, unless a different path is specified in the parameter file.
When a serious error is encountered, an error alert is displayed on the screen and the program is stopped. Sometimes less serious problems or inconsistencies cause a warning to be displayed on the screen or in a debugging file, and the program will continue. In particular, when a program is run, some error messages might appear in order, for example, to acknowledge if the input files are missing or to indicate if the parameter file does not respect the prefixed standard.
However, it is possible to make a mistake that has not been foreseen. In that case, check the input and compare it to the data. If you cannot solve the problem, send the parameter file and data file to sandra.deiaco@unisalento.it.