Tractable circula densities from Fourier series

This article proposes an approach, based on infinite Fourier series, to constructing tractable densities for the bivariate circular analogues of copulas recently coined ‘circulas’. As examples of the general approach, we consider circula densities generated by various patterns of nonzero Fourier coefficients. The shape and sparsity of such arrangements are found to play a key role in determining the properties of the resultant models. The special cases of the circula densities we consider all have simple closed-form expressions involving no computationally demanding normalizing constants and display wide-ranging distributional shapes. A highly successful model identification tool and methods for parameter estimation and goodness-of-fit testing are provided for the circula densities themselves and the bivariate circular densities obtained from them using a marginal specification construction. The modelling capabilities of such bivariate circular densities are compared with those of five existing models in a numerical experiment, and their application illustrated in an analysis of wind directions.


Introduction
A direction observed in the plane R 2 can be represented by an angle, , typically in [0, 2π) or [−π, π), measured in a specified direction from a specified origin. The natural support for such directions is the circumference of the unit circle, S 1 , data on them being referred to as circular. The natural support for a bivariate circular random vector with angular coordinates ( 1 , 2 ) is the unit torus T 2 = S 1 × S 1 , data on such vectors sometimes being referred to as 'toroidal'.
Early constructions used to obtain models for toroidal data included maximum (Shannon) entropy characterization (or, equivalently, exponential family distributions, Mardia and Jupp 1999, p. 43), wrapping (Johnson and Wehrly 1977) and projection (Saw 1983). The latter two approaches can be applied to any bivariate distribution defined on R 2 , although the resulting toroidal densities generally cannot be expressed in closed-form or are highly convoluted.
Perhaps the best-known exponential family toroidal model is the eight-parameter bivariate von Mises distribution of Mardia (1975). Four of its parameters control the dependence between 1 and 2 , but not in easily interpreted ways. Its submodels include the six-parameter model of Rivest (1988) and the latter's five-parameter sine, cosine and hybrid submodels, the properties of which are studied at length in Mardia and Frellsen (2012). Such five-parameter models are, to some extent, toroidal analogues of the bivariate normal distribution, with four of their parameters controlling the locations and concentrations of the marginal distributions, and the fifth the dependence between 1 and 2 . However, their concentration and dependence parameters also control the unimodality/bimodality of their densities. Like the bivariate normal, the densities of their unimodal cases have contours that are elliptical around the mode. Their normalizing constants must be computed numerically. Their conditional distributions are von Mises, but their marginal distributions are generally not and, for some parameter values, can be bimodal.
More recently, Navarro et al. (2017) conditioned a multivariate normal distribution to obtain the twelve-parameter bivariate generalized von Mises model of order 2 (BGvM 2 ), whose conditional distributions are second-order generalized von Mises (GvM 2 ), see, for example, Gatto (2008), while Hassanzadeh and Kalaylioglu (2018) used a conditional specification approach to obtain a toroidal model with one of its marginal distributions being GvM 2 . GvM 2 densities can be symmetric or asymmetric and unimodal or bimodal. A highly flexible family of toroidal models obtained by normalizing time series spectra was proposed by Taniguchi et al. (2020). In general, the interpretation of the parameters of these various models is difficult, and their normalizing constants must be computed numerically.
Ameijeiras-Alonso and Ley (2022) extended the sine-skewing approach of Umbach and Jammalamadaka (2009) and Abe and Pewsey (2011) to generate models for mildly asymmetric toroidal data, focussing explicitly on sine-skewed extensions of the bivariate uniform, sine, cosine and wrapped Cauchy models. An appealing property of such models is that their normalizing constants are the same as those of the base models to which sine-skewing is applied.
In addition, c • is usually assumed to be continuous at (2π k, ψ 2 ) and (ψ 1 , 2πl) for any ψ 1 , ψ 2 ∈ [0, 2π) as well as all other points on the torus. Jones et al. (2015) provided an in-depth treatment of the Wehrly and Johnson (1980) class of circulas, with density where g is a circular density. See Jones et al. (2015) for details of the evolution of this class, Shieh and Johnson (2005) and Kato and Pewsey (2015) for special cases, and Pewsey and Kato (2016) for work on goodness-of-fit testing. The wider copula-related literature is summarized in Jones et al. (2015), and generalizations of copulas to other compact Riemannian manifolds have been considered in Jupp (2015); see also Jupp and Kume (2020).
Here we propose a general, Fourier series-based, approach to constructing tractable circula densities for use within (1) that includes (2) as a special case and focus on c • generated using five basic patterns of Fourier coefficients. The construction is attractive because: (i) various new circula models can be generated having simple closed-form expressions for their densities and flexible forms of dependence structure controlled by relatively few parameters; (ii) three well-known circular dependence measures are simple functions of just two of the Fourier coefficients; (iii) the conditional mean directions and mean resultant lengths depend only on a limited number of Fourier coefficients; (iv) methods for simulation, model identification, parameter estimation and testing goodness-of-fit and independence are available for both the c • and the bivariate circular distributions generated from using them within (1).
Although any forms of f 1 and f 2 can be employed within (1) to generate bivariate circular distributions, throughout the paper we illustrate the construction's application using marginal circular densities from the Kato and Jones (2015) family. The latter is the most flexible four-parameter family of unimodal circular densities presently available, has closed-form expressions for both its density and distribution functions, and all four of its parameters have clear interpretations. As we shall see, most of the conditional distributions of the specific circula densities considered in Sect. 3 are also members of the Kato-Jones family. We note that the c • that we explore here might also be used in a semiparametric approach to modelling in which kernel density estimates and empirical distribution functions are used in place of f 1 , f 2 , F 1 and F 2 in (1). Such an approach would inherit the tractable and interpretable parametric modelling of dependence.
We stress that the approach based on infinite Fourier series introduced here differs in important ways from previously proposed constructions incorporating truncated (or partial) Fourier series. Fernández-Durán (2007) considered bivariate circular models based on (1), the Wehrly and Johnson (1980) circula density (2) and non-negative trigonometric sum densities (i.e. truncated Fourier series constrained to be non-negative) for f 1 , f 2 and g. More generally, Pertsemlidis et al. (2005) and Fernández-Durán and Gregorio-Domínguez (2014) proposed toroidal densities obtained from truncated bivariate Fourier series, the latter constrained to be nonnegative. When fitted to toroidal data, such models generally include large numbers of parameters that are difficult to interpret (Mardia et al. 2007) and, due to harmonic artefacts, manifest modes, sometimes multiple, unsupported by the data (see, for example, Fig. 4 of Fernández-Durán and Gregorio-Domínguez 2014). In contrast, our tractable c • , as well as the toroidal models generated by using them within (1) in combination with the highly flexible unimodal circular marginal densities of Kato and Jones (2015), have relatively few parameters, all of which have clear interpretations, and closedform expressions for their densities which involve no computationally demanding normalizing constants. As an important consequence of the latter property, numerically implemented maximum likelihood estimation is swift. While truncated bivariate Fourier series with Fourier coefficients satisfying the constraints in (6) and (4) to follow might be used as circulas, they will generally not be as tractable as the ones based on infinite Fourier series that we propose here.
The rest of the article is organized as follows. In Sect. 2, we provide the details of our proposed construction together with general results for three circular dependence measures and the conditional mean directions and mean resultant lengths of circulas generated using it. In Sect. 3, we consider c • generated using five basic patterns of Fourier coefficients within the proposed construction and provide details of their basic properties. In Sect. 4, we explain how model identification, parameter estimation and goodness-of-fit testing can be performed, for both the considered c • and the models for toroidal data derived from using them within (1). The results from a numerical experiment, designed to compare the large-sample modelling capabilities of such toroidal models with those of five existing bivariate circular models, are presented in Sect. 5. In Sect. 6, a new Wehrly and Johnson (1980) circula density and our proposed inferential methods are applied in the analysis of wind directions. Lastly, we offer some concluding remarks in Sect. 7. All equations, figures and tables with numbers preceded by the letter S are contained in an accompanying online supplementary materials document.

