Causal modelling of heavy-tailed variables and confounders with application to river flow

Confounding variables are a recurrent challenge for causal discovery and inference. In many situations, complex causal mechanisms only manifest themselves in extreme events, or take simpler forms in the extremes. Stimulated by data on extreme river flows and precipitation, we introduce a new causal discovery methodology for heavy-tailed variables that allows the effect of a known potential confounder to be almost entirely removed when the variables have comparable tails, and also decreases it sufficiently to enable correct causal inference when the confounder has a heavier tail. We also introduce a new parametric estimator for the existing causal tail coefficient and a permutation test. Simulations show that the methods work well and the ideas are applied to the motivating dataset. Supplementary Information The online version contains supplementary material available at 10.1007/s10687-022-00456-4.


Introduction
The field of causal inference has developed massively in recent decades (e.g., Pearl, 2009;Peters et al., 2017), with much recent work on the detection of causality from observational data (e.g., Maathuis and Nandy, 2016).Most of this literature concerns central quantities such as expectations, but certain causal mechanisms manifest themselves only in rare events and/or may simplify in distribution tails.Standard methods of causal inference are ill-suited for such situations, and recent work has begun to link causality and extreme value theory.Examples are Gissibl and Klüppelberg (2018), who define recursive max-linear models on directed acyclic graphs, Klüppelberg and Krali (2021), who propose a scaling technique to determine the causal order of the variables in such graphs, Kiriliouk and Naveau (2020), who use multivariate generalized Pareto distributions to study probabilities of necessary and sufficient causation as defined in the counterfactual theory of Pearl, and Mhalla et al. (2020), who construct a causal inference method for tail quantities relying on Kolmogorov complexity of extreme conditional quantiles.See surveys by Naveau et al. (2020) on extreme event attribution and by Engelke and Ivanovs (2021) on the detection and modeling of sparse patterns in extremes.
Our work stems from that of Gnecco et al. (2021), who propose an estimator of the causal tail coefficient and an algorithm that, under mild conditions, consistently retrieves a causal order on an underlying graph even in the presence of hidden confounders.Such an order helps to exclude some causal structures, but does not provide evidence for the existence of a specific structure, as in general a given order is causal for several possible graphs; in particular, all orders are causal for the empty graph corresponding to absence of causality.Although it is asymptotically invariant to hidden confounders, this estimator can suffer from confounding in finite samples when inference on the direct relationship between two variables is needed, when these effects are too strong or when the confounders have heavier tails than the two variables.
This paper addresses a central challenge in causal inference: the presence of confounders.In theoretical development it is often assumed that all the relevant variables are observed and can be included in the model, but in practice one can rarely be sure of this.The available variables are often subject to external influences, observed or unobserved, that affect the variables of interest and can make it harder or even impossible to infer a correct causal relationship.Our goals are to mitigate the effect of a set of known confounders on an extremal causal analysis by treating them as covariates, and to present a permutation test for direct causality between the two observed variables.Our approach relaxes the assumption of Gnecco et al. (2021) that the confounders have the same tail index as the two main variables of interest, and thus encompasses a much broader range of situations, such as that in our application.Such a model enables causal discovery and inference for a greater variety of situations.
Our work was stimulated by average daily discharge data from 68 gauging stations along the Rhine and Aare catchments in Switzerland, see Figure 1.The data were collected by the Swiss Federal Office for the Environment (hydrodaten.admin.ch),but were provided by the authors of Engelke and Ivanovs (2021), with some useful preliminary insights.We focus on the causal relationship between extreme discharges, for which precipitation is an obvious confounder, and use daily precipitation data from 105 meteorological stations, provided by the Swiss Federal Office of Meteorology and Climatology, MeteoSwiss (gate.meteoswiss.ch/idaweb).Unlike in our simulation experiments, we know neither the true tail properties of the discharges and precipitation nor the effect of the confounder.We use precipitation as a covariate in our test, allowing inference on the direct causal relationships between discharges for the majority of the station pairs, with at least 95% estimated confidence, which was impossible without our proposed approach.
The paper is organised as follows.Section 2 discusses the causal tail coefficient, its interpretation and its properties.Section 3 introduces a new parametric estimator for it based on generalized Pareto modelling of threshold excesses, which allows a known confounder to be used as a covariate.A simulation study in Section 4 underlines the strengths and limitations of the two estimators.Section 5 presents a permutation test intended to detect direct causality between two heavy-tailed variables, which is also assessed via simulation.Section 6 applies the methodology to the river discharges, and Section 7 gives a brief discussion. 2 Causal Tail Coefficient and its Estimation

