# Detecting causal covariates for extreme dependence structures

Juraj Bodík1111 Email<EMAIL_ADDRESS>Linda Mhalla1,2, Valérie Chavez-
Demoulin1

(1 Faculté des Hautes Études Commerciales, Université de Lausanne, Switzerland
2Institute of Mathematics, École polytechnique fédérale de Lausanne,
Switzerland
)

###### Abstract

Determining the causes of extreme events is a fundamental question in many
scientific fields. An important aspect when modelling multivariate extremes is
the tail dependence. In application, the extreme dependence structure may
significantly depend on covariates. As for the general case of modelling
including covariates, only some of the covariates are causal. In this paper,
we propose a methodology to discover the causal covariates explaining the tail
dependence structure between two variables. The proposed methodology for
discovering causal variables is based on comparing observations from different
environments or perturbations. It is a desired methodology for predicting
extremal behaviour in a new, unobserved environment. The methodology is
applied to a dataset of $\text{NO}_{2}$ concentration in the UK. Extreme
$\text{NO}_{2}$ levels can cause severe health problems, and understanding the
behaviour of concurrent severe levels is an important question. We focus on
revealing causal predictors for the dependence between extreme $\text{NO}_{2}$
observations at different sites.

## Introduction

Modelling extreme events is a crucial part of many real-world problems and the
focus is often on the extremal dependence structure between several variables
(Coles, 2001). In a univariate case, under very general assumptions, a
normalized maxima of random variables converges to a Generalized Extreme Value
(GEV) distribution (Gnedenko, 1943). In the multivariate case, however, the
limiting dependence structures for extreme values do not follow a simple
parametric form. Therefore, alternative joint tail representations should be
used, and extensive literature describes different approaches (Resnick, 1987;
Beirlant et al., 2006).

More recent literature describes the modelling and estimation of covariate
influence on the joint tail structure. Mhalla et al. (2019b) propose a
parametric framework for integrating covariate information through the angular
density model. Non-parametric modelling of the joint tail structure using GAM
modelling (Wood, 2017) has been developed in Mhalla et al. (2017). In this
paper, our main emphasis is on the causal properties of the extremal tail
coefficient between two random variables $Z_{1},Z_{2}$ (Joe, 1997).

A real-world problem of interest consists of $\text{NO}_{2}$ measurements in
the UK. Features of $\text{NO}_{2}$ concentration levels or other air
pollutants are of interest in environmental science (Shi et al., 2014) as
extreme values of $\text{NO}_{2}$ can be very dangerous and have been linked
to an increased risk of lung cancer (Hamra et al., 2018). The dataset taken
from Carslaw and Ropkins (2022) comprises 241 measurement stations in the UK
with hourly records of $\text{NO}_{2}$ [µg/m3]. We use only the stations with
at least 15 years of measurements (21 stations in total). $\text{NO}_{2}$ is
partially produced by burning fossil fuel and motor vehicle exhaust. Hence, we
distinguish traffic stations where the $\text{NO}_{2}$ level is predominantly
determined by nearby traffic (12 stations) from background stations (9
stations). Figure 1 shows $\text{NO}_{2}$ concentration (in $\mu$g/m3) for one
traffic station and one background station over time. We observe that the
magnitude is usually larger in traffic stations than in background stations.
Although a yearly seasonality seems clearer for the background station, both
series highlight a long term decreasing trend. The point of our interest is
the extremal dependence between a pair of stations explained by some
covariates. Some of these covariates may be causal, others not. We are
interested in answering the following question: What is the cause of the
extremal dependence between $\text{NO}_{2}$ concentrations at two stations?

Figure 1: An example of traffic and background stations in London.

Determining a system’s causal structure is an important problem in many
scientific fields. There is a growing literature on causal inference under
various types of assumptions and different frameworks, with applications
ranging from public health to biology and economics (Chang et al., 2015;
Varian, 2016). Often, the goal is to find the causes of some response variable
$Y$ among a given set of covariates X or to quantify the causal relationships
between a set of variables. In recent years, a causal discovery based on
conditional invariance has become an important growing research area. Starting
with a groundwork paper (Peters et al., 2016), connections between
distributional invariance and causality started to emerge. The reviewing paper
Bühlmann (2020) summarises the groundwork and highlights the importance of
causal prediction. The methodology based on invariant causal prediction (ICP)
relies on observing a system under different environments (Pearl, 2009). Then,
one can observe some invariant properties for causal variables through all
environments, but not for non-causal ones.

Here, we are interested in unveiling causal links explaining extremal
dependence between random variables. Precisely, we want to assess whether an
intervention on a covariate affects the tail dependence of the random
variables. Classical causal techniques that usually give a good description of
data tend to fail in the upper tail due to the scarcity of observations
(Engelke and Ivanovs, 2021). A combination of causality and extreme value
theory (EVT) is a developing research area. Mhalla et al. (2019a) use a
Kolmogorov complexity to detect causal direction using observational data.
Kiriliouk and Naveau (2020) studied probabilities of necessary and sufficient
causation as defined in the counterfactual theory using the multivariate
generalized Pareto distributions. Gnecco et al. (2020) define a causal tail
coefficient revealing the causal structure between two random variables by
capturing asymmetries in the extremal dependence of the two random variables
under the context of the linear structural causal model. Their methodology
holds also in the presence of latent common causes with the same tail
heaviness as the observed variables. Pasche et al. (2021) proposes a
methodology relying on the causal tail coefficient that can deal with latent
common causes that have different tail indices. Deuber et al. (2021) describe
a framework for estimating extremal quantile treatment effects. However, to
our knowledge, no theory combines EVT and causal inference using the ICP
framework and inferring causes of the extremal dependence structure.

A well-established approach in causal inference relies on the structural
causal model (SCM) of Pearl (2009). The ICP framework exploits invariance
properties in SCM for inference; see Haavelmo (1944); Dawid (2006); Schölkopf
et al. (2012) for similar notions. The main emphasis in the ICP setting has
been on linear regression models assuming the additive form
$Y=\textbf{X}_{S^{\star}}\beta+\varepsilon$, where $Y$ is the response
variable of interest, $\textbf{X}_{S^{\star}}$ is the set of direct causal
covariates among the set of covariates X and
$S^{\star}\subseteq\\{1,\dots,d\\}$ the set of indices of direct causes of
$Y$. Generalizations of this model have been given by Heinze-Deml et al.
(2018) with the more flexible form $Y=g(\textbf{X}_{S^{\star}})+\varepsilon$
where $g$ is a smooth, non-parametric function. We aim to use the ICP
methodology in the context of the joint tail structure.

The paper is organized as follows: Section 1 introduces EVT in the bivariate
context. Section 2 describes the ICP methodology. Section 3 proposes a
methodology to unveil causal covariates among a set of covariates explaining
the extremal dependence between two random variables. A simulation study is
presented in Section 4 and Section 5 illustrates the proposed methodology
through the dataset of $\text{NO}_{2}$ measurements in the UK.

##  1 Bivariate extreme value theory

Let $(Y_{i})_{i\geq 1}$ be a sequence of independent and identically
distributed (iid) random variables with common distribution $F$. We are
interested in the behaviour of the maximum of a sequence of $n$ such random
variables, that is $M_{n}=\max_{i=1,\dots,n}Y_{i}$. If there exist a non-
degenerate distribution $G$ and sequences of constants $\\{a_{n}>0\\}$ and
$\\{b_{n}\\}$ such that $\frac{M_{n}-b_{n}}{a_{n}}\overset{D}{\to}G$ as
$n\to\infty$, then $G$ belongs to the Generalized Extreme Value (GEV)
distributions family. The GEV distribution has three parameters, namely a
location $\mu\in\mathbb{R}$, a scale $\sigma>0$, and a shape
$\xi\in\mathbb{R}$. Its distribution function takes the form

$\displaystyle F(x)$
$\displaystyle=\exp\left\\{-\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-1/\xi}\right\\}\text{,
}\xi\neq 0,$ (1) $\displaystyle F(x)$
$\displaystyle=\exp\left\\{-\exp\left(-\frac{x-\mu}{\sigma}\right)\right\\}\text{,
}\xi=0,$

defined on the support
$[\mu-\sigma/\xi,\infty),(-\infty,\infty),(-\infty,\mu-\sigma/\xi]$ for cases
$\xi<0,\xi=0,\xi>0$ respectively. Cases when $\xi>0,\xi=0$, and $\xi<0$
correspond to the Fréchet, Gumbel, and Weibull distributions, respectively
(Fisher and Tippett, 1928; Gnedenko, 1943).

Suppose now $(\textbf{Y}_{i})_{i\geq 1}$ is a sequence of independent copies
of a two dimensional random vector $\textbf{Y}=(Y_{1},Y_{2})^{\top}$ with
marginal distributions $F_{j}$, $j=1,2$ and joint distribution $F$ and denote
by $\textbf{M}_{n}=(M_{n,1},M_{n,2})^{\top}$ their component-wise maxima. We
suppose that there exist sequences
$\\{\textbf{a}_{n}\\}\subset\mathbb{R}^{2}_{+}$ and
$\\{\textbf{b}_{n}\\}\subset\mathbb{R}^{2}$ such that the normalized vector
$(\textbf{M}_{n}-\textbf{b}_{n})/\textbf{a}_{n}$ converges in distribution to
a random vector Z with joint distribution $G$ and non-degenerate margins
$G_{j}$, $j=1,2$. Then, the limiting distribution $G$ is called the
multivariate extreme value distribution (MEVD) (Resnick, 1987, Chapter 5.3).
The class of MEVDs coincides with the class of max-stable distribution
functions with non-degenerate margins. The margins $G_{j}$ are max-stable, or
equivalently GEV. Without loss of generality, it is convenient to transform
the marginal distributions $G_{j}$ to unit Fréchet margins with distribution
function $P(Z_{j}\leq z)=\exp{(-1/z)}1_{[0,\infty)}(z)$. The transformation is
marginally performed using a probability integral transform with