Circula densities
It is well known (Mardia and Jupp 1999, Sect. 3.3.2) that any continuous circular density, f , can be expressed in the form of a Fourier series as where i = √ −1, for appropriately chosen Fourier coefficients φ(m) (m ∈ Z). Now consider the family of continuous distributions on the torus whose density can be expressed analogously as where the Fourier coefficients φ(m, n) ∈ C (m, n ∈ Z) are appropriately defined so that f (θ 1 , θ 2 ) ≥ 0 and π −π π −π f (θ 1 , θ 2 )dθ 1 dθ 2 = 1.

Proposition 1
The following hold for density (3).

(ii) A density in family (3) is a circula density if and only if
Proof (i) It is clear from the following well-known equation that (i) holds: (ii) If φ(m, n) satisfies (4), the marginal density f 1 (θ 1 ) can be expressed as The second equality follows from (5). Similarly, f 2 (θ 2 ) = 1/(2π), and thus (3) is a circula density. Next, assume that a density in family (3) is a circula density. Since this assumption implies that each marginal density of (3) is the circular uniform density, the equations in (4) follow from (i).
It follows from Proposition 1(i) that, for any φ(m, n) in (3), where z denotes the complex conjugate of z. The circula densities of Proposition 1(ii) can be expressed as If φ(m, n) = 0 for all m, n = 0, the density is that of the bivariate circular uniform distribution. Circula densities of the form (7) are the main focus of the paper.