Existing Work
We first give some basic notions needed to describe the setting in which causal relationships between random variables can be recovered.Definition 2.1.A linear structural causal model (LSCM) over a set of random variables X 1 , . . ., X p satisfies X j = k∈pa(j) where V := {1, . . ., p} is a set of nodes representing the corresponding random variables, pa(j) ⊆ V is the set of parents of j, β jk ∈ R \ {0} is called the causal weight of node k on node j, and ε 1 , . . ., ε p are jointly independent noise variables.We suppose that the associated graph G = (V, E), in which the directed edge (i, j) ∈ V × V belongs to E if and only if i ∈ pa(j), is a directed acyclic graph (DAG).
In a DAG G = (V, E), we say that i ∈ V is an ancestor of j ∈ V in G, if there exists a directed path from i to j.The set of the ancestors of j in G is denoted by An(j, G), and we define an(j, G) := An(j, G) \ {j}.In a LSCM over random variables X 1 , . . ., X p , with associated DAG G = (V, E), we say that X i causes X j , if i ∈ an(j, G).We call X i a confounder (or common cause) of X j and X k if there exist directed paths from i to j and from i to k in G that do not include k and j, respectively.We say that there is no causal link between X i and X j if An(i, G) ∩ An(j, G) = ∅.For any i, j ∈ V we let β i→j denote the sum of the products of the causal weights along the distinct directed paths from vertex i to vertex j; we set β j→j := 1 and Let X i and X j be random variables from a LSCM with respective distributions F i and F j .The causal (upper) tail coefficient of a random variable X i on another random variable X j is defined as (Gnecco et al., 2021) if the limit exists.This coefficient lies between zero and one and captures the causal influence of X i on X j in their upper tails: if X i has a linear causal effect on X j , Γ 1,2 will be close to unity.The coefficient is asymmetric, as extremes of X j need not lead to extremes of X i , and in that case, Γ ji will be appreciably smaller than Γ ij .As Γ ij only depends on the rescaled margins of the variables, it is invariant to monotone increasing marginal transformations.
If both tails are of interest, the causal tail coefficient can be generalized to capture the causal effects in both directions, by considering the symmetric causal tail coefficient of X i on X j , i.e., if the limit exists, where ρ : The interpretation and properties of Ψ ij are similar to those of Γ ij .The symmetric version captures the causal influence of X i on X j in both of their tails.
For simplicity we focus on Γ ij in this paper, though all of our results and methods can be generalized to both tails by considering Ψ ij instead, if the assumptions for the upper tails are also satisfied in the lower tails of the variables considered.
Before stating the theorem that describes how the underlying causal relationships in a set of random variables can be recovered, we define the concept of regular variation.Definition 2.2.A positive measurable function f is said to be regularly varying with index α ∈ R, Definition 2.3.The random variable X j is said to be regularly varying with index α > 0, if, for some ∈ RV 0 , P(X j > x) ∼ (x)x −α as x → ∞.
The following theorem describes how the causal relationships underlying a set of random variables can be recovered from their causal tail coefficients.
Theorem 2.4 (Gnecco et al., 2021).Let X 1 , . . ., X p be random variables from a LSCM, with associated directed acyclic graph G = (V, E) and suppose that (a) the coefficients β jk of the linear structural causal relationship X j = k∈pa(j,G) β jk X k + ε j are strictly positive for all j ∈ V and k ∈ pa(j, G), and (b) the real-valued noise variables ε 1 , . . ., ε p are independent and regularly varying with comparable upper tails.
Then the values of Γ ij and Γ ji allow one to distinguish between the different possible causal relationships between X i and X j summarized in Table 1.
Table 1: Equivalence of the possible values of Γ ij and Γ ji with the underlying causal relationship between X i and X j .
Under the theorem's assumptions, the blank entries in Table 1 cannot occur.Theorem 2.4 is generalizable to the Ψ ij variant of the coefficient and possibly negative β ij values if the assumptions are also satisfied in the lower tails of the variables.Gnecco et al. (2021) show that under the setup and assumptions of Theorem 2.4, the causal tail coefficient (1) for any distinct i, j ∈ V , and with A ij := An(i, G) ∩ An(j, G), is Without loss of generality we set i = 1 and j = 2 in what follows, and thus consider the causal effect of X 1 on X 2 .
If {(X i,1 , X i,2 )} n i=1 are independent replicates of (X 1 , X 2 ), with the random variables X i and X j from the LSCM, then the non-parametric estimator of Γ 1,2 is defined to be for some k ∈ {1, . . ., n − 1}, where 1(•) denotes the indicator function, X (h),1 denotes the h th order statistic and Fj is the empirical cumulative distribution function of X j , i.e., This estimator is the empirical counterpart to (1), as X (h),1 = F ← 1 (h/n) is a quantile of the corresponding empirical distribution.The value of k controls the number of data pairs in the upper tail of X 1 that contribute to the estimator.Under the assumptions of Theorem 2.4 and a "very mild assumption that is satisfied by most univariate regularly varying distributions of interest", estimator (3) is consistent as n → ∞, for a choice of k such that k → ∞ and k/n → 0 (Gnecco et al., 2021).