$t_{j}(z)=-1/\log{G_{j}(z)},z>0,$ (2)

applied to the $j$th GEV margin $G_{j}$. When the margins are transformed to
unit Fréchet, the dependence structure of the transformed random vector is the
same as the dependence structure of the vector Z. We therefore assume
hereafter that Z has unit Fréchet margins. In that case, $G$ is called a
simple max-stable distribution and has the representation
$G(\textbf{x})=e^{-V(\textbf{x})}$, where $V(\textbf{x})$ is called the
exponent function (Resnick, 1987, Chapter 5.4). It can be shown that $V$ is
$-1$ homogeneous, that is,
$V(t\textbf{x})=\frac{1}{t}V(\textbf{x}),t>0,\textbf{x}=(x_{1},x_{2})\in\mathbb{R}^{2}_{+}$
(Coles and Tawn, 1991). The Pickands dependence function (Pickands, 1981)
$A:\mathcal{S}_{2}\to[\frac{1}{2},1]$ is defined by

$V(\textbf{x})=\bigg{(}\frac{1}{x_{1}}+\frac{1}{x_{2}}\bigg{)}A(\boldsymbol{\omega}),$
(3)

where $\mathcal{S}_{2}$ is the unit simplex and
$\boldsymbol{\omega}=(\omega_{1},\omega_{2})=(\frac{x_{1}}{\|\textbf{x}\|},\frac{x_{2}}{\|\textbf{x}\|})\in\mathcal{S}_{2}$,
for any norm $\|\cdot\|$. The Pickands dependence function is the restriction
of the exponent function $V$ to $\mathcal{S}_{2}$, and it measures the
strength of extremal dependence along a direction
$\boldsymbol{\omega}\in\mathcal{S}_{2}$. This function must satisfy the
boundary restrictions $\max{(\omega_{1},\omega_{2})}\leq
A(\boldsymbol{\omega})\leq 1$ where the lower bound corresponds to perfect
extremal dependence and the upper bound to the independence case. It also
satisfies the border conditions $A(0,1)=A(1,0)=1$ and is a continuous and
convex function on $S_{2}$. The extremal coefficient (Smith, 1990), defined by

$\theta=V(1,1)=2A(1/2,1/2),$ (4)

is a summary of the strength of extremal dependence between $Z_{1}$ and
$Z_{2}$. Noting that $P(\textbf{Z}\leq z)=\exp(-1/z)^{\theta}$, the extremal
coefficient lies between 1 (perfect dependence) and 2 (independence).
Moreover,

$\theta=2-\lim_{u\to\infty}P(Z_{1}>u\mid Z_{2}>u).$ (5)

In the rest of the paper, we focus mainly on modelling $\theta$.

For modelling $G$, it is common to use parametric representations of the
exponent function $V$. An example is the logistic model, defined as

$P(Z_{1}\leq z_{1},Z_{2}\leq
z_{2})=\exp\\{-\big{(}z_{1}^{-1/\alpha}+z_{2}^{-1/\alpha}\big{)}^{\alpha}\\},\alpha\in(0,1].$
(6)

Independence and complete dependence correspond to $\alpha=1$, $\alpha\to 0$,
respectively (Tawn, 1988). It is a symmetric model, where $Z_{1},Z_{2}$ are
exchangeable and a simple relation $\theta=2^{\alpha}$ holds. Another example
is the Hüsler–Reiss model, defined as

$P(Z_{1}\leq z_{1},Z_{2}\leq
z_{2})=\exp\bigg{\\{}-\Phi\big{(}\lambda+\frac{z_{2}-z_{1}}{2\lambda}\big{)}e^{-z_{1}}-\Phi\big{(}\lambda+\frac{z_{1}-z_{2}}{2\lambda}\big{)}e^{-z_{2}}\bigg{\\}},$
(7)

where $\Phi$ is a Gaussian distribution function and $\lambda\in(0,\infty)$ is
a parameter. Independence and complete dependence correspond to
$\lambda\to\infty$, $\lambda\to 0$, respectively, and $\theta=2\Phi(\lambda)$
holds (Engelke et al., 2015). Other commonly used parametric models include,
for example, Dirichlet models (Tawn, 1990).

In many applications, the extremal dependence structure depends on covariates.
Consider the bivariate max-stable random vector
$\textbf{Z}(\textbf{x})=(Z_{1}(\textbf{x}),Z_{2}(\textbf{x}))^{\top}$ given a
set of $d$ observed covariates $\textbf{x}=(x_{1},\ldots,x_{d})^{\top}$ and
suppose that it has unit Fréchet marginal distributions. Both the related
Pickands function and extremal coefficient are therefore covariate-dependent
and can be written $A(\boldsymbol{\omega}|\textbf{x})$ for
$\boldsymbol{\omega}\in S_{2}$ and
$\theta(\textbf{x})=2A\left\\{(\frac{1}{2},\frac{1}{2})\mid\textbf{x}\right\\}$.
While our emphasis is on non-stationary extremal dependence structures, one
may also consider the margins of $\textbf{Z}(\textbf{x})$ to be covariate-
dependent, that is, $F_{j,\textbf{x}}(z)=P(Z_{j}(\textbf{x})\leq z)$ for
$j=1,2$. In that case, a regression model applied to the marginal GEV
distributions’ parameters (Pauli and Coles, 2001) is required to transform the
margins to a unit Fréchet scale.

In order to infer about the covariate-dependent Pickands function, Mhalla et
al. (2017) introduced, for a fixed direction in the unit simplex, the max-
projection

$Y(\omega,1-\omega)\mid\textbf{x}=\max[H^{\lambda_{1}}\\{Z_{1}(\textbf{x})\\},H^{\lambda_{2}}\\{Z_{2}(\textbf{x})\\}],\quad\omega\in(0,1),$

where
$\lambda_{1}=\omega(\frac{1}{\omega}+\frac{1}{1-\omega}),\lambda_{2}=(1-\omega)(\frac{1}{\omega}+\frac{1}{1-\omega})$,
and $H(x)=e^{-1/x}$. It was shown that the resulting random variable is Beta
distributed with parameters
$A\left\\{(\omega,1-\omega)\mid\textbf{x}\right\\}$ and $1$. Hence, for a
fixed value of $\omega$, statistical inference for
$A\left\\{(\omega,1-\omega)\mid\textbf{x}\right\\}$ can be performed through
likelihood procedures or regression techniques.

##  2 Causal modelling

First, we introduce a graphical causal notation (e.g., Spirtes et al. (2001)).
A DAG (directed acyclic graph) $\mathcal{G}=(V,E)$ contains a finite set of
vertices (nodes) $V$ and a set of directed edges $E$ between distinct
vertices. We assume that there exist no directed cycle or multiple edges. Let
$i,j\in V$ be two distinct nodes of $\mathcal{G}$. We say that $i$ is a parent
of $j$ if there exist an edge from $i$ to $j$ in $E$, and denote $i\in
pa_{j}(\mathcal{G})$. We usually omit $\mathcal{G}$ in the notation if it is
clear from context. Consider a random vector $(X_{i})_{i\in V}$ over a
probability space $(\Omega,\mathcal{A},P)$. With a slight abuse of notation,
we identify the vertices $j\in V$ with the variables $X_{j}$. We denote
$\textbf{X}_{S}=\\{X_{s}:s\in S\\}$ for $S\subseteq V$. Typically, we denote
$V=\\{0,\dots,d\\}$, where $X_{0}$ (usually denoted by $Y$) is a variable of
interest and $\textbf{X}=(X_{1},\dots,X_{d})^{\top}$ are covariates. We assume
that $(Y,\textbf{X})$ follow a structural causal model (SCM) (Pearl, 2009),
defined as follows.

######  Definition 1.

The random variables $X_{0},\dots,X_{d}$ (where $X_{0}$ plays the role of the
target variable $Y$) follow a SCM with DAG $\cal{G}$ representing the causal
structure (Verma and Pearl, 2013) if for each $i\leq d$ the variable $X_{i}$
arises from a structural equation

$X_{i}=g_{i}(\textbf{X}_{pa_{i}},\varepsilon_{i}),$

where $\varepsilon_{i}$ are jointly independent random variables. The
functions $g_{i}$ are the assignments (or link functions) and we say that
$X_{j}$ is a direct cause of $X_{i}$, if $j\in pa_{i}(\mathcal{G})$.

We are interested in estimating the set of direct causes of $Y$. Doing that
only from a random sample is often impossible, unless we assume some
identifiable structure and additional strong assumptions (Shimizu et al.,
2006; Peters et al., 2014).

######  Definition 2 (Intervention).

Consider an SCM as in Definition 1. We replace one of the structural equations
to obtain a new SCM. Assume that we replace the assignment for $X_{i}$ by
$X_{i}=\tilde{g}_{i}(\textbf{X}_{\tilde{pa}_{i}},\tilde{\varepsilon}_{i})$,
where the new $\tilde{\varepsilon_{i}}$ is independent of all
$\varepsilon_{j}$ and $\tilde{g}_{i},\tilde{pa_{i}}$ are potentially
arbitrary. We then call the entailed distribution of the new SCM an
intervention distribution and say that $X_{i}$ has been intervened on.