Circular dependence measures
Here we provide general results for three signed circular dependence measures when applied to circula densities of the form (7). Let ( 1 , 2 ) be a bivariate circular random vector. Then, the dependence measures of Fisher and Lee (1983), Jammalamadaka and Sarma (1988) and, for a circula, that of Rivest (1982) are defined as respectively, where, for j = 1, 2, X j = (cos j , sin j ) T , μ j is the mean direction of j , and λ 2 denotes the smallest singular value of E(X 1 X T 2 ). For circula densities (7), it can easily be shown that the dependence measures between the circular uniform random variables 1 and 2 reduce to Clearly, when only one of φ(1, −1) and φ(1, 1) is nonzero, then ρ R , ρ JS and ρ FL are simple functions of either φ(1, −1) or φ(1, 1).
Analogous results hold for the conditional distribution of 2 | 1 = ψ 1 .

Preliminaries
In this section we consider examples of (7) generated by five simple patterns of nonzero Fourier coefficients and provide details of the basic properties of certain specific models. Given the constraints on the Fourier coefficients in (6) and (4) , unless explicitly stated otherwise, we consider circula densities defined through nonzero Fourier coefficients on the Z + × (Z + ∪ Z − ) lattice. Figure 1 illustrates the five patterns of nonzero Fourier coefficients. We emphasize the parametric specification of these coefficients. In particular, the specific models generated by Patterns 2-5 are derived using geometric series of nonzero Fourier coefficients which have the added benefit of generating circulas with closed-form expressions for their densities. In order to obtain circula densities exhibiting pointwise symmetry about the origin, all of the nonzero Fourier coefficients are assumed to be real. Since all of the patterns we consider have at most one of φ(1, −1) and φ(1, 1) which is nonzero, we only quote values of ρ R , those for ρ JS and ρ FL following from the relations identified in Sect. 2.2. The conditional mean directions and resultant lengths of all of the models we consider can be easily calculated using the results presented in Sect. 2.3; see Appendix A of the supplementary materials document. It follows that the conditionals of the models considered afford uniform, linear-homoscedastic, or nonlinear-heteroscedastic circular-circular regression. These properties and expressions for modes and antimodes of the circula densities we consider can be found in Table 1.
Appendix B of the supplementary materials document explains how data from any c • or bivariate circular density of the form (1) can be simulated.

Pattern 1: single (nonzero) point
Let Linear-homoscedastic Linear-homoscedastic Equation (11) Along Linear-homoscedastic Linear-homoscedastic where 0 ≤ γ ≤ 1/2 and, here and henceforth, q ∈ {−1, 1}. Thus, this pattern includes just a single nonzero real-valued Fourier coefficient. In this case, the circula density (7) reduces to Figure 2 portrays a planar contour plot of density (9) when q = 1 and γ = 0.3. For this model, and all but one of the other models in this section, ρ R = qγ (see Table 1). The strength of dependence between 1 and 2 is thus controlled by γ , and its sign by q. The conditional distributions of 1 | 2 = ψ 2 and 2 | 1 = ψ 1 are cardioid distributions on the circle.