Practical Limitations
A strength of the causal tail coefficient approach is its asymptotic robustness to hidden confounders.Studies of causation frequently presuppose that all the relevant variables have been observed, which is usually moot, but Theorem 2.4 holds even when some variables in the underlying LSCM are unobserved.This capacity to deal with confounders both when studying the causal relationship between two variables and when retrieving a causal order is not generally shared by other approaches in causal inference, as argued by Gnecco et al. (2021, Section 4.2), but the unobserved variables must satisfy a regular variation assumption that is hard to check and may be unrealistic.In practice, moreover, the tail behaviour of the confounders may differ from that of X 1 and X 2 , violating assumption (b) of Theorem 2.4.In our motivating setting, for example, the tail of the confounder, precipitation, may not behave like the tails of the river discharges.This problem worsens when the confounder has a heavier tail than the variable of interest.Furthermore, distinguishing between different causal situations using empirical estimates may be difficult; an increase in the strength of the causal effect of a common confounder of X 1 and X 2 will increase Γ 1,2 , making it harder to tell whether a high value of Γ1,2 indicates that Γ 1,2 = 1 or that Γ 1,2 1, as we shall see in Section 4.
The discussion above suggests that conditioning on the values of known confounders might be valuable.In the presence of a vector H of potential confounders we therefore define unless there are no confounders.If X 1 causes X 2 , on the other hand, then Γ 1,2|H = Γ 1,2 = 1.In the presence of potential confounders, therefore, (4) seems preferable to Γ 1,2 .The difficulty is that the estimation of (4) requires the modelling of the dependence of both X 1 and X 2 on H.The first is more straightforward, because for large u only the upper tail of X 1 need be considered, whereas the second ostensibly requires a model for the entire distribution of X 2 , and this may be complex.We compromise by fitting similar models to both variables, letting the upper tails alone vary with H.As we shall see below, this can greatly improve estimation of the causal dependence structure relative to the original approach.Moreover fitting such a model should highlight simpler, potentially linear, structures in the tails, rather than more complex ones in the body of the data.This leads us to propose a peaks-over-threshold approach to estimating the conditional dependence of X 1 and X 2 on H (Section 3).Another useful tool, a reliable statistical test for direct causality, is discussed in Section 5.
3 Parametric Tail Causality and Confounder Dependence

Generalized Pareto Causal Tail Coefficient
As mentioned above, we use the generalized Pareto distribution (GPD) to model the tails of our variables (Coles, 2001, Chapter 4).For j = 1, 2, and under mild conditions on X j , for a large enough threshold u j large enough, we have with a scale parameter σ j > 0 and a shape parameter ξ j ∈ R: • ξ j = 0 corresponds to light-tailed distributions, and then X j lies in the maximum domain of attraction of the Gumbel distribution; • ξ j > 0 corresponds to heavy-tailed distributions, and then X j lies in the maximum domain of attraction of the Fréchet distribution; and • ξ j < 0 corresponds to distributions with bounded upper tails, and then X j lies in the maximum domain of attraction of the (reverse) Weibull distribution.
Any random variable satisfying the assumptions of Theorem 2.4 satisfies ( 5), as a regularly varying random variable with index α > 0 lies in the Fréchet maximum domain of attraction.If the threshold u j is chosen to be the q quantile of X j for some q ∈ (0, 1), then we can write and using the empirical distribution F (x) to estimate P(X j ≤ x) and maximum likelihood estimation using the excesses of u j to obtain σj and ξj yields a hybrid estimator of the distribution function F j (x) of X j , i.e., The choice of q involves a bias-variance trade-off: q should be chosen large enough for the tail to be well approximated by a GPD, thus reducing the bias, but small enough to have enough exceedances, thus reducing the variance of the estimator.Using hybrid estimators for F 1 and F 2 for an integer k ∈ {1, . . ., n − 1} yields the parametric GPD causal tail coefficient estimator for where k g := |{i ∈ {1, . . ., n} : F1 (X i,1 ; σ1 , ξ1 ) > 1 − k/n}|.Unlike with the non-parametric estimator (3), the number of data pairs k g used in (6) may not equal k, as it depends on the fit of F1 (X i,1 ; σ1 , ξ1 ).
The GPD model can be extended to allow dependence on covariates of interest by expressing its parameters in the form θ(i) = h{γ Z(i)}, where θ denotes one or both of σ and ξ, h is an inverse link function, γ is a vector of parameters and Z(i) is the vector of explanatory variables on which the model might depend (Davison and Smith, 1990).
We wish to reparametrise the model to reduce or remove the effect on Γ 1,2 of a vector of potential confounders H of X 1 and X 2 .If H is part of the LSCM then under the setup in Section 2 it is straightforward to show that H affects the scale parameters of the GPD model that applies to X 1 and X 2 above high thresholds, but not their shapes, so we write where H i is the replicate of H corresponding to the observations This yields, for k ∈ {1, . . ., n − 1}, the parametric H-conditional linear generalized Pareto distribution (LGPD) causal tail coefficient estimator, where k l := |{i ∈ {1, . . ., n} : F1 {X i,1 ; σ1 (i), ξ1 } > 1 − k/n}|.Estimation of σ 0 j , σ 1 j and ξ j is performed by maximum likelihood.In applications it is preferable to center and rescale each confounder in H componentwise to unit variance and zero mean, to avoid numerical issues.Although the confounder is here assumed to be part of the LSCM, this does not seem to be necessary in practice, as non-linear effects can be approximated linearly, especially in the tail region.We investigate the effect of varying the tail index in Section 4.2.

