ORIGINAL PAPER Changing correlations: a flexible definition of non-Gaussian multivariate dependence András Bárdossy1 Accepted: 13 February 2023 / Published online: 2 March 2023 � The Author(s) 2023 Abstract Dependencies between variables are often very complex, and may for high values, be different from that of the low values. As the normal distribution and the corresponding copula behave symmetrically for low and high values the frequent application of the normal copula for the description of the dependence may be inappropriate. In this contribution a new way of defining high dimensional multivariate distributions with changing correlations is presented. The method can also be used for a flexible definition of tail dependence. Examples of copulas with linear changing correlations illustrate the methodology. Parameter estimation methods and simulation procedures are discussed. A five dimensional example using groundwater quality data and another four dimensional one using air pollution data, are used to illustrate the methodology. Keywords Multivariate copulas � Tail dependence � Asymmetry 1 Introduction The dependencies between two or more variables, can be very complex. In the case of environmental variables physical, chemical and biological processes have a major influence on the variables of interest. These processes are often not explicitly understood and are non-linear. Statis- tical investigations of dependence are frequently based on correlations between the variables of interest. This however may lead to a sub-optimal recognition and description of dependencies. Multivariate distributions relating variables are often considered to be normal or normal after data transforma- tions. Popular transformations like the log transformation or Box–Cox transformations, transform the one dimen- sional marginals to normal, which does not imply that the bi- or multivariate marginal distributions also become normal. Relationships between variables are often more complex, even non-monotonic and frequently deviate from normal. An elegant possibility to describe complex relationships between continuous variables is possible by using copulas (Sklar 1959). Copulas enable the investigation of the dependence independently of the one dimensional marginal distributions. A large number of different theoretical cop- ulas are described in different publications such as (Joe 1997) or (Nelsen 1999). Copulas are frequently used in many different disci- plines such as finance, economy, environmental studies or hydrology. In hydrology they are used for describing multivariate flood frequencies (Gräler et al. 2013) rela- tionships between flood characteristics (Chen and Guo 2019) precipitation (Favre et al. 2018) or drought (Won et al. 2020) just to mention a few. In Bárdossy (2006) copulas were used for the spatial statistics for groundwater quality parameters. In Brunner et al. (2019) the Fisher copula was used to investigate the complex dependencies of flood occurrences. Some of the well known copulas are derived from known multivariate distributions such as the Gaussian or the t-distribution. Another way to construct complex dependencies is to use vine-copulas (Czado and Nagler 2022). While the vine copulas offer a very flexible way to describe dependence, their construction and application for high dimensional cases is relatively complicated. & András Bárdossy bardossy@iws.uni-stuttgart.de 1 Institute for Modelling Hydraulic and Environmental Systems, University of Stuttgart, Pfaffenwaldring 61, 70550 Stuttgart, Germany 123 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 https://doi.org/10.1007/s00477-023-02408-1(0123456789().,-volV)(0123456789().,- volV) http://crossmark.crossref.org/dialog/?doi=10.1007/s00477-023-02408-1&domain=pdf https://doi.org/10.1007/s00477-023-02408-1 The purpose of this paper is to introduce a very flexible family of distributions with value dependent (asymmetri- cal) dependence and varying tail dependence. The paper is divided into 8 sections. After the introduction the definition of copulas is presented, and a new family of multivariate distributions with non-Gaussian copulas is introduced. Parameter estimation issues and simulation procedures are presented. Theoretical examples illustrate the flexibility of the methodology. Two data sets, one five dimensional on groundwater, the other four dimensional on air pollution are used as real-life examples to demonstrate the method- ology. A short discussion and conclusions section com- pletes the paper. 2 Methodology 2.1 Copulas A copula C is defined as a multivariate distribution on the n dimensional unit cube: C : ½0; 1�n ! ½0; 1� ð1Þ which has to have uniform marginals: CðuðiÞÞ ¼ ui if uðiÞ ¼ ð1; . . .; 1; ui; 1; . . .; 1Þ Copulas are related to multivariate distributions through Sklar’s theorem (Sklar 1959): Each multivariate distribu- tion Fðt1; . . .; tnÞ can be represented with the help of a copula: Fðt1; . . .; tnÞ ¼ CðFt1ðt1Þ; . . .;FtnðtnÞÞ ð2Þ where FtiðtÞ represents the i-th one dimensional marginal distribution of the multivariate distribution. If the marginal distributions are continuous then the copula C in (2) is unique. Hence, copulas can be regarded as the pure expression of the dependence without the influence of the marginal distributions. Copulas can be defined explicitly using their density functions or another possibility is to define copulas using multivariate distributions by inverting (2). This means the copula is defined as: Cðu1; . . .; unÞ ¼ FðF�1 t1 ðu1Þ; . . .;F�1 tn ðunÞÞ ð3Þ Many well known copulas such as the normal, the t-copula and the skew-normal copula are constructed using this method. 2.2 Distributions with value dependent correlations In this paper a very flexible way of defining the multi- variate distributions which define copulas with interesting non-Gaussian properties is presented. The basic idea behind the construction of the distribution is to gradually change the dependence structure depending on the values of the variables. As described in Guthke and Bárdossy (2012) one can obtain very similar spatial fields if one uses the same set of random numbers to generate them. A similar methodology was used in Bardossy and Pegram (2012) to exchange correlation structures of simulated precipitation. This idea combined with continuity can be used to define random variables with changing dependence structure. Formally: Definition A matrix valued function F : R ! Rk�k is called continuous if each element i, j of the matrix FðsÞi;j is a continuous function Definition Let RðsÞ be a in s continuous function of cor- relation matrices of dimension k � k. Let Xs be a k-di- mensional normal random variable with the correlation matrix RðsÞ and standard normal marginal distributions. In this case XðsÞ can be coupled with the help of independent standard normal variables U ¼ ðU1; . . .;UkÞ in the form: XðsÞ ¼ RðsÞ 1 2U ð4Þ If the correlation matrices RðsÞ are continuous in s than U defines a set of interdependent XðsÞ random variables which are continuous in s. A random vector Z ¼ ðZ1; . . .; ZkÞ is defined as: Zi ¼ minfU�1ðsÞ ; ðRðsÞ 1 2UÞi ¼ U�1ðsÞg ð5Þ Note that as U�1ð0Þ ¼ �1\Xð0Þi and U�1ð1Þ ¼ þ1[Xð1Þi and XðsÞi is by definition a continuous function for any i and U, the above definition leads to well defined Zi-s. This k-dimensional vector variable has, by this defini- tion a dependence structure which is different for small and large values, resulting in a non-Gaussian multivariate distribution. Figure 1 explains the construction. For a given vector of random numbers ðu1; . . .; ukÞ and for two selected pairs of variables ði1; i2Þ and ðj1; j2Þ the corresponding ðXi1ðsÞ;Xi2ðsÞÞ and ðXj1ðsÞ;Xj2ðsÞÞ are plotted as a function of U�1ðsÞ. The ðZi1 ; Zi2Þ and ðZj1 ; Zj2Þ values correspond to the first intersection of the diagonal and the ðXi1ðsÞ;Xi2ðsÞÞ 2620 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 123 and ðXj1ðsÞ;Xj2ðsÞÞ functions. The ðX�ðsÞ functions are continuous due to the continuity of the function of corre- lation matrices. The intersection with the diagonal repre- sents the value assigned to Z1 and Z2. As the Figure shows the two different pairs get values from different XðsÞs. If RðsÞ ¼ R is constant for all s values, then the resulting random variable is multivariate normal with the correlation matrix R. 3 Construction of continuous R(s) functions For the definition of the random variables in equation (5) a continuous correlation matrix valued function is required. The construction of such a function is a non trivial task as these matrices should not have negative eigenvalues to be valid correlation matrices. 3.1 Construction using the square roots of the matrices In order to define such matrices one can use the fact that C is a covariance matrix if and only if it can be written as a product of a matrix Y and its transpose YT : C ¼ YT � Y . (Y is the square root of the covariance matrix.) If YðsÞ is a continuous matrix valued function then CðsÞ ¼ YðsÞT � YðsÞ is also a continuous matrix valued function. All these matrices are valid covariance matrices. By defining RðsÞ ¼ ðrðsÞi;jÞ from CðsÞ ¼ ðcðsÞi;jÞ as: ri;jðsÞ ¼ ci;jðsÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ci;iðsÞcj;jðsÞ q ð6Þ one obtains a continuous matrix valued function of corre- lation matrices due to the continuity of the transformation in (6). The general construction uses infinitely many matrices. For practical use one has to find simplifications. The most simple construct is to use two parameters for each pair of variables - a starting and a final correlation, called two- correlations model in the subsequent text. Formally: Let Rð0Þ and Rð1Þ two valid correlation matrices. If RðiÞ ¼ YðiÞT � YðiÞ for i ¼ 0; 1 then Y�ðsÞ ¼ sYð1Þ þ ð1� sÞYð0Þ ð7Þ can be used to define CðsÞ ¼ Y�ðsÞTY�ðsÞ covariance matrices. The correlation matrices RðsÞ corresponding to these covariance matrices form a continuous function of correlation matrices connecting Rð0Þ and Rð1Þ. This model is called the two correlation linear model. This construction can be generalized to produce a set of mþ 1 correlation matrices, the mþ 1 correlation model: Rðs0Þ;Rðs1Þ; . . .;Rðsm�1Þ;RðsmÞ with 0 ¼ s0\s1\. . .\sm�1\sm ¼ 1 . Then: Y�ðsÞ ¼ s� sj sjþ1 � sj Yðsjþ1Þ þ 1� s� sj sjþ1 � sj � � YðsjÞ if sj � s� sjþ1 ð8Þ defines the sequence of covariance, and a subsequent continuous sequence of correlation matrices. For any continuous function F of matrices the and continuous function h : ½0; 1� ! ½0; 1�, F(h(t)) is also a continuous function of matrices. This new matrix function defines a different random variable ZðhÞ with a different copula. The choice of the function h(t) can change the shape of the corresponding copula even if all correlation matrices remain the same. The left and right hand deriva- tives of the transformation function h at 0 h0ð0�Þ and at 1 h0ð1þÞ are responsible for the tail behavior of the corre- sponding copula. For those pairs where the starting Rð0Þ and or the ending Rð1Þ correlation matrix contain non diagonal values equal to 1, upper and/or lower tail dependence may occur. The value of the tail dependence can be anything between 0 and 1 by adjusting the speed of convergence of s to 1 or 0 respectively. The corresponding proof is in the Appendix of this paper. Fig. 1 Definition of the random variable Z: XiðsÞ functions for two given u1; u2 pairs and their intersection with the diagonal defining the values Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 2621 123 If all correlations are positive then another geometric construction could be used for the definition of the matrix function. 3.2 Construction using spatial statistics In this case the correlation matrices are constructed using curves in an m dimensional space. All curves are param- eterized with s in [0, 1] and the correlation matrices are defined using a stationary spatial covariance function. Let ðy1ðsÞ; . . .; ðykðsÞÞ be such that yiðsÞ is a point in an m\k dimensional space, and yiðsÞ is a line continuous in s. If C(h) is a valid continuous correlation function in space then: RðsÞ ¼ Cðy1ðsÞ � y1ðsÞÞ � � � Cðy1ðsÞ � ykðsÞÞ .. . . . . .. . CðykðsÞ � y1ðsÞÞ � � � CðykðsÞ � ykðsÞÞ 0 B B @ 1 C C A ð9Þ is a valid correlation matrix, and due to the continuity of the yiðsÞ-s the corresponding matrices are continuous in s. 4 Simulation The simulation of a realization of the model is quite simple. One only has to know the correlation matrices for each s and has to simulate k independent standard normal variables. The exact simulation can be done using the following procedure: 1. Select the starting Rð0Þ and ending Rð1Þ correlation matrices. 2. Draw k independent normally distributed random numbers ðu1; . . .; ukÞ 3. Solve the linear equation (refeq:regu) for si for each i ¼ 1; . . .k. Assign xi ¼ U�1ðsiÞ. The vector ðx1; . . .; xkÞ is a simulated member. 4. Repeat steps 2–3 N times This procedure is slightly slower than an approximate simulation using a discrete set of possible s values. For the simulation of the linear model (7) the following step-by- step procedure can be used. 1. Select the starting Rð0Þ and ending Rð1Þ correlation matrices. 2. Calculate the square roots of the starting and end correlation matrices. 3. Select a set of s values 0 ¼ sð0Þ\s1\. . .\sm ¼ 1. 4. Calculate the in-between correlation matrices for each si using linear interpolation of the square roots and renorming. 5. Draw k independent normally distributed random numbers ðu1; . . .; ukÞ 6. For each variable i find the simulated value is the xi for which: xi ¼ U�1ðsjÞ kU�1ðsjÞ � ðRðsjÞ 1 2ðu1; . . .; ukÞik minimal the vector ðx1; . . .; xkÞ is the simulated member. 7. Repeat steps 5–6 N times Both algorithms are very simple and large samples can be generated with little computational effort. 5 Examples As a first example consider a bi-variate distribution with correlation matrices Rð0Þ ¼ 1 0:2 0:2 1 � � Rð1Þ ¼ 1 0:9 0:9 1 � � ð10Þ with normal marginals and linearly changing correlation matrices according to (7). The Pearson correlation of the two variables (with normal marginals) is 0.60. For comparison the multivariate normal distribution with the same Pearson correlation (0.60) is considered. For both bi-variate distributions a sample of N ¼ 2000 was generated. Figure 2 shows the corresponding empirical copulas. One can see the effect of changing correlations leading to weak dependence for low values and a strong dependence for high values. Another example is a tri-variate distribution with two variables being independent and the third having a lower tail dependence with one of the variables and an upper tail dependence with the other one. The value of both tail dependencies can be arbitrary. This distribution can be constructed using the correlation matrices: Rð0Þ ¼ 1 0 1 0 1 0 1 0 1 0 B @ 1 C A Rð1Þ ¼ 1 0 0 0 1 1 0 1 1 0 B @ 1 C A ð11Þ Figure 3 shows the result of a simulation with N ¼ 3000 points. One can see that variables Z1 and Z2 are indepen- dent, while Z1 and Z3 have some kind of upper tail dependence and Z2 and Z3 have some kind of lower tail dependence. Such copulas are difficult to construct with other methods. Note that the exact value of the tail 2622 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 123 dependence is determined by the speed of convergence to the matrices Rð0Þ and Rð1Þ. To show the high flexibility of the copulas obtained by this construction two three-correlations examples in the sense of (8) are shown. They correspond to the three cor- relation matrices: Rð0Þ ¼ 1 0:95 0:95 1 � � R 1 2 � � ¼ 1 0:65 0:65 1 � � Rð1Þ ¼ 1 0:95 0:95 1 � � ð12Þ and Rð0Þ ¼ 1 0:65 0:65 1 � � R 1 2 � � ¼ 1 0:95 0:95 1 � � Rð1Þ ¼ 1 0:65 0:65 1 � � ð13Þ The first shows a relationship weaker than a normal dependence for the medium values. For the second the relationship is reversed; the medium values show a stron- ger dependence than the extremes (Figure 4). The construction of similar examples in higher dimen- sions is not difficult, one only has to be careful that the Fig. 2 Simulated empirical copulas corresponding to the two-correlations model (q(0)=0.2?q(1)=0.9) with Pearson correlation = 0.60 left, and corresponding to the normal copula with the same Pearson correlation (0.60) Fig. 3 Simulated empirical copulas corresponding to the 3 dimensional two-correlations model with Rð0Þ ¼ 1 0 1 0 1 0 1 0 1 0 @ 1 A Rð1Þ ¼ 1 0 0 0 1 1 0 1 1 0 @ 1 A Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 2623 123 correlation matrices corresponding to the m breakpoints are chosen so that they are positive semidefinite. 6 Parameter estimation The parameters of the distributions depend on how the changing correlation structure of the random variable is chosen, as there is no general analytic expression for the density and distribution function of the corresponding random variable, allowing a straightforward maximum likelihood estimation. Instead specific procedures have to be used. For the simplest two-correlations linear model two matrices Rð0Þ and Rð1Þ have to be estimated. This can be done in a pairwise manner as is done in the case of the multivariate normal distribution. The random variable inherits the pairwise definition of the dependence. For each pair of variables (i, j) the two correlation coefficients qi;jð0Þ and qi;jð1Þ have to be estimated. For this purpose two parameters describing properties of the dependence are needed. One possibility is to pick two classical dependence measures from Pearson’s correlation of the normal score transformed data, Spearman’s rank correlation or Kendall’s s. However these parameters all focus on the strength of dependence and not on value-related differences in dependencies. As an alternative, one may use the measure of dependence asymmetry as introduced for spatial statis- tics in Bárdossy (2008) ai;j ¼ X n m¼0 ðFiðxm;iÞ � 1 2 Þ2ðFjðxm;jÞ � � 1 2 Þ þ ðFiðxm;iÞ � 1 2 ÞðFjðxm;jÞ � 1 2 Þ2 � ð14Þ where n is the number of samples and Fi and Fj are the empirical distribution functions of variables i and j. This measure describes the difference between the dependence of high values and that of the low values. It is well suited for variables with positive dependence, but for negative dependence it is better to first invert one of the variables. The bi-variate two correlation linear model is the sim- plest model. Its parameters can be estimated using the Pearson correlation and the dependence asymmetry. The estimation can be done by using large simulated samples. The algorithm is as follows: 1. The observed data are transformed to normal using the normal score transformation. 2. The Pearson correlation q and the asymmetry a of the transformed data are calculated. 3. A pair of correlations ðrð0Þ; rð1ÞÞ is selected 4. A random sample of the size N of the two correlations model is simulated with the correlations ðrð0Þ; rð1ÞÞ denoted as SZ ¼ fðz�1ðnÞ; z�2ðnÞÞ; n ¼ 1; . . .;Ng. 5. The correlation and the asymmetry q� and a� of the simulated data are calculated. 6. The difference function: Dðq�; a�Þ ¼ ðq� q�Þ2 þ ða� a�Þ2 is then minimized using an appropriate algorithm (for example steepest descent method). In the minimization at each new evaluation an new random sample of the size N is generated. The optimal ðrð0Þ; rð1ÞÞ is taken as parameters of the model. As the simulation of the two correlations model is very simple and fast very large sample sizes N can be generated to assure that the parameters become stable. The Chebysev Fig. 4 Simulated empirical copulas corresponding to the three-correlations model with correlation matrices defined in (12) and (13) 2624 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 123 inequality can be used to estimate the appropriate sample size. As the above procedure is computationally intensive a numerical approximation of the function Gðrð0Þ; rð1ÞÞ ! ðq; aÞ can be done using a regular grid of ðrð0Þ; rð1ÞÞ values via simulation. The obtained table can then be used to estimate the two parameters of the two correlations model. Figure 5 shows the Pearson correlation and the asym- metry as a function of the lower (rð0Þ) and the upper (rð1Þ) correlation for the two-correlation linear model. The two measures are kind of orthogonal allowing a simple esti- mation of qð0Þ and qð1Þ from the observed Pearson cor- relation and asymmetry. For the multivariate case the parameters of the two correlations model are estimated in a pairwise manner as in the case of the multivariate normal distribution. An interesting case is where correlations change their sign - for example low values show a negative correlation and high values are positively correlated. The overall Pearson correlation of such variables is usually close to zero, but the dependence is present and can be detected by other measures such as entropy. Note that the number of parameters for the simplest two correlations linear model is only twice as much as for the multivariate normal distribution. For the estimation of the parameters of a more complex dependencies additional measures have to be used. These could correspond to higher moments. As an alternative one may also use a set of indicator correlations. These are defined for different thresholds 0\h\1: IhðZkÞ ¼ 1 if FkðZkÞ[ h 0 else � ð15Þ The indicator correlations qIðhÞi;j ¼ Corr ðIhðZiÞ; IhðZjÞÞ ð16Þ considered as a function of h can also be used to to see if a dependence is symmetrical or not. For copulas like the Gaussian or the t these functions are symmetrical around 0.5, the indicator correlations qIðhÞ and qIðð1� hÞ should be equal. The parameter estimation of the more complex models is also based on large simulated random samples. In this case the difference function of the form: DðwÞ ¼ ðq� q�Þ2 þ ða� a�Þ2 þ Z 1 0 ðqIðhÞ � q�I ðhÞÞ 2 dh should be minimized. For the mþ 1 correlation model for m[ 1 the break- points si also have to be estimated. This changes the parameter estimation problem - as a pairwise estimation is not possible due to the common breakpoints. In order to simplify this problem it is reasonable to select for a given m the sj ¼ j m for j ¼ 0; . . .;m. Note that as the third example in Sect. 5 defines random variables symmetrical dependence with respect to the values, thus one cannot use the asymmetry measure defined in (14) for the estimation of the parameters. Instead one could use indicator correlations for the parameter estimation. Fig. 5 Pearson correlation (left) and asymmetry (right) as a function of the lower q(0) and the upper q(1) correlation for the two-correlations model Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 2625 123 7 Application As an example groundwater quality data from the state of Baden-Württemberg in Germany are considered. In the framework of regular groundwater observation samples of more than 2300 wells were collected. These were analysed for their chemical composition. The parameters selected are chloride, nitrate, pH, sulfate and oxygen. The basic statistics of the data are listed in Table 1. The data are highly skewed, with the exception of pH, which has a small negative skew. Due to the high skew the Pearson correlation of the data is strongly influenced by the large values of the data, and thus it is an imperfect description of the strength of the pairwise dependence. In order to investigate the interde- pendence between the different parameters the data were transformed to normal using the normal score transforma- tion. As a next step the Pearson correlations (Table 2) and the dependence asymmetry according to equation (14) were calculated for each pair of parameters. In order to investigate the Gaussianity of the dependence, the asym- metry of the Gaussian distribution was calculated for simulated samples with the same size and Pearson corre- lations N ¼ 1000 times. The corresponding asymmetries were calculated and the 95 % confidence interval was identified. For all 10 pairs of parameters the asymmetry was outside the confidence interval of the Gaussian, with the same correlation, indicating that the dependence is not normal. A two-correlations linear model was fitted to the data using the Pearson correlations and the asymmetries as parameters. The two correlation matrices Rð0Þ and Rð1Þ were estimated in a pairwise manner (Table 3). It is interesting to observe that the Pearson correlation between pH and Sulfate is very low, while the two-corre- lation model shows a relationship with changing sign. The lower correlation (corresponding to s ¼ 0) is positive 0.70 while the upper corresponding to s ¼ 1 is negative �0.52. This changing relationship cannot be captured by the normal copula (and also not by any other commonly used copulas). The indicator correlations shown on Fig. 6 con- firms the changing relationship between the two variables, which is reasonably well captured by the two-correlations model, and are not captured by the normal model. This leads to a loss of information when applying the normal copula based model. A three correlation model was also fitted to the data, such that only the relationship between pH and Sulfate was altered. The three correlations were assessed using the indicator correlations. The new rela- tionship is now 0.85 ! 0.00 ! �0.3 (corresponding to si =0, 0.5 and 1). Note that the correlation for the low values corresponding to s ¼ 0 increased, but decreases faster. Figure 6 shows the improvement of the fit to the indicator correlations. A comparison of the two dimensional marginals of the simulated data with the observed ones shows that in all cases the Kolmogorov distance between the simulated two- correlations model and the observations, is lower than for the normal copula based simulations. For the second example, air quality measurements taken near Zurich-Schimmelstrasse in Switzerland were used. Four parameters NOx, SO2, NO2 and PM10 were selected. These data are publicly available for the time period of 2011–2021 on the internet under opendata.swiss. Basic statistics of the data are listed in Table 4. The Pearson correlations of the normal score transformed data are given in Table 5. The correlations are higher than those of the groundwater example. The matrices of the two-correlations model were fitted to the data by using the Pearson correlation and the asymmetry as measures using the algorithm described in Sect. 6. The fit was done numerically by using large sim- ulated samples for the estimation of the parameters. The pairwise calculated asymmetries showed that out of the 6 pairs, 4 are significantly different from the normal. The theoretical model is in this case contains two pairs for which ri;jð0Þ ¼ ri;jð1Þ. In order to investigate the appro- priateness of the model, N ¼ 1000 simulations of the 4 dimensional distributions were considered using the two- correlation linear model with parameters listed in Table 6. Table 1 Basic statistics of the observed groundwater quality data Chloride Nitrate pH Oxygen Sulfate mg/l mg/l [-] mg/l mg/l Mean 32.37 31.97 7.34 8.03 67.07 Standard deviation 41.10 27.98 0.41 3.00 120.68 Skewness 7.11 2.30 2.28 - 0.53 6.93 Minimum 0.60 0.20 5.10 0.20 0.70 Maximum 759.00 265.00 13.10 16.70 1628.00 Sample size 2537 2537 2537 2537 2537 Table 2 Pearson correlations of the normal score transformed observed groundwater quality data Chloride Nitrate pH Oxygen Sulfate Chloride 1.00 0.45 0.04 0.35 0.63 Nitrate 0.45 1.00 0.14 � 0.12 0.37 pH 0.04 0.14 1.00 � 0.16 0.12 Oxygen 0.35 � 0.12 � 0.16 1.00 0.35 Sulfate 0.63 0.37 0.12 0.35 1.00 2626 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 123 The same number of simulations were carried out for the normal copula case. As an example Fig. 7 shows the empirical copula a simulated normal and a simulated changing correlation copula. The sample sizes are in all cases the same. As one can see, the strong dependence of the low values is not captured by the normal simulation, but well captured by the two-correlations model. A com- parison of the two dimensional marginals shows that the two-correlation linear model captures the asymmetrical dependence well. Note that in contrast to the previous example in this case all correlations are positive. Deviations from the Gaussian dependence are very fre- quent, and both examples include cases where the high values have stronger dependence than the low ones and the reverse case too. 8 Discussion and conclusions In this paper a method to construct multivariate non- Gaussian distributions was presented. The construction is very general and special cases with 2 or more parameters for the description of dependence of the pairs can be used. The copulas obtained via these distributions are not only useful for monotonic dependence but can also represent dependencies with changing character, for example nega- tive dependence for small values and positive dependence for high values. Copulas were most frequently used for the investigation and description of the dependence of extremes. However, the dependence of the variables might be applied to non- extreme values and also deviations from the normal. The copulas defined in this paper can describe arbitrary upper and lower tail dependence, and can also be used for the description of asymmetrical dependence even without focusing on the extremes. A disadvantage of this construction is that the distribu- tion functions do not have a general closed analytical form. This makes the estimation of the parameters difficult. The distributions defined using this construction may have many parameters, but due to the increasing data volumes of Fig. 6 Indicator correlations for pH and Sulfate - blue = observed, orange = two-correlations model, red = three-correlations model and green dashed = normal model Table 4 Basic statistics of the observed air quality data NOx SO2 NO2 PM10 measurement unit ppb lg/m3 lg/m3 lg/m3 Mean 39.65 18.61 43.92 28.93 Standard deviation 15.12 10.84 27.44 26.18 Skewness 0.56 1.77 1.99 2.64 Minimum 6.17 1.94 4.87 1.77 Maximum 109.44 109.14 272.87 271.14 Sample size 3903 3903 3903 3903 Table 5 Pearson correlations of the normal score transformed air quality data NOx SO2 NO2 PM10 NOx 1.00 0.70 0.93 0.85 SO2 0.70 1.00 0.66 0.59 NO2 0.93 0.66 1.00 0.98 PM10 0.85 0.59 0.98 1.00 Table 3 Changing correlations of the normal score transformed observed groundwater quality data with the two-correlations model. (qijð0Þ ! qijð1Þ) Chloride Nitrate pH Oxygen Sulfate Chloride 1.00 0.68 ! 0.22 0.58 ! -0.54 0.38 ! 0.36 0.84 ! 0.40 Nitrate 0.68 ! 0.22 1.00 0.52 ! -0.26 0.14 ! -0.38 0.56 ! 0.16 pH 0.58 ! -0.54 0.52 ! -0.26 1.00 0.08 ! -0.42 0.70 ! -0.52 Oxygen 0.38 ! 0.36 0.14 ! -0.38 0.08 ! -0.42 1.00 0.38 ! 0.36 Sulfate 0.84 ! 0.40 0.56 ! 0.16 0.70 ! -0.52 0.38 ! 0.36 1.00 Pairs with non-normal dependence at 95 % significancance level are in boldface Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 2627 123 environmental variables may be well used for the investi- gation of big data. The development of parameter estimation methods for more complex models requires further research. The presented construction yields bi-variate marginal copulas which are symmetrical with respect the main axis addedof the copula. Non-symmetrical dependence can be achieved by using the same idea with other multivariate distributions, such as the skew-normal distribution. In this case both the parameters of the correlation matrix as the k parameter defining the skewness can be varied in the same way as done in equation (5). The use of these copulas as an alternative for multi- variate linear regression is also possible, but goes beyond the scope of this contribution. The methodology can be extended to time series and spatial random fields. Appendix Proposition For any 0� a� 1 there is a bi-variate construction such that the upper tail dependence is exactly a Proof Let U2ðz1; z2; qÞ be the bi-variate normal distribution function with standard normal marginals and q correlation. For any 0� a\1 and s there is a correlation q� such that the exceedence probability 1� 2sþ U2ðU�1ðsÞ;U�1ðsÞ; q�Þ 1� s ¼ a As for q ! 1[ a the left hand side converges to 1 and for q ¼ 0 the left hand side is 0 by the continuity there must be a q� which fulfills the equation. Defining qðsÞ ¼ q� leads to the construction. For a ¼ 1 a function converging to 1 has to be selected at the right hand side. The above construction can also be used to obtain variables with a lower tail dependence a. This way a random variable with a lower tail dependence a and an upper tail dependence b 6¼ a can be constructed. Acknowledgements Research leading to this paper has partly been supported by the German Research Foundation (DFG), by funding of the research group FOR 2416‘‘Space-Time Dynamics of Extreme Floods (SPATE)’’.The author thanks the anonymous reviewer and G. Pegram for their useful comments, and dedicates this paper to his memory. Author Contributions AB did everything. Funding Open Access funding enabled and organized by Projekt DEAL. Declarations Conflict of interest The authors declare no competing interests 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 Fig. 7 Empirical and simulated copulas for the air quality parameters NOx and NO2 for Zürich: left simulated using normal bivariate copula, middle observed, right simulated using the two-correlations model Table 6 Changing correlations of the normal score transformed air quality data with the two- correlations model (qijð0Þ ! qijð1Þ) NOx SO2 NO2 PM10 NOx 1.00 0.76 ! 0.63 0.99 ! 0.86 0.95 ! 0.77 SO2 0.76 ! 0.63 1.00 0.69 ! 0.63 0.54 ! 0.65 NO2 0.99 ! 0.86 0.69 ! 0.63 1.00 0.99 ! 0.98 PM10 0.95 ! 0.77 0.54 ! 0.65 0.99 ! 0.98 1.00 Pairs with non-normal dependence at 95 % significancance level are in boldface 2628 Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 123 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/. References Bárdossy A (2006) Copula-based geostatistical models for ground- water quality parameters. Water Resour Res. https://doi.org/10. 1029/2005WR004754 Bárdossy A, Li J (2008) Geostatistical interpolation using copulas. Water Resour Res. https://doi.org/10.1029/2007WR006115 Bardossy A, Pegram G (2012) Multiscale spatial recorrelation of rcm precipitation to produce unbiased climate change scenarios over large areas and small. Water Resour Res. https://doi.org/10.1029/ 2011WR011524 Brunner MI, Furrer R, Favre A-C (2019) Modeling the spatial dependence of floods using the fisher copula. Hydrol Earth Syst Sci 23(1):107–124. https://doi.org/10.5194/hess-23-107-2019 Chen L, Guo S (2019) Copula-based flood frequency analysis. Copulas Appl Hydrol Water Resour. https://doi.org/10.1007/ 978-981-13-0574-0_3 Czado C, Nagler T (2022) Vine copula based modeling. Annu Rev Stat Appl 9:453–477. https://doi.org/10.1146/annurev-statistics- 040220-101153 Favre A-C, Quessy J-F, Toupin M-H (2018) The new family of fisher copulas to model upper tail dependence and radial asymmetry: Properties and application to high-dimensional rainfall data. Environmetrics 29(3):e2494 Gräler B, van den Berg MJ, Vandenberghe S, Petroselli A, Grimaldi S, De Baets B, Verhoest NEC (2013) Multivariate return periods in hydrology: a critical and practical review focusing on synthetic design hydrograph estimation. Hydrol Earth Syst Sci 17(4):1281–1296. https://doi.org/10.5194/hess-17-1281-2013 Guthke P, Bárdossy A (2012) Reducing the number of MC runs with antithetic and common random fields. Adv Water Resour 43:1–13. https://doi.org/10.1016/j.advwatres.2012.03.014 Joe H (1997) Multivariate models and dependence concepts. Chap- man Hall, Boca Raton Nelsen RB (1999) An introduction to copulas. Springer, New York Sklar A (1959) Fonctions de répartition à n dimensions et leurs marges. Publ Inst Stat Paris 8:229–231 Won J, Choi J, Lee O, Kim S (2020) Copula-based joint drought index using SPI and EDDI and its application to climate change. Sci Total Environ 744:140701 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Stochastic Environmental Research and Risk Assessment (2023) 37:2619–2629 2629 123 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ https://doi.org/10.1029/2005WR004754 https://doi.org/10.1029/2005WR004754 https://doi.org/10.1029/2007WR006115 https://doi.org/10.1029/2011WR011524 https://doi.org/10.1029/2011WR011524 https://doi.org/10.5194/hess-23-107-2019 https://doi.org/10.1007/978-981-13-0574-0_3 https://doi.org/10.1007/978-981-13-0574-0_3 https://doi.org/10.1146/annurev-statistics-040220-101153 https://doi.org/10.1146/annurev-statistics-040220-101153 https://doi.org/10.5194/hess-17-1281-2013 https://doi.org/10.1016/j.advwatres.2012.03.014 Changing correlations: a flexible definition of non-Gaussian multivariate dependence Abstract Introduction Methodology Copulas Distributions with value dependent correlations Construction of continuous \Sigma \lpar \tau \rpar functions Construction using the square roots of the matrices Construction using spatial statistics Simulation Examples Parameter estimation Application Discussion and conclusions Appendix Author Contributions Funding References