Pattern 2: diagonal line
Consider circula densities generated using nonzero Fourier coefficients on a diagonal of the where the ϕ(m) are the Fourier coefficients of any circular distribution, for which ϕ(0) = 1. Then, the circula density is of the form where ] is a density on the unit circle. This is the Wehrly and Johnson (1980) class of circula densities discussed in Jones et al. (2015). For all but the case when g is circular uniform, (10) has linear contours parallel to the qπ/4 diagonal. Density (10) reduces to (9) when ϕ(1) = γ and ϕ(m) = 0 otherwise, and 0 ≤ γ ≤ 1/2. (10) can be expressed in closed form as For this model, all three parameters, γ = φ(1, −q) = ϕ(1), ρ and q, affect the dependence between ψ 1 and ψ 2 although, as elsewhere, ρ R = qγ . Figure S8 presents planar contour plots of density (11) for q = −1, γ = 0.7 and three values of ρ.
The parameter ρ regulates the concentration of the wrapped Cauchy-like distribution. Both conditional distributions are special cases of the unimodal circular distributions of Kato and Jones (2015) which are two component mixtures with circular uniform and wrapped Cauchy components.

Pattern 3: vertical line
Consider the pattern formed by nonzero Fourier coefficients on the vertical line of the Panels (a) and (b) of Fig. S3 present planar contour plots of density (12) with q = 1 and two (γ , ρ) combinations. The dependence between 1 and 2 is clearly regulated by all three parameters: q, γ and ρ. In particular, γ regulates the strength of the dependence between 1 and 2 , and ρ the degree of deformation of the density's shape around the main diagonal. The conditional distribution of 1 | 2 = ψ 2 is a cardioid distribution and that of 2 | 1 = ψ 1 is a special case of the Kato-Jones family.
Generalizations of Patterns 1-3, for which the nonzero Fourier coefficients are positioned more generically, are considered in Appendix C of the supplementary materials document.
Planar contour plots of circula density (13) designed to illustrate the roles of ρ 1 and ρ 2 are displayed in Fig. 3. When ρ 1 = ρ 2 , the density is symmetric about the main diagonal and increasingly concentrated in the neighbourhood of the origin as ρ 1 = ρ 2 increases. For a fixed value of ρ 1 , as ρ 2 increases the main axis of the central elliptical contour tilts increasingly away from the main diagonal towards ψ 2 = 0 and the dispersion increases in the neighbourhood of (−π, −π) = (π, π ). Due to the symmetry of (13), for a fixed value of ρ 2 the main axis tilts increasingly towards ψ 1 = 0 as ρ 1 increases.
Appendix D of the supplementary materials document briefly discusses the shapes of bivariate circular densities obtained using (1) with circula densities generated by Patterns 2 and 5.

Circula densities
Let {(ψ 1k , ψ 2k ); k = 1, . . . , n} denote an i.i.d. sample of random vectors from a circula density c • where, from here onwards, n denotes sample size. If the form of c • is unknown and n is moderate to large, an inspection of a planar scatterplot of the data will usually be sufficient to identify q and provide insight as to the form of c • . However, the absolute values of the empirical Fourier coefficients, φ(r , −qs) = 1 n n k=1 e i(r ψ 1k −qsψ 2k ) , for r , s = 1, 2, . . . , 6, say, represented graphically in a level plot, will generally prove more revealing, showing patterns, like those in Fig. 1, indicative of the structure of the Fourier coefficients of c • . Computation of the empirical Fourier coefficients is extremely fast. In practice, a range of potential c • 's might be explored and the best fitting model established using information criteria, model reduction based on likelihood ratio testing and formal goodness-of-fit testing.
All of the densities of the five specific circula models in Sect. 3 have, by design, relatively simple closed-form expressions involving no computationally demanding normalizing constants, so computation of their log-likelihood functions is straightforward. Maximum likelihood (ML) estimation is then conducted using standard constrained optimization techniques. Method of moments (MM) estimates, calculated sequentially using the relations in Table 2, q = sgn(|φ(1, −1)| − |φ(1, 1)|), and φ(r , −qs) substituted for φ(r , −qs), can be used as starting values.