A classical example of an intervention is a hard (atomic) intervention on
$X_{i}$ (Peters et al., 2017, Definition 6.8.), replacing the structural
equation in the SCM by $X_{i}=const$. Note the difference between a simple
conditioning $Y\mid X_{i}=const$ and a hard intervention, where we assume that
the joint distribution $P_{(Y^{e},\textbf{X}^{e})}$ results from a different
data-generating process.

More generally, suppose that we observe $(Y,\textbf{X})$ in different
environments (Peters et al., 2017, Section 7.1.6). In each environment
$e\in\mathcal{E}$, there is a target variable $Y^{e}$ and covariates
$\textbf{X}^{e}=(X_{1}^{e},\dots,X_{d}^{e})^{\top}$. The number of
environments $|\mathcal{E}|$ is assumed finite. The set of environments
$\mathcal{E}$ is defined in a general way by assuming that for each
$e\in\mathcal{E}$, $(Y^{e},\textbf{X}^{e})$ follows (a different) SCM with
distribution $P^{e}=P_{(Y^{e},\textbf{X}^{e})}$. A different environment can
be generated by (possibly several) interventions (Eberhardt and Scheines,
2007) on the covariates. However, we do not assume that any information about
the interventions is known.

Note that interventions on X do not change the conditional distribution of
$Y\mid\textbf{X}$. We assume that the parents of $Y^{e}$ remain the same in
all environments. Consider the set $S^{\star}\subseteq\\{1,\dots,d\\}$ of
indices determining the covariates of X that are true direct causes of $Y$.
Invariance based causal discovery (Peters et al., 2016) relies on the fact
that $Y^{e}\mid\textbf{X}^{e}_{S^{\star}}$ is invariant for each
$e\in\mathcal{E}$. We can put it more formally by assuming the following:

######  Assumption 1.

There exist $S^{\star}\subseteq\\{1,\dots,d\\}$, a function
$g:\mathbb{R}^{|S^{\star}|+1}\to\mathbb{R}$ and a distribution
$F_{\varepsilon}$, such that for all $e\in\mathcal{E}$ we have

$Y^{e}=g(\textbf{X}^{e}_{S^{\star}},\varepsilon^{e}),\,\,\varepsilon^{e}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\textbf{X}^{e}_{S^{\star}},\,\,\varepsilon^{e}\sim
F_{\varepsilon}.$ (8)

In other words, the structural equation for the target variable can be written
as $Y^{e}=g(\textbf{X}^{e}_{S^{\star}},\varepsilon^{e})$ for all
$e\in\mathcal{E}$. The link function $g$ remains the same in all environments,
and noise variables $\varepsilon^{e}$ are iid. As discussed in Bühlmann
(2020), Assumption 1 guarantees that on one hand $e$ does not change the
mechanism between $\textbf{X}^{e}$ and $Y^{e}$, and on the other hand that $e$
does not act directly on $Y$ and ideally, $e$ should change the distribution
of $\textbf{X}^{e}$.

In what follows, the formulation “Assumption 1 is valid for set $S$”,
$S\subseteq\\{1,\dots,d\\}$ means that there exist $g$ and $F_{\varepsilon}$
such that (8) holds with set $S$. It is important to observe the following
(Peters et al., 2016, Section 6.1):

######  Lemma 1.

Assumption 1 is valid for the set S if and only if for all $e,f\in\mathcal{E}$
we have
$Y^{e}\mid(\textbf{X}^{e}_{S}=\textbf{x})\overset{D}{=}Y^{f}\mid(\textbf{X}^{f}_{S}=\textbf{x})$
for all x.

The ICP methodology stems from testing for each $S\subseteq\\{1,\dots,d\\}$
the following hypothesis:

$H_{0,S}(\mathcal{E}):\text{Assumption 1 is valid for the set S}.$

Therefore,

$\hat{S}^{\star}=\bigcap_{S:H_{0,S}(\mathcal{E})\text{ is not rejected}}S$ (9)

is the proposed estimator of the set of direct causes $S^{\star}$.

If Assumption 1 holds and we have a proper test for $H_{0,S}(\mathcal{E})$,
then $H_{0,S^{\star}}$ will not be rejected on a given confidence level, and
therefore $\hat{S}^{\star}\subseteq S^{\star}$. On the other hand, if we want
to guarantee $\hat{S}^{\star}=S^{\star}$, we need to observe “enough”
environments. It is a big challenge to show some necessary and sufficient
conditions for the full identifiability of the causal variables. Some results
are known in only special cases (Pfister and Peters, 2022). It is easy to see
that the more environments we observe, the larger $\hat{S}^{\star}$ is.

As for testing $H_{0,S}(\mathcal{E})$, we need to assume some form of the link
function $g$. Usually, classical regression methods are used to infer about
the link function $g$ and testing $H_{0,S}(\mathcal{E})$ is performed through
a testing procedure such as a $k$-sample test. However, if we have a test
$H_{0,S}(\mathcal{E})$ having a proper level $\alpha$, the Theorem 1
guarantees that each element of $\hat{S}^{\star}$ indeed lies in $S^{\star}$
with probability $1-\alpha$.

######  Theorem 1.

(Peters et al., 2016) Let Assumption 1 be valid. Consider a test
$H_{0,S}(\mathcal{E})$ that has a level $\alpha$ for
$S\subseteq\\{1,\dots,d\\}$. Then, the previously defined estimate
$\hat{S}^{\star}$ fulfills

$P(\hat{S}^{\star}\subseteq S^{\star})\geq 1-\alpha.$

Hence, a proper test $H_{0,S}(\mathcal{E})$ that has a level $\alpha$ (at
least asymptotically) is needed. One of the goals of the work by Heinze-Deml
et al. (2018) was to conduct an appropriate $H_{0,S}(\mathcal{E})$ test under
different assumptions on the link function. Here, we use the invariant
residual distribution test explained below.

###  2.1 Invariant residual distribution test

We assume that the link function $g$ from Assumption 1 is additive, that is,
the structural equation is of the form

$Y^{e}=g(\textbf{X}^{e}_{S^{\star}})+\varepsilon^{e},\,\,\varepsilon^{e}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\textbf{X}^{e}_{S^{\star}},\,\,\varepsilon^{e}\sim
F_{\varepsilon}$

for all $e\in\mathcal{E}$. The core of the proposed test for
$H_{0,S}(\mathcal{E})$ has two main steps:

  1. 1.

Pool the data together from all environments and fit a model to predict $Y$
from $\textbf{X}_{S}$. Stated differently, estimate $g$ relying on data from
all environments.

  2. 2.

Use a $k$-sample test to assess whether the residuals obtained from different
environments have the same distribution.

A more detailed description can be found in Appendix II.5 in Heinze-Deml et
al. (2018).

##  3 Causal modelling of extremal dependence

In order to frame the modelling of extremes into a structural causal setup, we
start by introducing the LogMax-projection.

######  Lemma 2.

Assume that $Z=(Z_{1},Z_{2})^{\top}$ is a bivariate simple max-stable random
vector with a Pickands dependence function $A_{Z_{1},Z_{2}}(\cdot)$ and
extremal coefficient $\theta$. Denote
$A:=A_{Z_{1},Z_{2}}(\frac{1}{2},\frac{1}{2})$. Let

$\displaystyle Y^{\prime}:=$
$\displaystyle\min\bigg{(}\frac{2}{Z_{1}},\frac{2}{Z_{2}}\bigg{)},$
$\displaystyle Y:=$ $\displaystyle\log\\{\max(Z_{1},Z_{2})\\}-\gamma,$

where $\gamma\approx 0.57$ is the Euler constant. Then,

  * •

$Y^{\prime}$ is exponentially distributed with parameter $A$. In particular,
$\mathbb{E}(Y^{\prime})=\frac{1}{A}$.

  * •

$Y\overset{D}{=}\log(\theta)+\varepsilon$, where $\varepsilon\sim
Gumbel(-\gamma,1)$222$Gumbel(\mu,\beta)$ distribution function has the form
$P(\varepsilon\leq x)=e^{-e^{-(x-\mu)/\beta}}$. Specifically,
$Gumbel(-\gamma,1$) is a centered Gumbel distribution.. In particular,
$\mathbb{E}(Y)=\log(\theta)$.

###### Proof.

From $P\\{\max(Z_{1},Z_{2})\leq x\\}=e^{-\frac{\theta}{x}}$ and $\theta=2A$,
we get

$e^{-Ax}=P\\{\max(Z_{1},Z_{2})\leq\frac{2}{x}\\}=P\\{\min\bigg{(}\frac{2}{Z_{1}},\frac{2}{Z_{2}}\bigg{)}\leq
x\\}.$

Similar results can be found in Mhalla et al. (2017) using the beta
distribution. The results for $Y$ follow easily when connecting
$Y=-\log\big{(}\frac{Y^{\prime}}{2}\big{)}-\gamma$. We use the fact that the
minus logarithm of an exponential distribution is a Gumbel distribution and
$\theta=2A$. The mean of the standard Gumbel distribution is $\gamma$, so we
obtain centred Gumbel noise by subtracting $\gamma$. ∎

We will refer to $Y$ defined in Lemma 2 as a LogMax-projection. This
transformation is very useful when the interest is to make inference on
$\theta$. This can be done for instance through classical regression tools.
More importantly, it has the desired property that the noise acts additively.
Hence, it can be fitted into the causal framework of Heinze-Deml et al.
(2018).

######  Remark 1.