The Positive Linear Scale Issue
Linear modelling of the GPD scale parameter may not yield positive scale estimates σj (i) > 0 for each i = 1, . . ., n and j = 1, 2. The use of a nonlinear link function to ensure that the scale estimates were positive would not agree with the assumption of extremal linearity of the causal relationships, as the effect of H on the scale is also necessarily linear.We now describe two different solutions to this problem, which we compare by simulation in Section 4.
The first solution, post-fit correction, replaces σj (i) in (8) by max{σ j (i), } for some arbitrary but small positive .The second solution, the constrained approach, applies the following linear constraints to the estimates when maximizing the likelihood where min i=1,...,n H i and max i=1,...,n H i represent the vectors of componentwise minima and maxima.When the data have a known distribution, box constraints can be used instead of (9).For example, in the case of a single confounder H and if X 1 , X 2 and where the lower and upper bounds are needed for positive and negative H i , respectively.

Simulation Study
Here we perform a simulation study using the Student t, Pareto and log-normal noise distributions.
The first two lie in the Fréchet maximum domain of attraction and are regularly varying with index α = 1/ξ > 0. We write Pareto(a, α) for the Pareto model with scale parameter a and tail index α; recal that lower values of α indicate heavier tails.This distribution satisfies Definition 2.3 exactly, so one might expect Pareto data to show better behaviour than Student data.The log-normal distribution, LogN(µ, σ 2 ) lies in the maximum domain of attraction of the Gumbel distribution and is not regularly varying, but finite samples from it can appear to be heavy-tailed.
We focus on the behaviour of the causal tail coefficient estimators (3) and (8) between two variables X 1 and X 2 in their causal configurations, as shown in Figure 2. As we study the estimators of causal effects of both X 1 on X 2 and of X 2 on X 1 , we generated simulations only for the four causal cases, A, B, C and D. The LSCM causal weights β 2,1 , β 1h and β 2h were chosen to equal 1.0, by default, for each existing edge in all four cases.Hence, in D, X 2 is caused by X 1 and the single confounder H with equal strength, even though H has another effect on X 2 through X 1 .
Unless stated otherwise, each estimate is based on a random sample of n = 10 6 triples (X 1 , X 2 , H), of which k = 2 n 0.4 = 502 were chosen - Gnecco et al. (2021) found that the optimal fractional exponent of n for choosing k seems to lie between 0.3 and 0.4.The factor 2 doubles the number of data pairs used in the estimator, thus decreasing its variability, but does not introduce much bias for such large values of n.The GPD-based estimators are based on the top (1 − q)n observations, where we take q = 0.9, though only around k of the largest observations are used to estimate the coefficients Γ ij .Setting q = 0.95 yields similar results.One thousand independent replicates were generated for each of the four causal configurations and three distributions.
We present only the highlights of the study; the code and all the results are available from github.com/opasche/ExtremalCausalModelling.