Bivariate circular densities
The shapes of the bivariate circular densities obtained using density (1) depend on the reference points from which the marginal densities are integrated in the definitions of the marginal distribution functions F 1 and F 2 . Traditionally, the reference point used has been the origin, 0. However, for this choice, changes in the location parameters of the marginal distributions result in shape changes, not just location shifts, in the densities obtained using Eq. (1). To avoid such shape changes, we define F j (θ ) ( j = 1, 2) as where ω A j denotes the antimode of the circular density f j .
We advocate the following sequential approach to modelling i.i.d. samples of random vectors {(θ 1k , θ 2k ); k = 1, . . . , n} exploiting the three-component structure of density (1). If the distributional forms of the marginal densities, f 1 and f 2 , are not specified beforehand, histograms and/or kernel density estimates can provide insight into forms for them. Then, their parameters are initially estimated separately using ML. Denoting the marginal distribution functions corresponding to those parameter estimates byF s 1 andF s 2 , next the 'pseudo-sample' {(ˆ s 1k ,ˆ s 2k ); k = 1, . . . , n} = {(2πF s 1 (θ 1k ), 2πF s 2 (θ 2k )); k = 1, . . . , n} is computed. The procedures described in Sect. 4.1 are then applied to this pseudo-sample to obtain initial estimates of the form of the underlying circula density and its parameters. Finally, the estimates from the previous two stages are used as starting values in the maximization of the full loglikelihood function derived from (1). We denote the marginal distribution functions corresponding to the parameter estimates obtained in this final estimation stage byF 1 andF 2 .
Appendix E of the supplementary materials document explains how data from a circula density or a bivariate circular density of the form (1) can be tested for independence using the permutation approach proposed in Sect. 3.3 of Kato and Pewsey (2015).

Numerical experiment
As part of our investigations into the performance of the modelling approach based on (1) and our proposed circula densities, we performed an experiment designed to compare its large-sample modelling capabilities with those of five existing bivariate circular models. We simulated single samples of size n = 2000 from each of the six models in Table 3, ranging from the 6-parameter pointwise symmetric model of Rivest (1988) to models capable of describing very varied distributional shapes. The parameter values in the same table were selected to produce representative unimodal cases of each model.
For each simulated sample we fitted all six models by ML, without pursuing model reduction, using R code developed by us and the CircNNTSR package of Fernández-Durán and Gregorio-Domínguez (2016) to fit the bivariate non-negative trigonometric sum (BNNTS) models. Table 4 contains the results obtained, the AIC values for the BNNTS models being those for the BIC-adjudged best fits, BIC being the model selection criterion advocated by Fernández-Durán and Gregorio-Domínguez (2014).
Using the column-wise sum of the ranks of the off-diagonal row-wise BIC values in Table 4 as a simple measure of overall performance, the following ordering, from best to worst, is obtained, with the sums of the ranks appearing between square brackets:  [20]. Both orderings identify the BGvM 2 model as having the best overall performance. The BIC-based measure identifies the KJ 2 (13) model as having an overall performance similar to those of the BNNTS and FBvM models, with a marginally better overall performance than the R and SsS models. According to the AIC measure, BNNTS models perform second best, with the KJ 2 (13) model performing similarly to the FBvM model and, again, marginally better than the R and SsS models.
We note that, for the samples from all five alternative models, the best-fitting case of (1) with Kato-Jones marginal densities and one of our proposed circula densities (KJ 2 Circ) was the one incorporating circula model (13) (KJ 2 (13)). This fact indicates that the circula structure underpinning the five alternative models can be approximately matched using just one of the circula models considered. Figures S11 and S12 present planar scatterplots of the six simulated data sets with contour plots of the best-fitting alternative and KJ 2 Circ densities superimposed. They provide insight into the ability of the KJ 2 Circ models to mimic the alternative models, and vice versa.
As the number of parameters of the underlying model increases, the performance of the KJ 2 Circ and BNNTS models becomes increasingly competitive, although the latter generally have far more parameters than their KJ 2 Circ counterparts. Also, while Table 3 Densities and parameter values of the six models used in simulations: Rivest (R); sine-skewed Sine (SsS); full bivariate von Mises (FBvM); bivariate generalized von Mises of order 2 (BGvM 2 ); (1) with Kato-Jones f j (θ j ) and (14) as circula density (KJ 2 (14)); bivariate non-negative trigonometric sum of order (1,2) (BNNTS (1,2)). Throughout, c j = cos(θ j − μ j ), s j = sin(θ j − μ j ) and j = 1, 2. The models are ordered, from top to bottom, by increasing number of parameters