Lemma 2 can be used to model the behaviour of the Pickands dependence function
evaluated at the specific value $(\frac{1}{2},\frac{1}{2})$. Arguably, this
value is the most important in most applications due to its interpretability
and connections to the tail dependence measures such as (5). However, we can
easily modify Lemma 2 for any value $\boldsymbol{\omega}=(\omega,1-\omega)$,
$\omega\in(0,1)$, by taking $Y=\log\\{\max((1-\omega)Z_{1},\omega
Z_{2})\\}-\gamma$.

###  3.1 Covariate-dependent structures and GAM

Let us move now to the case where the extremal coefficient $\theta$ depends on
a set of covariates $\textbf{X}\in\mathcal{X}\subseteq\mathbb{R}^{d}$.
Consider a bivariate simple max-stable random vector
$Z(\textbf{x})=(Z_{1}(\textbf{x}),Z_{2}(\textbf{x}))^{\top}$with a covariate-
dependent extremal coefficient. For any $\textbf{x}\in\mathcal{X}$, we denote
by $\theta(\textbf{x})$, the extremal coefficient between $Z_{1},Z_{2}$
conditioned on $\textbf{X}=\textbf{x}$.

From Lemma 2, we have that the LogMax-projection is defined as

$Y(\textbf{x}):=\log\\{\max(Z_{1}(\textbf{x}),Z_{2}(\textbf{x}))\\}-\gamma,$

and satisfies

$Y(\textbf{x})\overset{D}{=}\log\\{\theta(\textbf{x})\\}+\varepsilon,\text{
where }\varepsilon\sim Gumbel(-\gamma,1).$

As mentioned above, one can use this model formulation to make inference on
$\theta(\textbf{x})$. We propose using the flexible generalized additive model
(GAM) (Wood, 2017) for the estimation of $\theta(\textbf{x})$. Under the GAM
framework, we make the following assumptions:

######  Assumption 2.

Let
$\theta(\textbf{x})=h\\{\textbf{u}^{\top}\boldsymbol{\beta}+\sum_{k=1}^{K}h_{k}(t_{k})\\}$,
where:

  * •

$h:\mathbb{R}\to[1,2]$ is a link function,

  * •

$(u_{1},\dots,u_{s})$ and $(t_{1},\dots,t_{K})$ are subsets of x, and

  * •

$h_{k}:\mathbb{T}_{k}\to\mathbb{R}$ are functions admitting a finite quadratic
penalty representation (Green and Silverman (1994)) supported on a closed set
$\mathbb{T}_{k}\subset\mathbb{R}$ for all $k$.

Note that in application, ignoring the fact that the extremal coefficient lies
in $[1,2]$ and considering that it is real-valued instead, has no impact on
the set of causal covariates discovery. For instance, we do not aim to
estimate the extremal coefficient, per se, but rather compare the effect of
the covariates under different environments. This does not affect the
consistency. The estimation of $\theta(\textbf{x})$ is based on maximizing the
penalized log likelihood:

$l(\Gamma,\gamma)=l(\Gamma)-\sum_{k=1}^{K}\gamma_{k}\int_{\mathbb{T}_{k}}h_{k}^{\prime\prime}(t_{k})^{2}dt_{k},$
(10)

where $\Gamma$ represents the set of unknown parameters ($\beta$ and the
linear basis coefficients of each of the smooth functions $h_{k}$) and
$l(\Gamma)$ is the log likelihood associated with the $n$ Gumbel observations
$\\{(y_{i},\textbf{x}_{i})_{i=1}^{n}\\}$, that is,
$l(\Gamma)=\sum_{i=1}^{n}\theta(\textbf{x}_{i})-y_{i}-e^{\theta(\textbf{x}_{i})-y_{i}}$.
The integrated square second derivative penalties prevent interpolating
estimate functions of $h_{k}$ and the smoothing parameters $\gamma_{k}$
control the trade-off between model goodness of fit and model smoothness. When
$\gamma_{k}\to\infty$ the estimate of $h_{k}(x_{k})$ tends to be a straight
line.

The minimization of (10) is performed based on an outer-iteration procedure.
At each iteration, $\Gamma$ and $\gamma$ are estimated separately by penalized
iteratively re-weighted least squares and a prediction error method from
generalized cross validation, respectively (for details, see Chapter 5 in Wood
(2017)). Under appropriate regularity assumptions on the functions $h_{k}$ and
the smoothing parameters, the GAM resulting estimator
$\hat{\theta}(\textbf{x})=h\\{\textbf{u}^{\top}\hat{\boldsymbol{\beta}}+\sum_{k=1}^{K}\hat{h}_{k}(t_{k})\\}$
is consistent and asymptotically normal (Yoshida and Naito, 2012; Wood, 2017).
More GAM-related theory is available, for example, in Wood et al. (2016).

###  3.2 LogMax-projection in SCM

Let us move to the causal setup. We assume that we observe variables from
different environments $\mathcal{E}$, with the number of environments
$|\mathcal{E}|$ being finite. In environment $e\in\mathcal{E}$, let
$\textbf{X}^{e}=(X_{1}^{e},\dots,X_{d}^{e})^{\top}$ be our covariates,
$\textbf{Z}^{e}=(Z_{1}^{e},Z_{2}^{e})$ be a simple max-stable random vector
with covariate-dependent extremal coefficient $\theta(\textbf{X}^{e})$. Let
$Y^{e}$ be the corresponding LogMax-projection in environment
$e\in\mathcal{E}$. In each environment, we observe $n_{e}\in\mathbb{N}$ iid
observations $(\textbf{X}^{e}_{i},\textbf{Z}^{e}_{i})_{i=1}^{n_{e}}$. For
$S\subseteq\\{1,\dots,d\\}$, we denote $\theta(\textbf{X}^{e}_{S})$ the
extremal coefficient between $Z_{1},Z_{2}$ conditioned on
$\textbf{X}^{e}_{S}=\textbf{x}^{e}_{S}$ (as a reminder, we denote
$\textbf{X}_{S}=(X_{i_{1}},\dots,X_{i_{|S|}})$ for a set
$S=\\{i_{1},\dots,i_{|S|}\\}\subseteq\\{1,\dots,d\\}$).

Combining Assumption 1 and Lemma 2, we assume that the structural equation is
of the form

$Y^{e}=\log\\{\theta(\textbf{X}^{e}_{S^{\star}})\\}+\varepsilon^{e},\varepsilon^{e}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\textbf{X}^{e}_{S^{\star}},\,\,\varepsilon^{e}\sim
Gumbel(-\gamma,1),$

where $S^{\star}$ is the true set of parents. The following proposition states
that the residuals are asymptotically invariant.

######  Proposition 1.

Suppose Assumption 1 holds. Let $\hat{\theta}$ be an estimator of $\theta$
that is consistent in the sense that
$\hat{\theta}(\textbf{X}_{S^{\star}}^{e})-\theta(\textbf{X}_{S^{\star}}^{e})\overset{P}{\to}0$
for all $e\in\mathcal{E}$. Define the residuals
$\hat{\eta}_{S^{\star}}^{e}:=Y^{e}-\log\\{\hat{\theta}(\textbf{X}_{S^{\star}}^{e})\\}$.
Take $f\neq e\in\mathcal{E}$. Then, as the number of data (in both
environments) tends to infinity, we have
$\hat{\eta}^{e}_{S^{\star}}\overset{D}{=}\hat{\eta}^{f}_{S^{\star}}$.

The proof is based on the fact that
$Y^{e}-\log\\{\theta(\textbf{X}_{S^{\star}}^{e})\\}\overset{D}{=}Y^{f}-\log\\{\theta(\textbf{X}_{S^{\star}}^{f})\\}$
from Assumption 1 and that the estimator is consistent. Hence, the invariant
residual distribution test (2.1) can be directly used. Note that Proposition 1
is formulated only for the set of true direct causes $S^{\star}$.

We can directly test the invariance property
$\hat{\eta}^{e}_{S}\overset{D}{=}\hat{\eta}^{f}_{S}$ for all
$S\subseteq\\{1,\dots,d\\}$. This invariance should hold for the true set
$S^{\star}$ and hence, the same argument as in Theorem 1 can be used. In the
following, we propose a methodology for testing this invariance property.

###  3.3 Generic algorithm

The proposed methodology can be summarized in the following algorithm.

  1. 1.

For each margin separately:

     * •

Transform original data to unit Fréchet scale. This can be done by using the
block maxima approach to fit a GEV distribution to the univariate series.

     * •

The GEV parameters can also depend on other covariates.

  2. 2.

In each environment $e\in\mathcal{E}$ separately:

     * •

Compute the LogMax-projection $Y^{e}$,

  3. 3.

For each $S\subseteq\\{1,\dots,d\\}$:

     * •

Pool the data from all environments and calculate
$\hat{\theta}(\mathbf{X}_{S})$,

     * •

Compute the residuals
$\hat{\eta}_{S}^{e}=Y^{e}-\log\\{\hat{\theta}(\mathbf{X}^{e}_{S})\\}$,

     * •

Use a $k$-sample test, to test if the residuals from different environments
$\hat{\eta}_{S}^{e}$ are equally distributed.

  4. 4.

Get $\hat{S}^{\star}$ as the intersection of the subsets
$S\subseteq\\{1,\dots,d\\}$ where the $k$-sample test was not rejected.

In practice, we do not have to repeat the algorithm for all
$S\subseteq\\{1,\dots,d\\}$; if we do not reject the $k$-sample test for a set
$S$, testing for a superset of $S$, is redundant.

Several approaches can be used for transforming the data to unit Fréchet
margins. In our application, we fit a covariate-dependent GEV distribution to
block maxima. Alternatively, a peak over threshold approach (Pickands, 1975)
can be used as well.