Variables with Comparable Tails
Detailed results for variables with comparable tails may be found in Section S.1 of the Supplementary Material.In this case it is essentially always possible to infer the existence and direction of any causality between X 1 and X 2 , based on the non-parametric or H-conditional LGPD estimators, (3) or (8), of Γ 1,2 and Γ 2,1 alone.When the causal effects of H on X 1 and X 2 , i.e., β 1h and β 2h , are increased relative to the noise variance and any causal effect β 2,1 of X 1 on X 2 , both Γ 1,2 and Γ 2,1 increase in configuration B, and Γ 2,1 increases in configurations C and D. This increase is larger with the non-parametric estimators of Γ 1,2 and Γ 2,1 , which are biased upwards in these configurations.When the confounder has a high causal impact, inference based on the non-parametric estimator (3) for direct causal link between X 1 and X 2 can fail, as Γ1,2 , Γ2,1 ≈ 1 and hence | Γ1,2 − Γ2,1 | ≈ 0 in configurations B and D.
Use of the H-conditional LGPD estimator (8) greatly reduces the effect of H on the coefficient estimates in configurations B and D. For Pareto and log-normal data, the results are indistinguishable from those without the confounder, both in terms of location and variability, as if the effect of H had been entirely removed.The estimates based on Student data are also shifted to around the same values as in the corresponding confounder-free configurations, though their upper tails are marginally heavier.These few greater values remain appreciably lower than without H as a covariate.For configurations A and C, unlike for B and D, the estimator is almost unaffected by the addition of H as a covariate when it is not a confounder.This is also a useful property, as it could allow tests of whether a specific covariate is a confounder of two variables, based on changes to the estimated coefficients.

Confounder with a Different Tail
One generalisation allows the tail of the distribution of H to be heavier or lighter than those of X 1 and X 2 .A lighter tail does not negatively affect whether the non-parametric and H-conditional LGPD estimators can infer a direct causal relationship between X 1 and X 2 , as the tails of X 1 and X 2 then dominate.Figure 3 shows the sampling distributions of Γ1,2 and Γ2,1 for all four causal structures when the tail of H is heavier than those of X 1 and X 2 .The true coefficient values are unknown, as assumption (b) of Theorem 2.4 is not satisfied, though the coefficient for comparable tails, (2), is shown for comparison.
When H has a heavier tail than X 1 and X 2 , the non-parametric estimators Γ1,2 and Γ2,1 in configuration B and Γ2,1 in configuration D are shifted well towards unity.With an even heavier-tailed, Student t 2 , distribution for H (not shown here), the Student results resemble those for the Pareto and log-normal distributions.In all these cases it becomes impossible to infer a direct causal relationship between X 1 and X 2 , owing to the effect of the heavier confounder tail on the non-parametric estimators.
Figure 3 shows that in configurations B and D the non-parametric estimator is badly affected by the heavier tail of H. Figure 4, which displays the sample distributions of ΓGPD 1,2|H and ΓGPD 2,1|H with post-fit correction when the tail of H is heavier than those of X 1 and X 2 , shows that the use of H as a covariate solves this problem: the estimates shift towards the coefficient values in the corresponding confounder-free cases, and consistently yield positive values of the difference of estimates ΓGPD 1,2|H − ΓGPD 2,1|H for configuration D and differences centred at zero for configuration B; see also Section S.1 of the Supplementary Material.The estimates in configurations A and C, without the confounder causal effect, are barely changed by using H as a covariate.

Simulation results for ΓGPD
1,2|H and ΓGPD 2,1|H with the constrained fit are very similar to those for post-fit correction for the Pareto and log-normal distributions, but not for the Student distribution.Figure 5 shows the sample distribution of ΓGPD 1,2|H and ΓGPD 2,1|H with the constrained fit, for a heavier confounder tail.For the Student distribution, the confounder affects the estimator appreciably more for the constrained fit than for post-fit correction, compared to the non-parametric results.As the Student distribution is heavy in both tails, the lower constraint in (9) forces σj (i) (j = 1, 2) to have an appreciably smaller slope, explaining this reduced effect.In configurations with a confounder, the absolute values of the constrained σ1 j may be up ten times smaller than those for post-fit correction.With both approaches σ1 j rarely differs greatly from zero for configurations without a confounder.Both types of constraint yield very similar estimates for the Student distribution; see github.com/opasche/ExtremalCausalModelling.
To summarize, the simulations show that both the non-parametric estimator (3) and the Hconditional LGPD estimator (8) perform well when the theoretical assumptions are met and the influence of a hidden confounder is limited.When this influence grows, it becomes increasingly difficult to confidently infer the causal relationship between the variables using the non-parametric estimator, but the H-conditional LGPD estimator allows us to detect this relationship by reducing the effect of the confounding.

