# Causal Modelling of Heavy-Tailed Variables and Confounders with Application
to River Flow

Olivier C. Pasche
Research Center for Statistics, University of Geneva, Switzerland,
Institute of Mathematics, EPFL, 1015 Lausanne, Switzerland,
<EMAIL_ADDRESS>
&Valérie Chavez-Demoulin
Faculty of Business and Economics, University of Lausanne, Switzerland,
<EMAIL_ADDRESS>
&Anthony C. Davison
Institute of Mathematics, EPFL, 1015 Lausanne, Switzerland,
<EMAIL_ADDRESS>

###### Abstract

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.

_Keywords_ Causation $\cdot$ Causal tail coefficient $\cdot$ Confounder
$\cdot$ Extreme value statistics $\cdot$ Generalized Pareto distribution

##  1 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.

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).

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

###  2.1 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},\dots,X_{p}$ satisfies

$X_{j}=\sum_{k\in\operatorname*{pa}(j)}\beta_{jk}X_{k}+\varepsilon_{j},\quad
j\in V,$

where $V:=\\{1,\dots,p\\}$ is a set of nodes representing the corresponding
random variables, $\operatorname*{pa}(j)\subseteq V$ is the set of parents of
$j$, $\beta_{jk}\in\mathbb{R}\setminus\\{0\\}$ is called the _causal weight_
of node $k$ on node $j$, and $\varepsilon_{1},\dots,\varepsilon_{p}$ are
jointly independent noise variables. We suppose that the _associated graph_
$G=(V,E)$, in which the directed edge $(i,j)\in V\times V$ belongs to $E$ if
and only if $i\in\operatorname*{pa}(j)$, is a directed acyclic graph (DAG).

In a DAG $G=(V,E)$, we say that $i\in V$ _is an ancestor of_ $j\in 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 $\operatorname*{An}(j,G)$, and we define
$\operatorname*{an}(j,G):=\operatorname*{An}(j,G)\setminus\\{j\\}$. In a LSCM
over random variables $X_{1},\dots,X_{p}$, with associated DAG $G=(V,E)$, we
say that $X_{i}$ _causes_ $X_{j}$, if $i\in\operatorname*{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
$\operatorname*{An}(i,G)\cap\operatorname*{An}(j,G)=\emptyset$. For any
$i,j\in V$ we let $\beta_{i\rightarrow j}$ denote the sum of the products of
the causal weights along the distinct directed paths from vertex $i$ to vertex
$j$; we set $\beta_{j\rightarrow j}:=1$ and $\beta_{i\rightarrow j}:=0$ if
$i\notin\operatorname*{An}(j,G)$.

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)