In our application, we use GAM to estimate $\theta(\textbf{X}_{S})$ based on
minimizing (10). This choice was also used in Mhalla et al. (2017) and is a
default method in non-linear ICP methodology in Heinze-Deml et al. (2018).

As for the $k$-sample test, we use the Anderson-Darling $k$-sample test (AD-
test). In dimension two, the distance between the two samples is

$AD_{n,m}=\frac{mn}{n+m}\int_{-\infty}^{\infty}\frac{(F_{n}(x)-G_{m}(x))^{2}}{H_{n+m}(x)(1-H_{n+m}(x))}dH_{n+m}(x),$

where $H_{n+m}$ is an empirical distribution function of the pooled datasets.
The $k$-sample AD-test is the natural generalization of its two-sample
counterpart (Scholz and Stephens, 1987). It has typically better power than
the Kolmogorov-Smirnov test (Razali and Yap, 2011). This test is used at a
level $\alpha=0.05$ level. Note that an estimation error can arise from two
sources: from the estimation of the margins, and from the estimation of
$\hat{\theta}$. Although the equality
$\hat{\eta}_{S^{\star}}^{e}\overset{D}{=}\hat{\eta}_{S^{\star}}^{f}$ should
asymptotically hold, we do not claim that a proper $\alpha$ level holds for a
finite number of observations. In the following, we illustrate the methodology
on a series of simulations.

##  4 Simulations

We proceed with a simulation study. Section 4.1 describes the study and
results based on a specific example, discussing power with different causal
strength, intervention strength and multicollinearity. Section 4.2 discusses a
hidden confounder case, when Assumption 1 is not fulfilled. Section 4.3 deals
with a case when we start with non-extreme data; we compare the choices of
different thresholds and different marginal transformations and we discuss a
case when we do not “clean the margins” properly.

We want to stress that we are only interested in inferring causal relations,
not an inference on the extremal dependence structure itself.

###  4.1 Intervention, causal strength, and multicollinearity

Let $\textbf{Z}=(Z_{1},Z_{2})^{\top}$ have a simple bivariate max-stable
distribution and $\mathbf{X}=(X_{1},X_{2},X_{3})^{\top}$ be a set of
covariates. We assume the distribution of $\mathbf{Z}$ to be parameterized by
a logistic copula (6) with parameter $\alpha=\alpha(X_{1})=\frac{e^{\beta
X_{1}}}{1+e^{\beta X_{1}}}$. Since $\alpha$ is a function of $X_{1}$ only, the
true set of direct causes of our extremal dependence is therefore
$S^{\star}=\\{1\\}$. We consider $X_{2},X_{3}$ as other covariates,
graphically represented by a SCM in the form $X_{2}\leftarrow X_{1}\to Y\to
X_{3}$, where $Y=\log\\{\max(Z_{1},Z_{2})\\}-\gamma$ is the LogMax-projection.
We observe the covariates $X_{1},X_{2},X_{3}$ and the change of environments
are interventions on $X_{1},X_{2}$, and $X_{3}$. Details about the generation
of the random variables can be found in Table 1 and the three parameters
$p,q$, and $\beta$ are as follows

  * •

the parameter $p$ describes the strength of intervention. If $p=0$ then there
is no intervention on $X_{1}$,

  * •

the parameter $\beta$ represents the causal strength of $X_{1}$ on $Y$. If
$\beta=0$ then $Y$ is independent of $X_{1}$,

  * •

the parameter $q$ describes the dependence between $X_{1}$ and $X_{2}$. If
$q=0$ then $cor(X_{1},X_{2})=0$. Moreover, if $q=0.25$, the correlation
$cor(X_{1},X_{2})\approx 0.25$ and if $q=1$ then $cor(X_{1},X_{2})\approx
0.75$.

We proceed with our causal methodology using the algorithm described in
Section 3.3. We use $50$ repetitions for different choices of $q,p,\beta$ and
for two sample sizes $n=500$ and $n=1500$. Figure 2 describes how often we get
$\hat{S}^{\star}=S^{\star}=\\{1\\}$ using our methodology. The upper axis in
Figure 2 represents the expectation $\mathbb{E}(\alpha)$ in the second
environment under the corresponding intervention strength $p$ (note that
$\mathbb{E}(\alpha)=1/2$ in the first environment). Obviously, the larger the
intervention $p$, the higher the heterogeneity between the environments and
hence the easier it is to find the correct causal set. The reason why we
obtain worse results with a larger absolute value $q$ is that it becomes
statistically harder to distinguish between correlated $X_{1},X_{2}$ since
both are good predictors (with an exception in the third environment). Not
surprisingly, the larger the causal strength $\beta$, the better the results
we obtain.

Table 1: Definitions of environments for the simulation study in Section 4.1 with graphical representation $X_{2}\leftarrow X_{1}\to Y\to X_{3}$. The second and the third environments consist of interventions on $X_{1},X_{2}$. The fourth environment changes $X_{3}$ by adding an additional relation $X_{1}\to X_{3}$. Note that Assumption 1 is valid since $Y\mid X_{1}$ is the same in all environments.  First environment | Second environment
---|---
$X_{1}=\varepsilon_{1}$, where $\varepsilon_{1}\sim N(0,1)$ |  $X_{1}={\color[rgb]{1,0,0}p}+\varepsilon_{1}$, where $\varepsilon_{1}\sim N(0,1)$
$X_{2}=qX_{1}+\nu,$ where $\nu\sim Exp(1)$ |  $X_{2}=qX_{1}+\nu,$ where $\nu\sim Exp(1)$
$\alpha=\alpha(X_{1})=\frac{e^{\beta X_{1}}}{1+e^{\beta X_{1}}}$ | $\alpha=\alpha(X_{1})=\frac{e^{\beta X_{1}}}{1+e^{\beta X_{1}}}$
$X_{3}=Y+\varepsilon_{3}$, where $\varepsilon_{3}\sim N(0,1)$ |  $X_{3}=Y+\varepsilon_{3}$, where $\varepsilon_{3}\sim N(0,1)$
Third environment | Fourth environment
$X_{1}=\varepsilon_{1}$, where $\varepsilon_{1}\sim N(0,1)$ |  $X_{1}=\varepsilon_{1}$, where $\varepsilon_{1}\sim N(0,1)$
$X_{2}={\color[rgb]{1,0,0}-10}$ |  $X_{2}=qX_{1}+\nu,$ where $\nu\sim Exp(1)$
$\alpha=\alpha(X_{1})=\frac{e^{\beta X_{1}}}{1+e^{\beta X_{1}}}$ | $\alpha=\alpha(X_{1})=\frac{e^{\beta X_{1}}}{1+e^{\beta X_{1}}}$
$X_{3}=Y+\varepsilon_{3}$, where $\varepsilon_{3}\sim N(0,1)$ |  $X_{3}=Y+{\color[rgb]{1,0,0}X_{1}}+\varepsilon_{3}$, where $\varepsilon_{3}\sim N(0,1)$

Figure 2: Ratio of correct results for the different choices of parameters
$p,\beta,q$ as described in Table 1 and for two sample sizes $n=500,1500$. The
upper axis represents the expectation $\mathbb{E}(\alpha)$ in the second
environment.

###  4.2 Hidden confounder

We present a situation in which a hidden confounder of a causal covariate and
the extremal dependence coefficient are present. The graphical representation
of the model can be found in Figure 3.

Figure 3: Graphical representation of causal relations in the simulation study
with a hidden confounder $H$. The green arrow from the environment represents
the interventions described in subsection 4.2.

Let $\textbf{Z}=(Z_{1},Z_{2})^{\top}$ have a simple bivariate max-stable
distribution, $X_{1},X_{2}$ are the observed covariates and $H$ is an
unobserved covariate. We assume the distribution of $\mathbf{Z}$ to be
parameterized by a Hüsler–Reiss copula (7) with parameter
$\lambda=\lambda(X_{2},H)=|X_{2}+3H|/5$.

We intervene on $X_{1}$ and generate the variables in the corresponding
environments:

  1. 1.

$H+1\sim Exp(1)$, $X_{1}=\varepsilon_{1},X_{2}=X_{1}+H+\varepsilon_{2}$ where
$\varepsilon_{1},\varepsilon_{2}\overset{iid}{\sim}N(0,1)$,

  2. 2.

$H+1\sim Exp(1)$, $X_{1}=\varepsilon_{1}+\color[rgb]{1,0,0}{5}$,
$X_{2}=X_{1}+H+\varepsilon_{2}$ where
$\varepsilon_{1},\varepsilon_{2}\overset{iid}{\sim}N(0,1)$.

Moreover, in both environments, we compute the LogMax-projection $Y$.

Finally, we compute $\hat{S}^{\star}$ based on $X_{1},X_{2},Y$ only. Note that
the only observed direct cause of $Y$ is $X_{2}$. However, due to the hidden
variable $H$, Assumption 1 no longer holds. It is easy to see that $Y\mid
X_{1}$ is invariant through environments, however, $Y\mid X_{2}$ is not due to
the hidden confounder. Simulations show how much this heuristic is justified
in practice.

For a range of sample sizes $n$, we generate the variables as described above,
compute $\hat{S}^{\star}$ and repeat $50$ times. The results are described in
Figure 4. We can see that for smaller $n$, we tend to not reject $H_{0,S}$ for
both $S=\\{2\\},S=\\{1\\}$. Hence, the resulting estimate is an empty set
(intersection of $\\{1\\}\cap\\{2\\}$). For a larger $n$, the test can
recognise the effect of the hidden confounder and we tend to reject $H_{0,S}$
for $S=\\{2\\}$, and hence, the resulting estimate is
$\hat{S^{\star}}=\\{1\\}$.