Permutation Test
In situations such as the causal analysis presented in Section 6, the distributions of the Γ 1,2 and Γ 2,1 estimators must be estimated to be used for inference.One way to obtain such distributions would be bootstrap resampling, but the extremal nature of the causal tail coefficient would require an unrealistically large sample size for its bootstrap distributions to be trustworthy, as these distributions tend to be too discrete in the extremes.
We therefore propose a permutation test (Davison and Hinkley, 1997, Chapter 4) for direct causality between two observed variables, measuring the asymmetry in their direct causal relationship.Suppose we have a sample {(X i,1 , X i,2 )} n i=1 from a LSCM and wish to test the null hypothesis of no direct causal relationship between X 1 and X 2 , H 0 : β 2,1 = 0, versus the alternative that X 1 causes X 2 , H A : β 2,1 > 0. Our proposed procedure is as follows:  1. Rescale values Xi,j = Fj (X i,j ) (i = 1, . . ., n, j = 1, 2), where known confounders can be used in the distribution estimator Fj , as for ΓGPD 1,2|H .

Compute
4. Obtain the Monte Carlo p-value, by comparing the value of the test statistic on the original rescaled data with the permutation distribution, If there are no asymmetric confounding effects on the two variables, i.e. 2) and Theorem 2.4.This does not hold generally with asymmetric confounding.The direct causal relationship is symmetric under H 0 , i.e., X 2 is as likely to take extreme values when X 1 is extreme as is X 1 when X 2 is extreme.If so, then permutations such as those performed in step 2. are equally likely, so ∆1,2 , ∆ * 1 1,2 , . . ., ∆ * R 1,2 have a common distribution centered around zero, and p mc will be uniformly distributed.Under the alternative, the direct causal relationship is "asymmetric", as X 2 is more likely to be extreme when X 1 is extreme than conversely; then ∆1,2 is more likely to lie in the upper tail of ∆ * 1 1,2 , . . ., ∆ * R 1,2 .Thus the distribution of p mc will become increasingly skewed towards zero as the causal strength of X 1 on X 2 increases.If all asymmetric confounding effects are captured in Fj by estimating the distribution conditionally, X 1 and X 2 have comparable tails and causal effects behave linearly in the extremes, then the proposed procedure should provide a reliable p-value for testing direct causality of X 1 on X 2 .

Simulations
We used simulation from different data distributions and for different causal configurations involving X 1 , X 2 and a potential confounder H to assess our proposed test.We used values of 0, 0.01, 0.05, 0.1, 0.2 for the causal strength β 2,1 of X 1 on X 2 , with confounding effects both present and absent.Symmetric (β 1H = β 2H = 1) and asymmetric (β 1H = 0.8 and β 2H = 1, or β 1H = 1 and β 2H = 0.8) confounding effects were considered, and the noise variable were Pareto, Student t and log-normal.We generated m = 10 3 replicate samples of n = 10 4 independent triples (X i,1 , X i,2 , H i ) for each causal configuration and noise distribution.The sample size n was chosen closer to practical orders of magnitude, compared to our large-sample study in Section 4. Three versions of the permutation test were performed for each sample, corresponding to the causal tail coefficient estimators discussed in Sections 2 and 3: the non-parametric (3), and H-conditional LGPD ( 8) with either post-fit correction or constrained fit.Each used R = 10 3 permutations and the estimator hyper-parameters were set to k = 2 n 0.4 = 78 and q = 0.9.
Figure 6 shows uniform QQ-plots of p mc for the Pareto and Student distributions, in the case of heavier confounder tail, with symmetric effects.In the absence of confounding the test behaves as expected in both cases, and adding dependence on the independent H variable in the modelling through the parametric estimators has no visible effect on the distribution of p mc compared to the non-parametric approach.For the Pareto distribution, the test has a power of almost 0.9 for a direct causal strength of 0.01, and it behaves perfectly for higher causal strengths.For the Student distribution, the test reaches a power of 0.3 for a direct causal strength of 0.05, of 0.7 for causal strength of 0.1 and of near 1.0 for a causal strength of 0.2.
When the confounding effects are added, the test based on the non-parametric estimator fails for the Pareto distribution, as most of the p mc then lie outside the 95% confidence bands, indicating that the distribution of p mc is highly non-uniform.This is corrected when the value of the confounder is taken into account using the parametric approaches, with power 0.9 for a direct causal strength of only one twentieth of the confounder's marginal effects.In the Student case, p mc seems to be close to uniformity in the absence of direct causality (the difference in tail shape is much greater in the Pareto case), but post-fit correction increases the power from below 0.2 to above 0.4 for a direct causal strength of one fifth of the confounder's marginal effects.Similar conclusions to those of Section 4.2 about the constrained fit for distributions with both tails heavy apply, as the constrained fit estimator is not significantly better than the non-parametric estimator compared to post-fit correction.
Figure 7 shows the uniform QQ-plots with asymmetric confounding effects for the Pareto distribution with comparable tails.Unlike in the corresponding symmetric case, the test here fails when using the non-parametric estimator owing to the asymmetry induced by the confounder, but both parametric approaches remove this unwanted effect by enough that p mc nearly has a uniform distribution, with almost perfect power, for a causal strength of one sixteenth and one twentieth of the marginal confounding effects.