Illustrative application
In contrast with the large sample size employed in Sect. 5, our illustrative analysis considers a much smaller sample of n = 80 pairs of wind directions, {(θ 1k , θ 2k ), k = 1, . . . , n}, recorded daily at a Milwaukee meteorological station at 4am and 6am, respectively, during the last three months of 2020. (This analysis was inspired by a related data analysis in Wehrly and Johnson (1980), which we have updated and expanded.) The data were extracted from the vast Milwaukee Met Data Archive series for 2020 available at https://www.glerl.noaa.gov/metdata/mil/archive/mil2020. 01t.txt, and files containing them are included in the zip file linked to the supplementary Given the two-hour separation between the measurements in each pair, one would expect them to be correlated. However, as the pairs span a trimester, it is not necessarily obvious what form of relation might exist between the pairs, nor what shapes the marginal distributions might exhibit. For the univariate time series, the only sample autocorrelation coefficient identified as being significantly different from 0 using the randomization version of the approach of Fisher (1993, Sect. 7.2.2) was the lag 1 coefficient for the 6am series, with an estimated p-value of 0.03 based on 9999 randomizations. However, the value of that autocorrelation coefficient is low (0.11) and in our analysis of the pairs of observations we treat them as forming an i.i.d. sample of toroidal data. Figure 5 includes a planar scatterplot of the data converted to radians in [−π, π), and Fig. 6a linear histograms of the θ 1 and θ 2 values. In the latter, the superimposed circular densities are those of Kato-Jones ML fits.
The pseudo-sample {(2πF s 1 (θ 1k ), 2πF s 2 (θ 2k )); k = 1, . . . , n} is portrayed in Fig. 6b. From an inspection of it, the dependence between the pseudo-circular uniform variates is positive, i.e. q = 1. Figure 6c is a level plot of the absolute values of the empirical Fourier coefficients for that pseudo-sample and q = 1. The largest absolute values can be seen to form a diagonal pattern, indicative of circula density (11) as a potential model for the underlying c • . Figure 6d presents a planar scatterplot of {(2πF 1 (θ 1k ), 2πF 2 (θ 2k )); k = 1, . . . , n} superimposed upon a contour plot of the circula density from the fitted KJ 2 (11) model.   The point estimates from the different stages of the fitting process are given in Table 5. Table 6 contains the AIC and BIC values for the KJ 2 (11) model and the five alternative bivariate circular models employed in Sect. 5. Both criteria select the KJ 2 (11) model, a member of the Wehrly and Johnson (1980) class with Kato and Jones (2015) circular marginal distributions and circula density (11), as being best. Contour plots of all six fitted densities are presented in Fig. S13. The densities of the models with the highest BIC-values, SsS and BNNTS, are bimodal.

Conclusion
As the specific models introduced in Sect. 3 illustrate, the proposed general Fourier series construction provides a means of generating tractable parametric circula densities with simple closed-form expressions and wide-ranging distributional shapes. Moreover, when combined with Kato and Jones (2015) circular densities through (1) they provide highly flexible models for toroidal data for which the marginal distributions are unimodal. Multimodal toroidal data can be modelled using mixtures of such models.
The flexibility of the circula densities can be further enhanced by allowing the Fourier series coefficients to be complex. The main effect of such an extension is to skew the circula distributions in various ways.
Stationary Markov models for circular time series can be defined from our circula densities using an analogous approach to that of Wehrly and Johnson (1980).
In principle, our bivariate circula construction can be extended to produce ddimensional circula densities using the multivariate analogue of (3) and patterns of nonzero Fourier coefficients distributed in d dimensions. Another possibility is, as mentioned in Jones et al. (2015), to model multivariate circular data using the circular analogues of pair copulas.
Our proposed methodology makes use of level plots of the absolute values of empirical Fourier coefficients computed from pseudo-samples as a highly successful model identification tool. We stress again that such plots are generally easier to interpret than scatterplots of the pseudo-samples themselves because they usually provide insight into the structure of the Fourier coefficients of the underlying circula density.
The advantage of the modelling approach based on (1) is that it facilitates the separate modelling of the circular marginals and a circula density in a structured sequential way. Also, as we have illustrated, it accommodates formal goodness-of-fit testing, an issue neglected in the literature related to the application of existing models for toroidal data.

Electronic supplementary material
The online supplementary materials document contains the appendices, and the equations, tables and figures with numbers preceded by the letter S. It also includes a link to a zip file containing all the data sets and our R code developed to analyse them.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Conflict of interest
The authors declare that they have no conflict of interest.

Data and code availability
The data sets and the R code used to analyse them can be accessed via the link in Appendix G of the supplementary materials document.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.