$\Gamma_{ij}:=\lim_{u\to 1^{-}}\mathbb{E}\left\\{F_{j}(X_{j})\mid
F_{i}(X_{i})>u\right\\},$ (1)

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}$, $\Gamma_{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, $\Gamma_{ji}$ will be appreciably
smaller than $\Gamma_{ij}$. As $\Gamma_{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.,

$\Psi_{ij}:=\lim_{u\to
1^{-}}\mathbb{E}\left[\rho\left\\{F_{j}(X_{j})\right\\}\mid\rho\left\\{F_{i}(X_{i})\right\\}>u\right]$

if the limit exists, where $\rho:x\mapsto\lvert 2x-1\rvert$. As
$F_{i}(X_{i})\sim{\rm Unif}(0,1)$,

$\Psi_{ij}=\underbrace{\lim_{u\to
1^{-}}\dfrac{1}{2}\mathbb{E}\left[\rho\left\\{F_{j}(X_{j})\right\\}\mid
F_{i}(X_{i})>u\right]}_{=:\Psi_{ij}^{+}}+\underbrace{\lim_{u\to
0^{+}}\dfrac{1}{2}\mathbb{E}\left[\rho\left\\{F_{j}(X_{j})\right\\}\mid
F_{i}(X_{i})<u\right]}_{=:\Psi_{ij}^{-}}.$

The interpretation and properties of $\Psi_{ij}$ are similar to those of
$\Gamma_{ij}$. The symmetric version captures the causal influence of $X_{i}$
on $X_{j}$ in both of their tails.

For simplicity we focus on $\Gamma_{ij}$ in this paper, though all of our
results and methods can be generalized to both tails by considering
$\Psi_{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 $\alpha\in\mathbb{R}$, written $f\in\operatorname*{RV}_{\alpha}$, if for
all $c>0$, $\lim_{x\to\infty}f(cx)/f(x)=c^{\alpha}$. If
$f\in\operatorname*{RV}_{0}$, then $f$ is said to be _slowly varying_.

######  Definition 2.3.

The random variable $X_{j}$ is said to be _regularly varying_ with index
$\alpha>0$, if, for some $\ell\in\operatorname*{RV}_{0}$,
$\mathbb{P}(X_{j}>x)\sim\ell(x)x^{-\alpha}$ as $x\to\infty$.

Independent regularly varying random variables $X_{1},\dots,X_{p}$ are said to
have _comparable upper tails_ if there exist $c_{1},\dots,c_{p}>0$, $\alpha>0$
and $\ell\in\operatorname*{RV}_{0}$ such that, for each $j\in\\{1,\dots,p\\}$,
$\mathbb{P}(X_{j}>x)\sim c_{j}\ell(x)x^{-\alpha}$ as $x\to\infty$.

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},\dots,X_{p}$ be random variables from a LSCM, with associated
directed acyclic graph $G=(V,E)$ and suppose that

  1. (a)

the coefficients $\beta_{jk}$ of the linear structural causal relationship
$X_{j}=\sum_{k\in\operatorname*{pa}(j,G)}\beta_{jk}X_{k}+\varepsilon_{j}$ are
strictly positive for all $j\in V$ and $k\in\operatorname*{pa}(j,G)$, and

  2. (b)

the real-valued noise variables $\varepsilon_{1},\dots,\varepsilon_{p}$ are
independent and regularly varying with comparable upper tails.

Then the values of $\Gamma_{ij}$ and $\Gamma_{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 $\Gamma_{ij}$ and $\Gamma_{ji}$ with the underlying causal relationship between $X_{i}$ and $X_{j}$. | $\Gamma_{ji}=1$ | $\Gamma_{ji}\in(1/2,1)$ | $\Gamma_{ji}=1/2$
---|---|---|---
$\Gamma_{ij}=1$ |  |  $X_{i}$ causes $X_{j}$ |
$\Gamma_{ij}\in(1/2,1)$ |  $X_{j}$ causes $X_{i}$ | common cause only |
$\Gamma_{ij}=1/2$ |  |  | no causal link

Under the theorem’s assumptions, the blank entries in Table 1 cannot occur.
Theorem 2.4 is generalizable to the $\Psi_{ij}$ variant of the coefficient and
possibly negative $\beta_{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\in V$, and with
$A_{ij}:=\operatorname*{An}(i,G)\cap\operatorname*{An}(j,G)$, is

$\Gamma_{ij}=\dfrac{1}{2}+\dfrac{1}{2}\dfrac{\sum_{h\in
A_{ij}}\beta_{h\rightarrow
i}^{\alpha}}{\sum_{h\in\operatorname*{An}(i,G)}\beta_{h\rightarrow
i}^{\alpha}}.$ (2)

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 $\left\\{(X_{i,1},X_{i,2})\right\\}_{i=1}^{n}$ 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 $\Gamma_{1,2}$ is defined to be

$\hat{\Gamma}_{1,2}=\dfrac{1}{k}\sum_{i=1}^{n}\hat{F}_{2}(X_{i,2})\mathds{1}(X_{i,1}>X_{(n-k),1})$
(3)

for some $k\in\\{1,\ldots,n-1\\}$, where $\mathds{1}(\cdot)$ denotes the
indicator function, $X_{(h),1}$ denotes the $h^{\text{th}}$ order statistic
and $\hat{F_{j}}$ is the empirical cumulative distribution function of
$X_{j}$, i.e.,

$\hat{F}_{j}(x)=\dfrac{1}{n}\sum_{i=1}^{n}\mathds{1}(X_{i,j}\leq x),\quad
j=1,2.$

This estimator is the empirical counterpart to (1), as
$X_{(h),1}=\hat{F}_{1}^{\leftarrow}(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\to\infty$, for a choice of $k$ such that $k\to\infty$ and
$k/n\to 0$ (Gnecco et al., 2021).

###  2.2 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
$\Gamma_{1,2}$, making it harder to tell whether a high value of
$\hat{\Gamma}_{1,2}$ indicates that $\Gamma_{1,2}=1$ or that
$\Gamma_{1,2}\lesssim 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 $\mathbf{H}$ of
potential confounders we therefore define

$\Gamma_{1,2\mid\mathbf{H}}:=\lim_{u\to
1^{-}}\mathbb{E}_{(X_{1},X_{2},\mathbf{H})}\left\\{F_{2}(X_{2}\mid\mathbf{H})\mid
F_{1}(X_{1}\mid\mathbf{H})>u\right\\}.$ (4)

If there is no direct dependence of $X_{2}$ on $X_{1}$, then $X_{2}$ is
independent of $X_{1}$ conditional on $\mathbf{H}$, so
$\Gamma_{1,2\mid\mathbf{H}}=1/2$, whereas $\Gamma_{1,2}$ lies in $[1/2,1)$ but
might be close to unity. Thus $\Gamma_{1,2\mid\mathbf{H}}<\Gamma_{1,2}$ unless
there are no confounders. If $X_{1}$ causes $X_{2}$, on the other hand, then
$\Gamma_{1,2\mid\mathbf{H}}=\Gamma_{1,2}=1$. In the presence of potential
confounders, therefore, (4) seems preferable to $\Gamma_{1,2}$. The difficulty
is that the estimation of (4) requires the modelling of the dependence of both
$X_{1}$ and $X_{2}$ on $\mathbf{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 $\mathbf{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 $\mathbf{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

###  3.1 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

$\mathbb{P}(X_{j}-u_{j}\leq x\mid X_{j}>u_{j})\approx
G(x;\sigma_{j},\xi_{j})=1-\left(1+\xi_{j}x/\sigma_{j}\right)^{-1/\xi_{j}}_{+},\quad
x>0,$ (5)

with a scale parameter $\sigma_{j}>0$ and a shape parameter
$\xi_{j}\in\mathbb{R}$:

  * •

$\xi_{j}=0$ corresponds to light-tailed distributions, and then $X_{j}$ lies
in the maximum domain of attraction of the Gumbel distribution;

  * •

$\xi_{j}>0$ corresponds to heavy-tailed distributions, and then $X_{j}$ lies
in the maximum domain of attraction of the Fréchet distribution; and

  * •

$\xi_{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 $\alpha>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\in(0,1)$, then we can write

$\begin{split}\mathbb{P}(X_{j}\leq
x)&\approx\left\\{G(x-u_{j};\sigma_{j},\xi_{j})(1-q)+q\right\\}\mathds{1}(x>u_{j})+\mathbb{P}(X_{j}\leq
x)\mathds{1}(x\leq u_{j}),\end{split}$

and using the empirical distribution $\hat{F}(x)$ to estimate
$\mathbb{P}(X_{j}\leq x)$ and maximum likelihood estimation using the excesses
of $u_{j}$ to obtain $\hat{\sigma}_{j}$ and $\hat{\xi}_{j}$ yields a hybrid
estimator of the distribution function $F_{j}(x)$ of $X_{j}$, i.e.,

$\hat{F}_{j}(x;\hat{\sigma}_{j},\hat{\xi}_{j})=\hat{F}(x)\mathds{1}(x\leq
u_{j})+\left\\{G(x-u_{j};\hat{\sigma}_{j},\hat{\xi}_{j})(1-q)+q\right\\}\mathds{1}(x>u_{j}).$

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\in\\{1,\ldots,n-1\\}$ yields the parametric _GPD causal tail
coefficient_ estimator for $\Gamma_{1,2}$,

$\hat{\Gamma}_{1,2}^{\rm
GPD}=\dfrac{1}{k_{g}}\sum_{i=1}^{n}\hat{F}_{2}(X_{i,2};\hat{\sigma}_{2},\hat{\xi}_{2})\mathds{1}\left\\{\hat{F}_{1}(X_{i,1};\hat{\sigma}_{1},\hat{\xi}_{1})>1-k/n\right\\},$
(6)

where
$k_{g}:=\lvert\\{i\in\\{1,\ldots,n\\}:\;\hat{F}_{1}(X_{i,1};\hat{\sigma}_{1},\hat{\xi}_{1})>1-k/n\\}\rvert$.
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
$\hat{F}_{1}(X_{i,1};\hat{\sigma}_{1},\hat{\xi}_{1})$.

The GPD model can be extended to allow dependence on covariates of interest by
expressing its parameters in the form
$\theta(i)=h\\{\boldsymbol{\gamma}^{\top}\mathbf{Z}(i)\\}$, where $\theta$
denotes one or both of $\sigma$ and $\xi$, $h$ is an inverse link function,
$\boldsymbol{\gamma}$ is a vector of parameters and $\mathbf{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
$\Gamma_{1,2}$ of a vector of potential confounders $\mathbf{H}$ of $X_{1}$
and $X_{2}$. If $\mathbf{H}$ is part of the LSCM then under the setup in
Section 2 it is straightforward to show that $\mathbf{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

$\sigma_{j}(i):=\sigma_{j}^{0}+\boldsymbol{\sigma}_{j}^{1\top}\mathbf{H}_{i},\quad
i=1,\ldots,n,\;j=1,2,$ (7)

where $\mathbf{H}_{i}$ is the replicate of $\mathbf{H}$ corresponding to the
observations $(X_{i,1},X_{i,2})$ of $(X_{1},X_{2})$.

This yields, for $k\in\\{1,\ldots,n-1\\}$, the parametric _$\mathbf{H}$
-conditional linear generalized Pareto distribution (LGPD) causal tail
coefficient estimator_,

$\hat{\Gamma}_{1,2\mid\mathbf{H}}^{\rm
GPD}=\dfrac{1}{k_{l}}\sum_{i=1}^{n}\hat{F}_{2}\\{X_{i,2};\hat{\sigma}_{2}(i),\hat{\xi}_{2}\\}\mathds{1}\left[\hat{F}_{1}\\{X_{i,1};\hat{\sigma}_{1}(i),\hat{\xi}_{1}\\}>1-k/n\right].$
(8)

where
$k_{l}:=\lvert\\{i\in\\{1,\ldots,n\\}:\;\hat{F}_{1}\\{X_{i,1};\hat{\sigma}_{1}(i),\hat{\xi}_{1}\\}>1-k/n\\}\rvert$.
Estimation of $\sigma_{j}^{0}$, $\sigma_{j}^{1}$ and $\xi_{j}$ is performed by
maximum likelihood. In applications it is preferable to center and rescale
each confounder in $\mathbf{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.

###  3.2 The Positive Linear Scale Issue

Linear modelling of the GPD scale parameter may not yield positive scale
estimates $\hat{\sigma}_{j}(i)>0$ for each $i=1,\ldots,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 $\mathbf{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 $\hat{\sigma}_{j}(i)$ in
(8) by $\max\\{\hat{\sigma}_{j}(i),\epsilon\\}$ for some arbitrary but small
positive $\epsilon$. The second solution, the _constrained approach_ , applies
the following linear constraints to the estimates when maximizing the
likelihood

$\sigma_{j}^{0}+\boldsymbol{\sigma}_{j}^{1\top}\min_{i=1,\ldots,n}\mathbf{H}_{i}>0,\quad\sigma_{j}^{0}+\boldsymbol{\sigma}_{j}^{1\top}\max_{i=1,\ldots,n}\mathbf{H}_{i}>0,\quad
j=1,2,$ (9)

where $\min_{i=1,\ldots,n}\mathbf{H}_{i}$ and
$\max_{i=1,\ldots,n}\mathbf{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 $X_{h}=H$ have $t_{\nu}$ distributions, then
$\sigma_{j}^{0}=u_{j}/\nu$ and $\sigma_{j}^{1}=-\beta_{h\rightarrow i}/\nu$.
Thus, if $\sigma_{j}(i)=\sigma_{j}^{0}+\sigma_{j}^{1}H_{i}>0$
$(j=1,2;i=1,\ldots,n)$, then

$-\dfrac{u_{j}}{\nu\max_{i=1,\ldots,n}H_{i}}<\sigma_{j}^{1}<-\dfrac{u_{j}}{\nu\min_{i=1,\ldots,n}H_{i}},$
(10)

where the lower and upper bounds are needed for positive and negative $H_{i}$,
respectively.

##  4 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 $\alpha=1/\xi>0$. We write
${\rm Pareto}(a,\alpha)$ for the Pareto model with scale parameter $a$ and
tail index $\alpha$; recal that lower values of $\alpha$ 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, ${\rm LogN}(\mu,\sigma^{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
$\beta_{2,1}$, $\beta_{1h}$ and $\beta_{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}$.

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.

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\lfloor n^{0.4}\rfloor=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 $\Gamma_{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.

###  4.1 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 $\mathbf{H}$-conditional LGPD
estimators, (3) or (8), of $\Gamma_{1,2}$ and $\Gamma_{2,1}$ alone. When the
causal effects of $H$ on $X_{1}$ and $X_{2}$, i.e., $\beta_{1h}$ and
$\beta_{2h}$, are increased relative to the noise variance and any causal
effect $\beta_{2,1}$ of $X_{1}$ on $X_{2}$, both $\Gamma_{1,2}$ and
$\Gamma_{2,1}$ increase in configuration B, and $\Gamma_{2,1}$ increases in
configurations C and D. This increase is larger with the non-parametric
estimators of $\Gamma_{1,2}$ and $\Gamma_{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 $\hat{\Gamma}_{1,2},\hat{\Gamma}_{2,1}\approx
1$ and hence $|\hat{\Gamma}_{1,2}-\hat{\Gamma}_{2,1}|\approx 0$ in
configurations B and D.

Use of the $\mathbf{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.

###  4.2 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 $\mathbf{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 $\hat{\Gamma}_{1,2}$ and $\hat{\Gamma}_{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.

Figure 3: Histograms of $\hat{\Gamma}_{1,2}$ (turquoise) and
$\hat{\Gamma}_{2,1}$ (blue) for $t_{4}$-distributed $\varepsilon_{1}$ and
$\varepsilon_{2}$, and $t_{3}$-distributed $H$ (top four panels) and for ${\rm
LogN}(0,1)$-distributed $\varepsilon_{1}$ and $\varepsilon_{2}$, and ${\rm
LogN}(0,1.5)$-distributed $H$ (bottom four panels). Half-lines (black)
indicate $\Gamma_{1,2}$ and $\Gamma_{2,1}$ for comparable tails. The panels
for ${\rm Pareto}(1,3)$ distributed $\varepsilon_{1}$ and $\varepsilon_{2}$,
and ${\rm Pareto}(1,1.5)$ distributed $H$ are very similar to the lower four
panels.

When $H$ has a heavier tail than $X_{1}$ and $X_{2}$, the non-parametric
estimators $\hat{\Gamma}_{1,2}$ and $\hat{\Gamma}_{2,1}$ in configuration B
and $\hat{\Gamma}_{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 $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ 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 $\hat{\Gamma}_{1,2\mid H}^{\rm
GPD}-\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ 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.

Figure 4: Histograms of $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ (turquoise) and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ (blue) with post-fit correction for
$t_{4}$ distributed $\varepsilon_{1}$ and $\varepsilon_{2}$, and $t_{3}$
distributed $H$ (top four panels), for ${\rm Pareto}(1,3)$ distributed
$\varepsilon_{1}$ and $\varepsilon_{2}$, and ${\rm Pareto}(1,1.5)$ distributed
$H$ (middle four panels), and ${\rm LogN}(0,1)$ distributed $\varepsilon_{1}$
and $\varepsilon_{2}$, and ${\rm LogN}(0,1.5)$ distributed $H$ (lower four
panels). Half-lines (black) indicate $\Gamma_{1,2}$ and $\Gamma_{2,1}$ for
comparable tails.

Simulation results for $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ 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 $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ and $\hat{\Gamma}_{2,1\mid H}^{\rm
GPD}$ with the constrained fit, for a heavier confounder tail.

Figure 5: Histograms of $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ (turquoise) and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ (blue) with constrained fit for $t_{4}$
distributed $\varepsilon_{1}$ and $\varepsilon_{2}$, and $t_{3}$ distributed
$H$. Half-lines (black) indicate $\Gamma_{1,2}$ and $\Gamma_{2,1}$ for
comparable tails.

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 $\hat{\sigma}_{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
$\hat{\sigma}_{j}^{1}$ may be up ten times smaller than those for post-fit
correction. With both approaches $\hat{\sigma}_{j}^{1}$ 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 $\mathbf{H}$-conditional 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 $\mathbf{H}$-conditional LGPD estimator allows
us to detect this relationship by reducing the effect of the confounding.

##  5 Testing for Direct Causality

###  5.1 Permutation Test

In situations such as the causal analysis presented in Section 6, the
distributions of the $\Gamma_{1,2}$ and $\Gamma_{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
$\left\\{(X_{i,1},X_{i,2})\right\\}_{i=1}^{n}$ from a LSCM and wish to test
the null hypothesis of no direct causal relationship between $X_{1}$ and
$X_{2}$, $H_{0}:\beta_{2,1}=0$, versus the alternative that $X_{1}$ causes
$X_{2}$, $H_{A}:\beta_{2,1}>0$. Our proposed procedure is as follows:

  1. 1.

Rescale values $\tilde{X}_{i,j}=\tilde{F}_{j}(X_{i,j})$ ($i=1,\ldots,n$,
$j=1,2$), where known confounders can be used in the distribution estimator
$\tilde{F}_{j}$, as for $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$.

  2. 2.

For $r=1,\ldots,R$, obtain $\tilde{X}_{i,1}^{(r)}$ and $\tilde{X}_{i,2}^{(r)}$
by randomly permuting the indices $j=1,2$ for each pair
$(\tilde{X}_{i,1},\tilde{X}_{i,2})$ ($i=1,\ldots,n$).

  3. 3.

Compute $\tilde{\Delta}_{1,2}=\tilde{\Gamma}_{1,2}-\tilde{\Gamma}_{2,1}$ on
the transformed original data
$\\{(\tilde{X}_{i,1},\tilde{X}_{i,2})\\}_{i=1}^{n}$ and
$\tilde{\Delta}^{*r}_{1,2}=\tilde{\Gamma}_{1,2}^{*r}-\tilde{\Gamma}_{2,1}^{*r}$
on their bootstrapped values
$\\{(\tilde{X}_{i,1}^{(r)},\tilde{X}_{i,2}^{(r)})\\}_{i=1}^{n}$
$(r=1,\ldots,R)$.

  4. 4.

Obtain the Monte Carlo $p$-value, by comparing the value of the test statistic
on the original rescaled data with the permutation distribution,

$p_{\rm
mc}=\frac{1+\\#_{r}\\{\tilde{\Delta}^{*r}_{1,2}\geq\tilde{\Delta}_{1,2}\\}}{R+1}.$

If there are no asymmetric confounding effects on the two variables, i.e.
$\beta_{1h}=\beta_{2h}$ in the case of a single confounder, then
$\Delta_{1,2}:=\Gamma_{1,2}-\Gamma_{2,1}=0$ under $H_{0}$, whereas
$\Delta_{1,2}>0$ under $H_{A}$; see equation (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
$\tilde{\Delta}_{1,2},\tilde{\Delta}^{*1}_{1,2},\ldots,\tilde{\Delta}^{*R}_{1,2}$
have a common distribution centered around zero, and $p_{\rm 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 $\tilde{\Delta}_{1,2}$ is more likely to lie in
the upper tail of
$\tilde{\Delta}^{*1}_{1,2},\ldots,\tilde{\Delta}^{*R}_{1,2}$. Thus the
distribution of $p_{\rm 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 $\tilde{F}_{j}$ 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}$.

###  5.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 $\beta_{2,1}$ of $X_{1}$ on $X_{2}$, with confounding effects
both present and absent. Symmetric ($\beta_{1H}=\beta_{2H}=1$) and asymmetric
($\beta_{1H}=0.8$ and $\beta_{2H}=1$, or $\beta_{1H}=1$ and $\beta_{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 $\mathbf{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\lfloor
n^{0.4}\rfloor=78$ and $q=0.9$.

Figure 6 shows uniform QQ-plots of $p_{\rm 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_{\rm
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$.

Figure 6: Uniform QQ-plots of Monte Carlo $p$-values $p_{\rm mc}$, with
Kolmogorov–Smirnov confidence bands for different causal strengths
$\beta_{2,1}$ (colors), the three estimators (columns) and optional symmetric
confounding effects, $\beta_{1H}=\beta_{2H}=1$ (rows). Top six panels: ${\rm
Pareto}(1,2)$ distributed $\varepsilon_{1}$ and $\varepsilon_{2}$, and ${\rm
Pareto}(1,1)$ distributed $H$. Bottom six panels: $t_{4}$ distributed
$\varepsilon_{1}$ and $\varepsilon_{2}$, and $t_{3}$ distributed $H$.

When the confounding effects are added, the test based on the non-parametric
estimator fails for the Pareto distribution, as most of the $p_{\rm mc}$ then
lie outside the $95\%$ confidence bands, indicating that the distribution of
$p_{\rm 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_{\rm 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.

Figure 7: QQ-plot of the $p_{\rm mc}$ estimates against the standard uniform
distribution, with Kolmogorov–Smirnov confidence bands, for ${\rm
Pareto}(1,2)$ distributed $\varepsilon_{1}$, $\varepsilon_{2}$ and $H$, for
different causal strengths $\beta_{2,1}$ (colors), the three estimators
(columns) and optional asymmetric confounding effects, $\beta_{1H}=0.8$,
$\beta_{2H}=1$ (rows).

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_{\rm
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.

##  6 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.

###  6.1 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.

Figure 8: Relation between shape parameter estimates, scale parameter
estimates (log scale), station elevation and average discharge (log scale),
with standard errors ($\pm{\rm SE}$) shown as error bars.

###  6.2 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})\approx F_{2}(X_{i,2})$, for
all $i=1,\ldots,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.

As we shall use precipitation as the confounding covariate, the stations must
share likely meteorological effects and must lie in regions where
precipitation data is available. Based on these criteria, we chose seven
causal station pairs: $(43,62)$, $(42,63)$, $(36,63)$, $(24,61)$, $(44,61)$,
$(22,38)$, $(22,35)$, where the first station of each pair lies upstream from
the second.

The non-causal station pairs were selected to have similar average volume and
similar shape parameter estimates. Pairs with stations separated by long
distances and pairs relatively close to each other were both considered. The
$13$ pairs selected are $(30,45)$, $(36,39)$, $(42,34)$, $(32,33)$, $(62,63)$,
$(57,60)$, $(13,14)$, $(17,22)$, $(12,21)$, $(26,28)$, $(27,31)$, $(23,39)$,
$(23,35)$.

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.

###  6.3 Causal Analysis Results

For each station pair, the permutation test for direct causality was performed
using the non-parametric (3) and $\mathbf{H}$-conditional LGPD (8) estimators
with post-fit correction or constraints, with $R=10^{4}$ permutations and
estimator hyper-parameters $k=1.5\lfloor n^{0.4}\rfloor$ and $q=0.9$. Table
6.3 shows the values of $p_{\rm 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 $\numprint{2024}$ to $\numprint{8464}$, and $k$ lies between $31$ and
$55$. With precipitation covariates added, the number of common observations
ranges from $\numprint{1483}$ to $\numprint{7820}$, and $k$ lies between $27$
and $54$.

Table 2: Permutation $p$-values $p_{\rm mc}$ for station pairs using the non-
parametric approach (NP), the $\mathbf{H}$-conditional post-fit corrected
(PFC) and constrained fit (CF) LGPD approaches, and an
$\mathbf{H}$-conditional exponential inverse-link GPD approach (Exp). The
shape estimate $\hat{\xi}_{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).

Stations | Pair type |  NP |  PFC |  CF |  Exp |  ^ξH |  ^σ11 |  ^σ21
---|---|---|---|---|---|---|---|---
43-62 | causal | $\numprint{0.0080991901}$ | $\numprint{0.0079992001}$ | $\numprint{0.0074992501}$ | $\numprint{0.0080991901}$ | $\numprint{0.0566669337}$ |  $\numprint{0.8842009468}$(0.3) |  $\numprint{1.9066945453}$(1.3)
42-63 | causal | $\numprint{0.0279972003}$ | $\numprint{0.0211978802}$ | $\numprint{0.0245975402}$ | $\numprint{0.0351964804}$ | $\numprint{0.0601108523}$ |  $\numprint{6.4878543848}$(1.1) |  $\numprint{8.6026410091}$(2.2)
36-63 | causal | $\numprint{0.0274972503}$ | $\numprint{0.0226977302}$ | $\numprint{0.0198980102}$ | $\numprint{0.0295970403}$ | $\numprint{0.0601108523}$ |  $\numprint{5.027641044}$(1.1) |  $\numprint{7.2532593647}$(2.8)
24-61 | causal | $\numprint{0.0555944406}$ | $\numprint{0.0074992501}$ | $\numprint{0.0053994601}$ | $\numprint{0.0039996}$ | $\numprint{-0.014601569}$ |  $\numprint{3.4157640955}$(1.2) |  $\numprint{-2.3350882618}$(2.4)
44-61 | causal | $\numprint{0.0104989501}$ | $\numprint{0.00449955}$ | $\numprint{0.00459954}$ | $\numprint{0.0065993401}$ | $\numprint{0.0127824624}$ |  $\numprint{1.8939044269}$(0.7) |  $\numprint{-1.2100340565}$(2.0)
22-38 | causal | $\numprint{0.5802419758}$ | $\numprint{0.402859714}$ | $\numprint{0.397660234}$ | $\numprint{0.3305669433}$ | $\numprint{0.0659749896}$ |  $\numprint{3.425892068}$(0.8) |  $\numprint{8.000950484}$(2.0)
22-35 | causal | $\numprint{0.2177782222}$ | $\numprint{0.1730826917}$ | $\numprint{0.1745825417}$ | $\numprint{0.097190281}$ | $\numprint{0.0295179206}$ |  $\numprint{3.4292170544}$(0.9) |  $\numprint{11.666130149}$(3.0)
30-45 | non-caus. | $\numprint{0.5578442156}$ | $\numprint{0.4658534147}$ | $\numprint{0.4668533147}$ | $\numprint{0.4561543846}$ | $\numprint{0.0052839693}$ |  $\numprint{1.0132876404}$(0.4) |  $\numprint{0.8868363734}$(0.9)
36-39 | non-caus. | $\numprint{0.796620338}$ | $\numprint{0.699030097}$ | $\numprint{0.698330167}$ | $\numprint{0.6874312569}$ | $\numprint{0.0052839693}$ |  $\numprint{4.6059105148}$(1.1) |  $\numprint{4.168906454}$(1.6)
42-34 | non-caus. | $\numprint{0.2264773523}$ | $\numprint{0.0422957704}$ | $\numprint{0.0408959104}$ | $\numprint{0.099690031}$ | $\numprint{0.0052839693}$ |  $\numprint{5.967653909}$(1.2) |  $\numprint{0.4274070406}$(0.3)
42-34∗ | non-caus. | $\numprint{0.2264773523}$ | $\numprint{0.1303869613}$ | $\numprint{0.1250874913}$ | $\numprint{0.1125887411}$ | $\numprint{0.0543840826}$ |  $\numprint{6.2882944671}$(1.1) |  $\numprint{0.6556023806}$(0.3)
32-33 | non-caus. | $\numprint{0.0098990101}$ | $\numprint{0.0096990301}$ | $\numprint{0.0099990001}$ | $\numprint{0.00159984}$ | $\numprint{0.0052839693}$ |  $\numprint{0.6325946323}$(0.4) |  $\numprint{1.0012907206}$(0.3)
62-63 | non-caus. | $\numprint{0.101389861}$ | $\numprint{0.4889511049}$ | $\numprint{0.4753524648}$ | $\numprint{0.299970003}$ | $\numprint{0.0052839693}$ |  $\numprint{1.08066753}$(1.4) |  $\numprint{7.6742023274}$(2.1)
57-60 | non-caus. | $\numprint{0.9923007699}$ | $\numprint{0.99930007}$ | $\numprint{0.99910009}$ | $\numprint{0.9990001}$ | $\numprint{0.0052839693}$ |  $\numprint{6.3100949401}$(3.7) |  $\numprint{5.2283442628}$(1.8)
13-14 | non-caus. | $\numprint{0.3215678432}$ | $\numprint{0.5590440956}$ | $\numprint{0.5639436056}$ | $\numprint{0.5337466253}$ | $\numprint{0.0052839693}$ |  $\numprint{0.5925547817}$(0.2) |  $\numprint{1.1850691043}$(0.3)
17-22 | non-caus. | $\numprint{0.0050994901}$ | $\numprint{0.0547945205}$ | $\numprint{0.0565943406}$ | $\numprint{0.0478952105}$ | $\numprint{0.0052839693}$ |  $\numprint{0.776163342}$(0.5) |  $\numprint{2.1809668187}$(0.7)
12-21 | non-caus. | $\numprint{0.5146485351}$ | $\numprint{0.498550145}$ | $\numprint{0.500349965}$ | $\numprint{0.7209279072}$ | $\numprint{0.0052839693}$ |  $\numprint{0.7126162033}$(0.3) |  $\numprint{1.3336258551}$(0.4)
26-28 | non-caus. | $\numprint{0.6345365463}$ | $\numprint{0.900109989}$ | $\numprint{0.8927107289}$ | $\numprint{0.9220077992}$ | $\numprint{0.0052839693}$ |  $\numprint{1.9041701392}$(0.5) |  $\numprint{1.6324917792}$(0.4)
27-31 | non-caus. | $\numprint{0.404459554}$ | $\numprint{0.6259374063}$ | $\numprint{0.6215378462}$ | $\numprint{0.7521247875}$ | $\numprint{0.0052839693}$ |  $\numprint{1.7129163243}$(0.7) |  $\numprint{2.9055730916}$(1.1)
23-39 | non-caus. | $\numprint{0.802719728}$ | $\numprint{0.9134086591}$ | $\numprint{0.9184081592}$ | $\numprint{0.9323067693}$ | $\numprint{0.0052839693}$ |  $\numprint{2.4975701839}$(0.6) |  $\numprint{4.2681161898}$(1.5)
23-35 | non-caus. | $\numprint{0.6484351565}$ | $\numprint{0.8844115588}$ | $\numprint{0.8877112289}$ | $\numprint{0.8647135286}$ | $\numprint{0.0052839693}$ |  $\numprint{2.4975701839}$(0.6) |  $\numprint{6.6610856116}$(1.7)

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 level $2.5\%$. 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$.

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.

We also considered using an exponential rather than a linear inverse-link
function, i.e., taking $\log\sigma_{j}(i)=\sigma_{j}^{0}+\sigma_{j}^{1}H_{i}$
$(i=1,\ldots,n;j=1,2)$, to avoid any need for correction or constraints. The
resulting $p_{\rm mc}$ values, also shown in Table 6.3, lead to the same
conclusions as with the linear approaches.

Using the usual normal approximation, every $\hat{\sigma}_{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 $\hat{\sigma}_{2}^{1}$ are
systematically larger than those for $\hat{\sigma}_{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
$\hat{\sigma}_{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.

##  7 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.

## Acknowledgements

The work was supported by the Swiss National Science Foundation.

## References

  * Coles (2001) Coles, S.  _An Introduction to Statistical Modeling of Extreme Values_.  Springer London, 2001.  doi:10.1007/978-1-4471-3675-0.
  * Davison and Hinkley (1997) Davison, A. C. and Hinkley, D. V.  _Bootstrap Methods and their Application_.  Cambridge University Press, New York, 1997.  doi:10.1017/CBO9780511802843.
  * Davison and Smith (1990) Davison, A. C. and Smith, R. L.  Models for exceedances over high thresholds (with discussion).  _Journal of the Royal Statistical Society, Series B_ , 52:393–442, 1990.  doi:10.1111/j.2517-6161.1990.tb01796.x.
  * Engelke and Ivanovs (2021) Engelke, S. and Ivanovs, J.  Sparse structures for multivariate extremes.  _Annual Review of Statistics and its Application_ , 8:241–270, 2021.  doi:10.1146/annurev-statistics-040620-041554.
  * Gissibl and Klüppelberg (2018) Gissibl, N. and Klüppelberg, C.  Max-linear models on directed acyclic graphs.  _Bernoulli_ , 24:2693–2720, 2018.  doi:10.3150/17-BEJ941.
  * Gnecco et al. (2021) Gnecco, N., Meinshausen, N., Peters, J., and Engelke, S.  Causal discovery in heavy-tailed models.  _Annals of Statistics_ , 49(3):1755–1778, 2021\.  doi:10.1214/20-AOS2021.
  * Harris and Drton (2013) Harris, N. and Drton, M.  PC Algorithm for Nonparanormal Graphical Models.  _Journal of Machine Learning Research_ , 14:3365–3383, 2013.
  * Kiriliouk and Naveau (2020) Kiriliouk, A. and Naveau, P.  Climate extreme event attribution using multivariate peaks-over-thresholds modeling and counterfactual theory.  _Annals of Applied Statistics_ , 14(3):1342–1358, 2020.  doi:10.1214/20-AOAS1355.
  * Klüppelberg and Krali (2021) Klüppelberg, C. and Krali, M.  Estimating an extreme bayesian network via scalings.  _Journal of Multivariate Analysis_ , 181:104672, 2021.  doi:10.1016/j.jmva.2020.104672.
  * Maathuis and Nandy (2016) Maathuis, M. H. and Nandy, P.  _A review of some recent advances in causal inference. In P. Bhlmann, P. Drineas, M. Kane and M.J. van der Laan (Eds.), Handbook of Big Data_.  Chapman and Hall, 2016.
  * Mhalla et al. (2020) Mhalla, L., Chavez-Demoulin, V., and Dupuis, D.  Causal mechanism of extreme river discharges in the upper Danube basin network.  _Applied Statistics_ , 69:741–764, 2020.  doi:10.1111/rssc.12415.
  * Naveau et al. (2020) Naveau, P., Hannart, A., and Ribes, A.  Statistical methods for extreme event attribution in climate science.  _Annual Review of Statistics and Its Application_ , 7:89–110, 2020.  doi:10.1146/annurev-statistics-031219-041314.
  * Pearl (2009) Pearl, J.  _Causality: Models, Reasoning, and Inference_.  Cambridge University Press, New York, NY, USA, 2nd edition, 2009.
  * Peters et al. (2016) Peters, J., Bühlmann, P., and Meinshausen, N.  Causal inference by using invariant prediction: identification and confidence intervals (with Discussion).  _Journal of the Royal Statistical Society, Series B_ , 78(5):947–1012, 2016.  doi:10.1111/rssb.12167.
  * Peters et al. (2017) Peters, J., Janzing, D., and Schölkopf, B.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  MIT Press, Cambridge, MA, 2017.
  * Shimizu et al. (2006) Shimizu, S., Hoyer, P. O., Hyvärinen, A., and Kerminen, A.  A Linear Non-Gaussian Acyclic Model for Causal Discovery.  _Journal of Machine Learning Research_ , 7:2003–2030, 2006.
  * Spirtes et al. (2000) Spirtes, P., Glymour, C., and Scheines, R.  _Causation, Prediction, and Search_.  MIT press, Cambridge, MA, USA, 2000.

SUPPLEMENTARY MATERIAL TO

“Causal Modelling of Heavy-Tailed Variables and Confounders with Application
to River Flow”

Olivier C. Pasche
---
Research Center for Statistics, University of Geneva, Switzerland,
Institute of Mathematics, EPFL, 1015 Lausanne, Switzerland,
<EMAIL_ADDRESS>
Valérie Chavez-Demoulin
---
Faculty of Business and Economics, University of Lausanne, Switzerland,
<EMAIL_ADDRESS>
Anthony C. Davison
---
Institute of Mathematics, EPFL, 1015 Lausanne, Switzerland,
<EMAIL_ADDRESS>

##  Appendix S.1 Variables with Comparable Tails

###  S.1.1 Non-Parametric Causal Tail Coefficient Estimator

Figure 6.3 shows the sample distributions of the non-parametric estimators
$\hat{\Gamma}_{1,2}$ and $\hat{\Gamma}_{2,1}$ for all four causal structures,
for the $t_{4}$, ${\rm Pareto}(1,2)$ and ${\rm LogN}(0,1)$ noise
distributions, respectively. The true coefficient values $\Gamma_{1,2}$ and
$\Gamma_{2,1}$ are obtained using (2). Figure S.2 shows the sample
distribution of the coefficient difference estimator
$\hat{\Delta}_{1,2}:=\hat{\Gamma}_{1,2}-\hat{\Gamma}_{2,1}$ for the $t_{4}$
noise distribution.

Figure S.1: Histograms of $\hat{\Gamma}_{1,2}$ (turquoise) and
$\hat{\Gamma}_{2,1}$ (blue) for $t_{4}$ (top four panels), ${\rm Pareto}(1,2)$
(middle four panels) and ${\rm LogN}(0,1)$ (bottom four panels) distributed
noise variables, for the four causal configurations. Half-lines (black)
indicate $\Gamma_{1,2}$ and $\Gamma_{2,1}$.

Figure S.2: Histogram of $\hat{\Delta}_{1,2}$ for $t_{4}$ distributed noise
variables, for the four causal configurations. Lines indicate
$\Delta_{1,2}=\Gamma_{1,2}-\Gamma_{2,1}$.

###  S.1.2 LGPD Causal Tail Coefficient with Post-Fit and Constrained Fit
Corrections

Figure S.3 shows the sample distribution of $\hat{\Gamma}_{1,2\mid H}^{\rm
GPD}$ and $\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ with the constrained fit, for a
comparable confounder tail. Figure 6.3 shows the sample distribution of
$\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ and $\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$
with post-fit correction for all four causal configurations, for the $t_{4}$,
${\rm Pareto}(1,2)$ and ${\rm LogN}(0,1)$ noise distributions.

Figure S.3: Histograms of $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ (turquoise) and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ (blue) with constrained fit for $t_{4}$
distributed noise variables, for the four causal configurations. Half-lines
(black) indicate $\Gamma_{1,2}$ and $\Gamma_{2,1}$.

Figure S.4: Histograms of $\hat{\Gamma}_{1,2\mid H}^{\rm GPD}$ (turquoise) and
$\hat{\Gamma}_{2,1\mid H}^{\rm GPD}$ (blue) with post-fit correction for
$t_{4}$ (top four panels), ${\rm Pareto}(1,2)$ (middle four panels) and ${\rm
LogN}(0,1)$ (bottom four panels) distributed noise variables, for the four
causal configurations. Half-lines (black) indicate $\Gamma_{1,2}$ and
$\Gamma_{2,1}$.

##  Appendix 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.

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. Stations | Pair type |  LiNGAM |  LiNGAM-$H$
---|---|---|---
43-62 | causal |  $\numprint{1.91929474292931}$→ |  $\numprint{2.02147993659115}$→
42-63 | causal |  $\numprint{2.07879186262232}$→ |  $\numprint{2.20875242829659}$→
36-63 | causal |  $\numprint{3.28936394744567}$→ |  $\numprint{3.60555863168095}$→
24-61 | causal |  $\numprint{2.96189249513038}$→ |  $\numprint{3.02797522258647}$→
44-61 | causal |  $\numprint{2.66344886172686}$→ |  $\numprint{2.83230029304895}$→
22-38 | causal |  $\numprint{2.35242595575385}$→ |  $\numprint{2.35242595575385}$→
22-35 | causal |  $\numprint{2.54841558325143}$→ |  $\numprint{2.54841558325143}$→
30-45 | non-caus. |  $\numprint{0.844373945195069}$→ |  $\numprint{0.87138638378129}$→
36-39 | non-caus. |  $\numprint{0.6647466}$← |  $\numprint{0.6563216}$←
42-34 | non-caus. |  $\numprint{1.392225}$← |  $\numprint{1.2881381}$←
42-34∗ | non-caus. |  $\numprint{1.392225}$← |  $\numprint{1.392225}$←
32-33 | non-caus. |  $\numprint{0.593040513333397}$→ |  $\numprint{0.535594175755528}$→
62-63 | non-caus. |  $\numprint{1.02472643450335}$→ |  $\numprint{1.04797361573715}$→
57-60 | non-caus. |  $\numprint{0.68239595737648}$→ |  $\numprint{0.673070717007456}$→
13-14 | non-caus. |  $\numprint{0.5025991}$← |  $\numprint{1.10324474969347}$→
17-22 | non-caus. |  $\numprint{1.79725395000552}$→ |  $\numprint{1.69241757070694}$→
12-21 | non-caus. |  $\numprint{1.03944101513775}$→ |  $\numprint{1.07511701685805}$→
26-28 | non-caus. |  $\numprint{0.7487394}$← |  $\numprint{0.7193701}$←
27-31 | non-caus. |  $\numprint{0.535030227721156}$→ |  $\numprint{0.660891637427842}$→
23-39 | non-caus. |  $\numprint{0.254151192169057}$→ |  $\numprint{0.183805221682078}$→
23-35 | non-caus. |  $\numprint{0.423726060357545}$→ |  $\numprint{0.356553057370821}$→