Application to Swiss Rivers
We now illustrate how our method can discover direct causal relationships between the discharge extremes of pairs of river stations.This illustrates our method on a real example for which we know the 'ground truth' of extremal causality, but unlike in the simulations of Section 4, we cannot control and do not know the true tail behaviour of the station discharges and their potential confounders.Figure 7: QQ-plot of the p mc estimates against the standard uniform distribution, with Kolmogorov-Smirnov confidence bands, for Pareto(1, 2) distributed ε 1 , ε 2 and H, for different causal strengths β 2,1 (colors), the three estimators (columns) and optional asymmetric confounding effects, β 1H = 0.8, β 2H = 1 (rows).

Data Sources and Additional Collection
We use the average daily discharges between January 1913 and December 2014 at the 68 Swiss gauging stations shown in Figure 1, and add daily precipitation data from 105 meteorological stations during the same period.Some additional information, such as the station elevation, catchment surface area and mean elevation, glaciation percent and coordinates, was collected from the Federal Office for the Environment's website.To reduce any seasonal effects due to unobserved confounders, we only consider data during June, July and August, as the more extreme observations happen during this period when mountain rivers are less likely to be frozen.Temporal clustering is likely to appear for average daily discharge data but can be captured by considering the average catchment precipitation as a covariate in the model for the GPD scale parameter (7).
Figure 8 shows relationships between the estimates, station altitudes and average discharges.Altitude does not greatly affect the estimates, but the shape parameter estimates broadly decrease with increased average river discharge volume.

Choice of Stations and Comonotonicity
For the causal analysis, we consider pairs of stations with known direct causal relationships, and pairs with no direct causal relationship.Causal pairs are ordered by the flow of water, with one downstream of the other.The river volumes for the pairs should be as similar as possible, as our exploratory analysis indicated different tail behaviours for rivers with very different average discharges.There should also be enough confluences between the two stations, otherwise one would observe comonotonicity, i.e., almost perfect dependence, between their discharges.If there is comonotonicity between X 1 and X 2 , then F 1 (X i,1 ) ≈ F 2 (X i,2 ), for all i = 1, . . ., n, and it is impossible to know which variable causes which based on the data alone regardless of the approach, even if one is certain both of direct causality and of its direction.Confluences between the two stations reduce comonotonicity and make it possible to detect the direction of causality.
The choice of covariate for the causal pairs was the mean daily precipitation among the meteorological stations in the area and the catchment of the two stations.The choice of covariate was less meaningful for the non-causal pairs with large separating distances, which have different meteorological conditions, so the average daily precipitation over the whole country was used.For the pair (42, 34), which has the closest stations and local precipitation data available, the daily average in the local catchments was also considered.In the latter case, the pair will be highlighted with an asterisk to avoid confusion.

Causal Analysis Results
For each station pair, the permutation test for direct causality was performed using the nonparametric (3) and H-conditional LGPD (8) estimators with post-fit correction or constraints, with R = 10 4 permutations and estimator hyper-parameters k = 1.5 n 0.4 and q = 0.9.Table 2 shows the values of p mc , the covariate shape estimate and its estimated extremal linear effects for the two stations, the latter estimated without constraints.The number of common observations for the pairs varies from 2,024 to 8,464, and k lies between 31 and 55.With precipitation covariates added, the number of common observations ranges from 1,483 to 7,820, and k lies between 27 and 54.
With the non-parametric approach for the causal stations, the absence of direct causality was rejected for four of the seven station pairs at significance level 5%, and for two of these four at Table 2: Permutation p-values p mc for station pairs using the non-parametric approach (NP), the H-conditional post-fit corrected (PFC) and constrained fit (CF) LGPD approaches, and an H-conditional exponential inverse-link GPD approach (Exp).The shape estimate ξH for the precipitation covariate and the unconstrained scale slope estimates are also shown (with standard errors of at most 0.03 for the former and in parentheses for the latter).With the non-parametric approach, the absence of direct causality was not rejected for ten of the 13 non-causal station pairs.Adding precipitation as a covariate with the two parametric approaches 'corrected' the p-value for another station.For the pair (42, 34) using local instead of global precipitation as a covariate gave a higher p-value.

Stations
We also considered using an exponential rather than a linear inverse-link function, i.e., taking log σ j (i) = σ 0 j + σ 1 j H i (i = 1, . . ., n; j = 1, 2), to avoid any need for correction or constraints.The resulting p mc values, also shown in Table 2, lead to the same conclusions as with the linear approaches.