Note that if the effect of $H$ is small, we tend to accept $H_{0,S}$ for
$S=\\{2\\}$ and also for $S=\\{1\\}$ which results in
$\hat{S^{\star}}=\emptyset$. On the other hand, if the effect of $H$ is large
(larger than the effect of intervention), then the effect of $X_{1}$ on $Y$
will be (proportionally) small and we tend to accept $H_{0,S}$ for
$S=\emptyset$. The statistical guarantees such as Theorem 1 or Theorem 1 for
$\hat{S}^{\star}$ are no longer valid. However, following (Peters et al.,
2016)[Proposition 5], $\hat{S}^{\star}$ will estimate some subset of its
ancestors (instead of parents). For a more detailed discussion of hidden
confounders in the non-extremal case, see (Peters et al., 2016)[Appendix C].

Figure 4:  The figure corresponds to the simulation study from Section 4.2,
with a hidden confounder $H$. The figure represents the final estimate
$\hat{S}^{\star}$. For a smaller $n$, the estimate
$\hat{S}^{\star}=\emptyset$, and for a larger $n$, the estimate
$\hat{S}^{\star}=\\{1\\}$.

###  4.3 Effect on the margins

In this section, we consider the case when the observed covariate has an
effect on the margins as well as on the dependence structure. Suppose we
observe two covariates $X_{1},X_{2}$, where $X_{1}$ has an effect on the
extremal dependence of $\textbf{Z}=(Z_{1},Z_{2})$ and $X_{2}$ has an effect on
the margins $Z_{1}$ and $Z_{2}$. The graphical representation can be seen in
Figure 5.

Figure 5: Graphical representation of causal relations in the simulation study
with different marginals. Here, $C_{Z_{1},Z_{2}}$ represents the dependence
structure and $Z_{1},Z_{2}$ are marginals. Green arrows from the environment
represent the interventions described in subsection 4.3.

We simulate the following data: for $n,\tau\in\mathbb{N}$, we consider
$\tilde{\textbf{Z}}^{1},\tilde{\textbf{Z}}^{2},\dots,\tilde{\textbf{Z}}^{\tau\cdot
n}$ are bivariate iid random variables distributed as a generic element
$\tilde{\textbf{Z}}=(\tilde{Z}_{1},\tilde{Z}_{2})$, where margins
$\tilde{Z}_{1}=aX_{2}+\varepsilon_{1},\tilde{Z}_{2}=aX_{2}+\varepsilon_{2}$,
where $\varepsilon_{1},\varepsilon_{2}\overset{iid}{\sim}N(0,1)$ for a
parameter $a\geq 0$. The dependence structure between
$\tilde{Z}_{1},\tilde{Z}_{2}$ is modeled by a logistic copula with parameter
$\alpha=\alpha(X_{1})$. Hence, the dependence structure (as well as the
extremal dependence structure) is a function of $X_{1}$ only.

We transform
$\tilde{\textbf{Z}}^{1},\tilde{\textbf{Z}}^{2},\dots,\tilde{\textbf{Z}}^{\tau\cdot
n}$ into a $\tau$ period-maxima, that is, we define

$\textbf{Z}^{1}:=\max(\tilde{\textbf{Z}}^{1},\tilde{\textbf{Z}}^{2},\dots,\tilde{\textbf{Z}}^{\tau}),\textbf{Z}^{2}:=\max(\tilde{\textbf{Z}}^{\tau+1},\tilde{\textbf{Z}}^{\tau+2},\dots,\tilde{\textbf{Z}}^{2\tau}),\dots,$

where the maxima is component-wise maxima.

Hence, we obtain independent maxima
$\textbf{Z}^{i}=(Z^{i}_{1},Z^{i}_{2}),i=1,\dots,n$. We model their margins by
a GEV distribution with parameters depending on $X_{2}$. By fitting the
covariate-dependent GEV distribution, we transform the margins into unit
Frechet using (2). We model GEV parameters

$\mu=\mu_{0}+f_{1}(X_{2}),\sigma=\sigma_{0}+f_{2}(X_{2}),\xi=\xi_{0}$

and estimate them smoothly using the GAM framework. Finally, we compute the
LogMax-projections $Y^{1},\dots,Y^{{n}}$.

We intervene on the covariates and observe the system. More precisely, we
generate the variables in the corresponding environments as follows:

  1. 1.

$X_{1}=1$, $X_{2}=\varepsilon_{2}$ where
$\varepsilon_{2}\overset{iid}{\sim}Exp(1)$,

  2. 2.

$X_{1}=\color[rgb]{1,0,0}{5}$, $X_{2}=\varepsilon_{2}$ where
$\varepsilon_{2}\overset{iid}{\sim}Exp(1)$,

  3. 3.

$X_{1}=1$, $X_{2}\overset{}{=}\color[rgb]{1,0,0}{1}$.

We consider the copula parameter in the form $\alpha(X_{1})=\frac{1}{X_{1}}$.
In both environments, we generate
$\tilde{\textbf{Z}}^{1},\tilde{\textbf{Z}}^{2},\dots,\tilde{\textbf{Z}}^{\tau\cdot
n}$ and compute the LogMax-projections as mentioned before.

Finally, we estimate the set of direct causes of the extremal dependence
structure following points 3 and 4 in Section 3.3. We repeat $50$ times for a
range of parameters $a\in[0,2]$, $n=1000$ and $\tau\in\\{50,100\\}$, where $a$
represents the effect of $X_{2}$ on the margins, $n$ is the size of the final
dataset (of extremes) and $\tau$ is a size of the block maxima we take. The
blue lines in Figure 6 represent the ratio of correctly estimating
$\hat{S}^{\star}=\\{X_{1}\\}$.

Now, consider a different question: What happens if we do not include $X_{2}$
for the marginal modelling? We proceed with the same setup as before, but we
model GEV parameters $\mu=\mu_{0},\sigma=\sigma_{0},\xi=\xi_{0}$ and estimate
them using maximum likelihood.

Again, we estimate the set of direct causes of the extremal dependence
structure following points 3 and 4 in Section 3.3. We repeat $50$ times for a
range of parameters $a\in[0,2]$, $n=1000$ and $\tau\in\\{50,100\\}$. The red
lines in Figure 6 represent the ratio of correctly estimating
$\hat{S}^{\star}=\\{X_{1}\\}$.

We can see that ignoring the marginal effects of $X_{2}$ creates a spurious
causal effect on the extremal dependence structure. Indeed, if $X_{2}$ is
large, the margins of Z tend to be both very large. Hence, extremes occur
concurrently even though $X_{2}$ has no effect on the extremal copula, and our
results are compromised.

The same holds if $X_{2}$ is hidden; if there exists a hidden covariate that
affects both margins with “enough” strength, it can affect the inference about
the extremal dependence structure. Therefore, it is good practice to include
as many relevant covariates for the marginal inference as possible.

Figure 6:  Simulation study from Section 4.3, where $X_{2}$ has an effect on
the marginals with strength $a\geq 0$ and the marginal block-maxima are of
size $\tau$. The blue line correspond to the case in which we include $X_{2}$
in marginal modelling. The red line correspond to the case in which we do not
include $X_{2}$ into marginal modelling. The lines represent the ratio of
results when $\hat{S}^{\star}=\\{1\\}$. Each time when
$\hat{S}^{\star}\neq\\{1\\}$, then $\hat{S}^{\star}=\emptyset$.

##  5 NO2 Application

We illustrate the causal methodology for the tail dependence on a dataset of
$\text{NO}_{2}$ concentration levels measured in the UK (Carslaw and Ropkins,
2022). To explain the extremal coefficient $\theta$ between the
$\text{NO}_{2}$ levels of two different stations, we consider three different
potential causal covariates $\textbf{X}=(X_{1},X_{2},X_{3})^{\top}$. Namely,
the distance (in kilometres) between the stations, the time (year) of the
observation, and the type of area (traffic/non-traffic). A similar dataset is
studied in Mhalla et al. (2019c). We wonder which of these covariates have a
causal effect on the tail dependence. The choice of the time as a covariate
may sound unusual as we expect that the physical mechanism is inherently given
by laws of nature and does not change over time. However, it might be
considered as a proxy for potential confounders such as the increase in
temperature or $\text{CO}_{2}$.

###  5.1 Transforming margins and LogMax-projecting

We follow steps 1 and 2 from the algorithm (3.3). We transform our
observations into bivariate simple max-stable ones. The first step is
therefore transforming the margins to the unit Fréchet scale using
transformation (2). We use weekly maxima (around 900 weeks at each station;
recall that we have hourly measurements) and model these weekly maxima by a
GEV distribution (1) with time-varying parameters

$\displaystyle\mu(t_{1},t_{2})$
$\displaystyle=\mu_{0}+f_{1}(t_{1})+f_{2}(t_{2}),$
$\displaystyle\sigma(t_{1},t_{2})$
$\displaystyle=\sigma_{0}+g_{1}(t_{1})+f_{2}(t_{2}),$

where $t_{1}$ represents the year and $t_{2}$ represents the week of the year,
and $f_{1},f_{2},g_{1},g_{2}$ are smooth functions. Inference is performed
under the GAM framework (Youngman, 2020) and likelihood ratio tests are
conducted to assess whether we need smoothly varying terms or whether
parametric and sinusoidal terms already provide a good fit. The final fitted
model is then used to transform the data at each station to the unit Fréchet
distribution using the probability integral transform (2). Figure 7 shows a
QQ-plot of the resulting transformed observations for one of the stations
(RBKC Chelsea in London) against the unit Fréchet. It shows a reasonable fit
besides a few outliers.