Using the usual normal approximation, every σ1
1 is significantly positive for the causal pairs and 10 of the 14 estimates are positive for the non-causal pairs, with the highest confidence for the pair using local precipitation.Standard errors for σ1 2 are systematically larger than those for σ1 1 for the causal pairs, perhaps owing to the double causal effect of the covariate on the downstream station, both direct and indirect through the upstream station, as we do not observe this systematically for non-causal pairs.Consequently, the σ1 1 estimates are significantly positive for only four of the seven causal pairs, to be contrasted with 12 of the 14 estimates for the non-causal pairs.In particular, only the local precipitation effect is significant for the pair (42, 34).
We compare our results to two classical causal inference approaches appropriate to our problem.These are a non-Gaussian method for estimating causal linear structures based on results from independent component analysis, ICA-LiNGAM (Shimizu et al., 2006), and the PC algorithm, which retrieves the completed partially directed acyclic graph by performing conditional independence tests on the variables.For the latter, we consider both the classic PC algorithm (Spirtes et al., 2000), which uses Gaussian conditional independence tests, and the Rank PC algorithm (Harris and Drton, 2013), which uses rank-based Spearman correlation to perform the independence tests and thus is more robust to non-normal variables.The results for the ICA-LiNGAM method are presented in Table S.1 in the Supplementary Material, which shows the linear causal coefficients for the discharge station pairs estimated with the ICA-LiNGAM algorithm using either the station pair only (two variables) or the station pair and precipitation (three variables).Non-null values indicate significant causal effects.The upper-script arrows indicate the estimated direct causal direction between the station pair.Although in both cases of the two or three variables, ICA-LiNGAM retrieves all the correct causal pairs, with correct direction, all the non-causal pairs are indicated by non-null values as significantly causal.Both versions of the PC algorithm, once applied to our 21 pairs, provide existing direct causal links (without weights nor direction) between all the pairs of stations.Apparently both ICA-LiNGAM and PC methods are too eager to detect causality, unlike the tail coefficients.One explanation could be a set of unobserved confounders related to common global weather conditions triggering causal effects even between stations that are far apart.Extreme discharges depend more on local weather conditions, and particularly on heavy precipitation.Another explanation could be that causal effects are only linear in the tails, perhaps due to ground saturation by precipitation.

Discussion and Conclusion
This paper addresses the reduction or removal of the unwanted effect of known confounders from the extremal causal analysis between two variables and the discovery of extremal causal relationships using a parametric estimator of the causal tail coefficient, based on generalized Pareto modelling, and a permutation test for direct causality.Both allow the use of known confounders as covariates.
In our simulation study, the new estimator removed the confounder's unwanted effect almost entirely for variables with comparable tails, and reduced its effect enough to allow correct causal inference on the direct causal relationship in the case of a confounder with a heavier tail.The permutation test was shown to provide reliable p-values when all asymmetric confounding effects are captured in the model.
When applied to Swiss river discharge data, our methodology allowed correct inference on the direct causal relationships between discharges for the majority of the chosen station pairs, and the parametric approach captured the confounding effect of precipitation.
In many real-life situations, statistically significant covariates need not correspond to causal effects.Peters et al. (2016) have proposed a methodology for causal discovery, for when data from different settings or regimes are observed.Their method constructs invariant causal regression or classification models that should still make accurate predictions under interventions on the covariates or a change of environment.Adapting this approach to our setting would lead to a better understanding of causality of extremes.

S.2 Application Results for Competitors
Table S.1 shows the causal coefficients between the discharge station pairs estimated using ICA-LiNGAM, with and without considering the average catchment precipitation variable.

Figure 1 :
Figure 1: Topographic map of Switzerland showing the 68 gauging stations (red dots) along the Rhine, the Aare and their tributaries.Water flows towards station 68.Adapted from Engelke and Ivanovs (2021).

Figure 2 :
Figure 2: The six possible causal configurations between X 1 and X 2 with a possible confounder H, separated into the four cases studied in the simulations, and the two omitted by symmetry.

Figure 8 :
Figure 8: Relation between shape parameter estimates, scale parameter estimates (log scale), station elevation and average discharge (log scale), with standard errors (±SE) shown as error bars.
Adding daily precipitation as a covariate by either parametric approach decreases the p-values but two pairs remain non-significant; both lie in the same region and contain station 22.

Table S .
1: Linear causal coefficients for the discharge station pairs estimated with the ICA-LiNGAM algorithm using either the station pair only (LiNGAM, two variables) or the station pair and precipitation (LiNGAM-H, three variables).Non-null values indicate significant causal effects.The arrows indicate the estimated direct causal directions between the stations.