The next step of the algorithm (3.3) is to pair the stations and compute the
LogMax-projection. We pair only the stations with the same traffic type. We
end up with 106 pairs of stations. For each pair, we compute the LogMax-
projection transformation described in Lemma 2.

Figure 7: QQ-plot for the transformed values against the unit Fréchet
distribution corresponding to the station Newham Wren Close in London.

###  5.2 Defining environments

In this application, we define one environment as one pair of stations.
Consider two pairs of stations for which all covariates remain unchanged, but
$X_{3}$ changes from non-traffic to traffic. This can be considered as a "do-
intervention" on the covariate $X_{3}$ since we assume that the underlying
physical mechanism is the same regardless of the location of our stations. A
few remarks on the choice of such environments are in order.

######  Remark 2.

Notice that this is not simple conditioning on traffic/non-traffic stations.
Take measurements from one pair of stations. Imagine that we conduct an
experiment where we artificially add traffic to those stations to see if the
extremal coefficient changes. Instead of adding traffic to one station, we
measure a second (different) pair that has the same value of covariates,
except for traffic. The main difference here is that we assume that the data-
generating process remains unchanged, apart from the traffic covariate.

######  Remark 3 (Possible hidden variable).

The main drawback is the following: Is there a hidden variable that differs
between the pairs and causally affects our response? This is represented in
Figure 8. We need to assume that there is no such variable or if there is, it
has a negligible impact on the response variable.

Figure 8: Possible hidden variable in the application.

######  Remark 4 (Time).

Each station is measured over a different period. Some stations started in
2003, some in 1997. For each pair, we take only the intersection of measured
weeks. Hence, we can argue that by changing the pair, we also change the
period of observations. However, this change is minor and almost negligible.
Our findings support the fact that time is not a causal covariate. This might
of course be due to changes between environments not being strong enough to
find such evidence.

###  5.3 Details for the choices of environments and results for one choice of
stations

We take five different environments, that is, five pairs of stations (10
different stations in total). We want the interventions to be “strong”. For
$X_{1}$, we choose one pair from each of the bins of distances described in
Table 2. As for $X_{3}$, we choose three pairs with the type traffic and two
with the type non-traffic. In order to make an intervention on $X_{2}$, we
would need to divide the dataset into different periods and compare them.
However, due to the limited time span of our dataset, this is not an option,
so we cannot intervene on $X_{2}$.

<10km | (10km, 20km) | (20km, 50km) | (50km, 100km) | 100km<
---|---|---|---|---
Table 2: Definition of bins of distances. We aim to have one pair of stations
for each bin.

We explain now the results on one specific choice for the environments. This
choice is described in Table 3. In Section 5.4, we repeat the experiments with
50 different choices to get more robust results. We estimate $S^{\star}$ using
the methodology described in Section 3.3. For each $S\subseteq\\{1,2,3\\}$ we
test the invariance property using step 3 in algorithm 3.3. Results (for the
specific choice of environments described in Table 3) are described in Table
4. The hypotheses corresponding to sets
$\\{1\\},\\{1,2\\},\\{1,3\\},\\{1,2,3\\}$ were not rejected, and hence,
$\hat{S}^{\star}=\\{1\\}\cap\\{1,2\\}\cap\\{1,3\\}\cap\\{1,2,3\\}=\\{1\\}$.

Pairs of sites | Distances (km) | Range of years | Type
---|---|---|---
Knightsbridge n.3- Knightsbridge n.4 | 1.37 | 2000-2018 | Traffic
Slough Colnbrook n.3-Slough Colnbrook n.4 | 11.36 | 2000-2018 | Non-traffic
Watford Town Hall - Newham | 46.61 | 2008-2018 | Traffic
Hillingdon South Ruislip - Cambridge n.4 | 94.90 | 2002-2017 | Traffic
Trafford Moss Park - Hounslow Chiswick | 316.50 | 1999-2018 | Non-traffic
Table 3: Pairs defining the five environments.  Choice for a set S | P-values
---|---
Empty set | 1e-07
1 | 0.11
2 | 0
1,2 | 0.19
3 | 6e-07
1,3 | 0.30
2,3 | 0
1,2,3 | 0.41
Intersection | 1
Table 4: Resulting p-values for one specific choice of the stations from Table
3. The p-values are computed for each subset of our covariates. It seems that
$X_{1}$ is the only identified causal covariate.

###  5.4 Repeating results with different stations and discussing results

We repeat the same steps with different choices of stations in order to obtain
some confidence bounds for each p-value. In each repetition, we randomly
choose five pairs of stations satisfying the following:

  * •

There is one pair for each bin of distances from Table 2,

  * •

Each station is included at most once,

  * •

There are three pairs from traffic type and two from non-traffic type.

Then, we compute the corresponding p-values for each subset $S$. We use 50
repetitions, and boxplots for the p-values associated with each subset of
indices are drawn in Figure 9. Again, this analysis suggests that $X_{1}$ is
the only causal covariate.

Figure 9:  Boxplots of p-values for each subset of $\\{1,2,3\\}$ when randomly
choosing 10 stations 50 times (fulfilling the constraints defined in Section
5.4).

This result is not surprising. It is expected that the extremal dependence
decreases with increasing distance. The effect of $X_{1}$ on the extremal
dependence is shown in Figure 10. The fact that $X_{2}$ was not marked as a
cause is not surprising either. We would need much stronger interventions and
changes in time observations to detect anything. Additionally, it is unclear
how time can cause extremal dependence (there may be several time-changing
factors, but we do not know them explicitly). As for $X_{3}$, it seems that
the type of the site is not a direct cause of the extremal dependence.

Figure 10:  Smooth estimation of the effect of distance on the LogMax-
projection. It is based on all available pairs of stations pooled together
(With four exceptions. If there is a triplet of stations closer than 5km from
each other, we only take one pair). The estimate may be unreliable since the
same station is used in several pairs, which creates a dependence between some
pairs. We included only the causal covariate. The estimation and the
confidence intervals are computed using the GAM framework by regressing the
LogMax-projection on the distance. We can see that the larger the distance,
the larger the extremal coefficient (and smaller the extremal dependence).

##  6 Conclusion

This paper addresses the problem of causal discovery in the area of
multivariate extremes. For a covariate-dependent extremal coefficient
$\theta(\textbf{X})$ between two variables, we discover which of the
covariates are not only significant but have a causal influence on the joint
tail behavior. We introduced a LogMax-projection, which serves as a basis
variable for inference about the extremal coefficient. A causal methodology
from Peters et al. (2016) was adapted in our context to estimate the causal
covariates. Algorithm (3.3) provides a consistent methodology. We used a very
flexible GAM modelling for statistical inference.

We stress that the results rely on the absence of a hidden confounder, an
assumption that can and should be debated in many situations. The methodology
is robust against hidden covariates to some extent. However, if there is a
hidden covariate, it must be invariant between environments. We discussed in
simulations in Section 4 the case in which this assumption is violated.

We believe that this methodology can be very useful in many fields. More
practical examples come from econometrics. We often see an extremal dependence
between stocks. If one stock has a massive fall, other stocks seem to follow
this fall too. Exploiting the causes of such behavior can be very useful and
can broadly impact the field.

Possible future work can relate to the more general multivariate extensions of
our methodology. We introduced LogMax-projection to work only with a
univariate variable representing a complex multivariate tail dependence
structure. However, it can be interesting to work directly with a tail
dependence structure of multivariate random variables with dimensions more
than $2$ and infer which covariates are significant and causal.

## Conflict of interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The work was supported by the Swiss National Science Foundation.

## References

  * Beirlant et al. (2006) J. Beirlant, Y. Goegebeur, J. Segers, J.L. Teugels, D.D.Waal, and C. Ferro.  _Statistics of Extremes: Theory and Applications_.  Wiley Series in Probability and Statistics. Wiley, 2006.
  * Bühlmann (2020) P. Bühlmann.  Invariance, causality and robustness.  _Statistical Science_ , 35(3):404 – 426, 2020\.  doi: 10.1214/19-STS721.
  * Carslaw and Ropkins (2022) D. C. Carslaw and K. Ropkins.  Openair — an R package for air quality data analysis.  _Environmental Modelling and Software_ , 2022.  ISSN 1364-8152.  doi: 10.1016/j.envsoft.2011.09.008.
  * Chang et al. (2015) R. Chang, J.R. Karr, and E.E. Schadt.  Causal inference in biology networks with integrated belief propagation.  _Pac Symp Biocomput_ , pages 359–370, 2015.  doi: 10.2307/1968974.
  * Coles (2001) S. Coles.  _An Introduction to Statistical Modeling of Extreme Values_.  Springer Series in Statistics. Springer-Verlag, 2001.  ISBN 1-85233-459-2.
  * Coles and Tawn (1991) S.G. Coles and J. A. Tawn.  Modelling extreme multivariate events.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 1991.  ISSN 00359246.  URL http://www.jstor.org/stable/2345748.
  * Dawid (2006) A. P. Dawid.  Counterfactuals, hypotheticals and potential responses: a philosophical examination of statistical causality.  _Technical report, Department of Statistical Science,University College London, Research Report_ , 2006.
  * Deuber et al. (2021) D. Deuber, J. Li, S. Engelke, and M. Maathuis.  Estimation and inference of extremal quantile treatment effects for heavy-tailed distributions.  _preprint_ , 10 2021.  doi: 10.48550/arXiv.2110.06627.
  * Eberhardt and Scheines (2007) F. Eberhardt and R. Scheines.  Interventions and causal inference.  _Philosophy of Science_ , 2007.  URL http://www.jstor.org/stable/10.1086/525638.
  * Engelke and Ivanovs (2021) S. Engelke and J. Ivanovs.  Sparse structures for multivariate extremes.  _Annual Review of Statistics and Its Application_ , 8(1):241–270, 2021.  doi: 10.1146/annurev-statistics-040620-041554.
  * Engelke et al. (2015) S. Engelke, A. Malinowski, Z. Kabluchko, and M. Schlather.  Estimation of Hüsler–Reiss distributions and Brown–Resnick processes.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 2015.  URL http://www.jstor.org/stable/24774733.
  * Fisher and Tippett (1928) R.A. Fisher and L.H.C. Tippett.  Limiting forms of the frequency distribution of the largest or smallest member of a sample.  _Mathematical Proceedings of the Cambridge Philosophical Society_ , 24(2):180–190, 1928.  doi: 10.1017/S0305004100015681.
  * Gnecco et al. (2020) N. Gnecco, N. Meinshausen, J. Peters, and S. Engelke.  Causal discovery in heavy-tailed models.  _The Annals of Statistics_ , 49, 2020.  doi: 10.1214/20-AOS2021.
  * Gnedenko (1943) B. Gnedenko.  Sur la distribution limite du terme maximum d’une serie aleatoire.  _Annals of Mathematics_ , 44(3):423–453, 1943\.  doi: 10.2307/1968974.
  * Green and Silverman (1994) P. Green and B. Silverman.  _Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach_.  Chapman and Hall/CRC, 1994.  ISBN 9780412300400.  URL https://www.routledge.com/Nonparametric-Regression-and-Generalized-Linear-Models-A-roughness-penalty/Green-Silverman/p/book/9780412300400.
  * Haavelmo (1944) T. Haavelmo.  The probability approach in econometrics.  _Econometrica_ , 1944.  URL https://www.jstor.org/stable/pdf/1906935.pdf.
  * Hamra et al. (2018) G.B. Hamra, F. Laden, A.J. Cohen, O. Raaschou-Nielsen, M.Brauer, and D. Loomis.  Invariant causal prediction for nonlinear models.  _Journal of Causal Inference_ , 6(2), 2018.  doi: 10.1515/jci-2017-0016.
  * Heinze-Deml et al. (2018) C. Heinze-Deml, J. Peters, and N. Meinshausen.  Invariant causal prediction for nonlinear models.  _Journal of Causal Inference_ , 6(2):20170016, 2018.  doi: 10.1515/jci-2017-0016.
  * Joe (1997) H. Joe.  _Multivariate Models and Multivariate Dependence Concepts, 1st Edition_ , volume 1.  Chapman and Hall/CRC, 1 edition, 1997.  URL https://doi.org/10.1201/9780367803896.
  * Kiriliouk and Naveau (2020) A. Kiriliouk and P. Naveau.  Climate extreme event attribution using multivariate peaks-over-thresholds modeling and counterfactual theory.  _The Annals of Applied Statistics_ , 2020.  doi: 10.1214/20-AOAS1355.
  * Mhalla et al. (2017) L. Mhalla, V. Chavez-Demoulin, and P. Naveau.  Non-linear models for extremal dependence.  _Journal of Multivariate Analysis_ , 159:49–66, 2017.  ISSN 0047-259X.  doi: 10.1016/j.jmva.2017.04.006.
  * Mhalla et al. (2019a) L. Mhalla, V. Chavez-Demoulin, and D.J. Dupuis.  Causal mechanism of extreme river discharges in the upper Danube basin network.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 2019a.  doi: 10.48550/ARXIV.1907.03555.
  * Mhalla et al. (2019b) L. Mhalla, M. de Carvalho, and V. Chavez-Demoulin.  Regression-type models for extremal dependence.  _Scandinavian Journal of Statistics_ , 46(4):1141–1167, 2019b.  doi: 10.1111/sjos.12388.
  * Mhalla et al. (2019c) L. Mhalla, T. Opitz, and V. Chavez-Demoulin.  Exceedance-based nonlinear regression of tail dependence.  _Extremes_ , 22, 2019c.
  * Pasche et al. (2021) O. C. Pasche, V. Chavez-Demoulin, and A.C. Davison.  Causal modelling of heavy-tailed variables and confounders with application to river flow.  _Submitted_ , 2021.  doi: 10.48550/ARXIV.2110.06686.
  * Pauli and Coles (2001) F. Pauli and S. Coles.  Penalized likelihood inference in extreme value analyses.  _Journal of Applied Statistics_ , 28(5):547–560, 2001.  doi: 10.1080/02664760120047889.
  * Pearl (2009) J. Pearl.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, 2009.  ISBN 978-0521895606.
  * Peters et al. (2014) J. Peters, J.M. Mooij, and B. Schölkopf.  Causal discovery with continuous additive noise models.  _Journal of Machine Learning Research_ , 15:2009–2053, 2014\.
  * Peters et al. (2016) J. Peters, P. Bühlmann, and N. Meinshausen.  Causal inference by using invariant prediction: identification and confidence intervals.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 78(5):947–1012, 2016.  URL https://doi.org/10.1111/rssb.12167.
  * Peters et al. (2017) J. Peters, D. Janzing, and B. Schölkopf.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  The MIT Press, 2017.  ISBN 0262037319.  URL https://library.oapen.org/bitstream/id/056a11be-ce3a-44b9-8987-a6c68fce8d9b/11283.pdf.
  * Pfister and Peters (2022) N. Pfister and J. Peters.  Identifiability of sparse causal effects using instrumental variables.  _Conference on Uncertainty in Artificial Intelligence_ , 2022.  doi: 10.48550/ARXIV.2203.09380.
  * Pickands (1975) J. Pickands.  Statistical inference using extreme order statistics.  _The Annals of Statistics_ , 3(1):119 – 131, 1975\.  doi: 10.1214/aos/1176343003.
  * Pickands (1981) J. Pickands.  Multivariate extreme value distributions.  _Bull. Int. Statist. Inst. 49_ , pages 859–78, 1981.
  * Razali and Yap (2011) N.M. Razali and B.W. Yap.  Power comparisons of shapiro-wilk, kolmogorov-smirnov, lilliefors and anderson-darling tests.  _J. Stat. Model. Analytics_ , 2, 01 2011.
  * Resnick (1987) S.I Resnick.  _Extreme Values, Regular Variation and Point Processes_.  Springer, NY, 1987.  ISBN 978-0-387-75953-1.  URL https://link.springer.com/book/10.1007/978-0-387-75953-1.
  * Scholz and Stephens (1987) F.W. Scholz and M. A. Stephens.  K-sample Anderson-Darling tests.  _Journal of the American Statistical Association_ , 82(399):918–924, 1987.  doi: 10.2307/2288805.
  * Schölkopf et al. (2012) B. Schölkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and J. Mooij.  On causal and anticausal learning.  _In Proceedings of the 29th International Conference on Machine Learning (ICML)_ , 2012.  URL http://www.jstor.org/stable/10.1086/525638.
  * Shi et al. (2014) P. Shi, P. Xie, M. Qin, F. Si, K. Dou, and K. Du.  Cluster analysis for daily patterns of SO2 and NO2 measured by the DOAS system in Xiamen.  _Aerosol and Air Quality Research_ , 2014.  doi: 10.4209/aaqr.2013.05.0160.
  * Shimizu et al. (2006) S. Shimizu, P. Hoyer, A. Hyvärinen, and A. Kerminen.  A linear non-gaussian acyclic model for causal discovery.  _Journal of Machine Learning Research_ , 7:2003–2030, 10 2006.
  * Smith (1990) R.L. Smith.  Extreme value theory.  _Handbook of Applicable Mathematics_ , (7):437–447, 1990.
  * Spirtes et al. (2001) P. Spirtes, C. Glymour, and R. Scheines.  _Causation, Prediction, and Search, 2nd Edition_ , volume 1.  The MIT Press, 1 edition, 2001.  URL https://EconPapers.repec.org/RePEc:mtp:titles:0262194406.
  * Tawn (1988) J.A. Tawn.  Bivariate extreme value theory: Models and estimation.  _Biometrika_ , 75(3):397–415, 1988.  doi: 10.2307/2336591.
  * Tawn (1990) J.A. Tawn.  Modelling multivariate extreme value distributions.  _Biometrika_ , 77(2):245–253, 1990.  doi: 10.2307/2336802.
  * Varian (2016) H. R. Varian.  Causal inference in economics and marketing.  _Proceedings of the National Academy of Sciences_ , 113(27):7310–7315, 2016.  doi: 10.1073/pnas.1510479113.
  * Verma and Pearl (2013) T. S. Verma and J. Pearl.  On the equivalence of causal models.  _CoRR_ , 2013.  URL http://arxiv.org/abs/1304.1108.
  * Wood et al. (2016) S. N. Wood, N. Pya, and B. Säfken.  Smoothing parameter and model selection for general smooth models.  _Journal of the American Statistical Association_ , 111(516):1548–1563, 2016.  doi: 10.1080/01621459.2016.1180986.
  * Wood (2017) S.N. Wood.  _Generalized Additive Models: An Introduction with R_.  Chapman and Hall/CRC, 2 edition, 2017.
  * Yoshida and Naito (2012) T. Yoshida and K. Naito.  Asymptotics for penalized splines in generalized additive models.  _Journal of Nonparametric Statistics_ , 26, 08 2012.  doi: 10.1080/10485252.2014.899360.
  * Youngman (2020) B. Youngman.  Generalised additive extreme value models.  2020\.  URL https://CRAN.R-project.org/package=evgam.