Causal discovery in heavy-tailed models

TL;DR

This paper introduces a novel causal discovery method tailored for heavy-tailed data, leveraging extremal dependence asymmetries to infer causal structures even with latent variables, and demonstrates its effectiveness through theoretical guarantees and empirical tests.

Contribution

It proposes the causal tail coefficient for extremal dependence analysis and develops a consistent, efficient algorithm for causal inference in heavy-tailed models, extending causal discovery to extreme value contexts.

Findings

The causal tail coefficient reveals causal directions in heavy-tailed data.

The proposed algorithm consistently recovers causal order in simulations.

The method performs well on synthetic and real heavy-tailed datasets.

Abstract

Causal questions are omnipresent in many scientific problems. While much progress has been made in the analysis of causal relationships between random variables, these methods are not well suited if the causal mechanisms only manifest themselves in extremes. This work aims to connect the two fields of causal inference and extreme value theory. We define the causal tail coefficient that captures asymmetries in the extremal dependence of two random variables. In the population case, the causal tail coefficient is shown to reveal the causal structure if the distribution follows a linear structural causal model. This holds even in the presence of latent common causes that have the same tail index as the observed variables. Based on a consistent estimator of the causal tail coefficient, we propose a computationally highly efficient algorithm that estimates the causal structure. We prove that our method consistently recovers the causal order and we compare it to other well-established and non-extremal approaches in causal discovery on synthetic and real data. The code is available as an open-access R package.

Tables

Table 1. Table 1 : Summary of the possible values of Γ 12 subscript Γ 12 \Gamma_{12} and Γ 21 subscript Γ 21 \Gamma_{21} and the implications for causality.

	${\Gamma }_{21}=1$	${\Gamma }_{21}\in (1/2,1)$	${\Gamma }_{21}=1/2$
${\Gamma }_{12}=1$		($a$) $X_1$ causes $X_2$
${\Gamma }_{12}\in (1/2,1)$	($b$) $X_2$ causes $X_1$	($d$) common cause
${\Gamma }_{12}=1/2$			($c$) no causal link

Table 2. Table 2 : Summary of the possible values of Ψ 12 subscript Ψ 12 \Psi_{12} and Ψ 21 subscript Ψ 21 \Psi_{21} and the implications for causality.

	${\Psi }_{21}=1$	${\Psi }_{21}\in (1/2,1)$	${\Psi }_{21}=1/2$
${\Psi }_{12}=1$		($a$) $X_1$ causes $X_2$
${\Psi }_{12}\in (1/2,1)$	($b$) $X_2$ causes $X_1$	($d$) common cause
${\Psi }_{12}=1/2$			($c$) no causal link

Table 3. Table 3: The table refers to Section  5.1.2 . It displays the average SID and the standard error (SE) over 50 simulations, for ICA-LiNGAM and Pairwise LiNGAM. For each dimension p 𝑝 p we simulate n = 10000 𝑛 10000 n=10000 observations from a non-linear SCM with Student’s- t 𝑡 t noise with α = 3.5 𝛼 3.5 \alpha=3.5 degrees of freedom (see Section  5.1 ). We run each method both on the full dataset (‘Full dataset’) and on the partial dataset, where we keep only the observations in the tails of the distribution (‘Keep tails’). As dimension p 𝑝 p increases, the number of extreme observations in all their coordinates decreases exponentially. Therefore, for a given dimension p 𝑝 p , we say that an observation x ∈ ℝ p 𝑥 superscript ℝ 𝑝 x\in\mathbb{R}^{p} lies in the tails of the distribution if at least ⌊ p ⌋ 𝑝 \lfloor\sqrt{p}\rfloor of its coordinates are below (or above) the 10% (or 90%) quantile.

		Full dataset		Keep tails
		SID	SE	SID	SE
	$p=10$	0.236	0.029	0.110	0.021
	$p=20$	0.246	0.022	0.150	0.017
	$p=30$	0.231	0.017	0.148	0.014
ICA-LiNGAM	$p=50$	0.268	0.014	0.216	0.012
	$p=10$	0.090	0.024	0.003	0.001
	$p=20$	0.064	0.013	0.003	0.001
	$p=30$	0.026	0.007	0.005	0.003
Pairwise LiNGAM	$p=50$	0.019	0.004	0.006	0.002

Equations

$${\Gamma }_{12}:=\underset{u\to 1^{-}}{\lim }\text{mbox}E\left[F_2(X_2)\mid F_1(X_1)>u\right]\in [0,1],$$

$$X_j:=\sum _{k\in \mathrm{pa}(j)}{\beta }_{jk}{X}_k+{\epsilon }_j,j\in V,$$

$$\text{mbox}P(X>x)\sim \ell (x)x^{-\alpha },x\to \infty ,$$

$${\Gamma }_{jk}=\underset{u\to 1^{-}}{\lim }\text{mbox}E\left[F_k(X_k)\mid F_j(X_j)>u\right],$$

$$X_j:=\sum _{k\in \mathrm{pa}(j,G)}{\beta }_{jk}{X}_k+{\epsilon }_j,j\in V,$$

$$\text{mbox}P({\epsilon }_j>x)\sim c_j\ell (x)x^{-\alpha },x\to \infty .$$

$$X_j=\sum _{h\in \mathrm{An}(j,G)}{\beta }_{h\to j}{\epsilon }_h.$$

$$\text{mbox}P(X_j>x)\sim \sum _{h\in \mathrm{An}(j,G)}{\beta }_{h\to j}^{\alpha }\ell (x)x^{-\alpha },x\to \infty .$$

$${\Gamma }_{jk}=\frac{1}{2}+\frac{1}{2}\frac{{\sum }_{h\in A_{jk}}{\beta }_{h\to j}^{\alpha }}{{\sum }_{h\in \mathrm{An}(j,G)}{\beta }_{h\to j}^{\alpha }},$$

$${\Gamma }_{jk}=\frac{1}{2}+\frac{1}{2}\frac{e_j^T{\mathbf{H}}_{\alpha }{e}_{A_{jk}}}{e_j^T{\mathbf{H}}_{\alpha }{e}_{\mathrm{An}(j,G)}},$$

$${\beta }_{1\to 4}={\beta }_{42}{\beta }_{21}+{\beta }_{43}{\beta }_{31},{\beta }_{2\to 4}={\beta }_{42},{\beta }_{3\to 4}={\beta }_{43},{\beta }_{4\to 4}=1.$$

$${\Gamma }_{41}=\frac{1}{2}+\frac{1}{2}\frac{{\beta }_{1\to 4}^{\alpha }}{{\beta }_{1\to 4}^{\alpha }+{\beta }_{2\to 4}^{\alpha }+{\beta }_{3\to 4}^{\alpha }+{\beta }_{4\to 4}^{\alpha }}<1.$$

$${\hat{F}}_j(x)=\frac{1}{n}\sum _{i=1}^{n}1\left\{X_{ij}\le x\right\},x\in \mathbb{R},$$

$$g^{\leftarrow }(y)=\inf \left\{x\in \mathbb{R}:g(x)\ge y\right\},y\in \mathbb{R}.$$

$${\hat{\Gamma }}_{12}={\hat{\Gamma }}_{12}^{(n)}=\frac{1}{k}\sum _{i=1}^n{\hat{F}}_2(X_{i2})1\{X_{i1}>X_{(n-k),1}\}.$$

$$\underset{x\to \infty }{\lim }\frac{x{f}_j(x)}{1-F_j(x)}=\frac{1}{\gamma },\text{ for some }\gamma >0.$$

$$k_n\to \infty \text{and}{k}_n/n\to 0,n\to \infty .$$

$${\hat{\Gamma }}_{12}\stackrel{P}{⟶}{\Gamma }_{12}\text{ and }{\hat{\Gamma }}_{21}\stackrel{P}{⟶}{\Gamma }_{21}.$$

$$\text{mbox}P\left(\hat{\pi }\notin {\Pi }_G\right)\le p^2\underset{i,j\in V:i\ne j}{\max }\text{mbox}P\left(\left|{\hat{\Gamma }}_{ij}-{\Gamma }_{ij}\right|>\frac{1-\eta }{2}\right),$$

$$k_n\to \infty \text{and}{k}_n/n\to 0,n\to \infty .$$

$$\text{mbox}P\left(\hat{\pi }\notin {\Pi }_G\right)\to 0,\text{ as }n\to \infty .$$

$$\text{mbox}P({\epsilon }_j>x)\sim c_j^{+}\ell (x)x^{-\alpha },\text{mbox}P({\epsilon }_j<-x)\sim c_j^{-}\ell (x)x^{-\alpha },$$

$$\begin{array}{rl}{\Psi }_{jk}=& \underset{u\to 1^{-}}{\lim }\text{mbox}E\left[\sigma (F_k(X_k))\mid \sigma (F_j(X_j))>u\right],j,k\in V,\end{array}$$

$$\begin{array}{rl}{\Psi }_{jk}=& \underset{u\to 1^{-}}{\lim }\frac{1}{2}\text{mbox}E\left[\sigma (F_k(X_k))\mid F_j(X_j)>u\right]\\ & +\underset{u\to 0^{+}}{\lim }\frac{1}{2}\text{mbox}E\left[\sigma (F_k(X_k))\mid F_j(X_j)<u\right]\\ =& {\Psi }_{jk}^{+}+{\Psi }_{jk}^{-},\end{array}$$

$${\Psi }_{jk}=\frac{1}{2}+\frac{1}{4}\frac{{\sum }_{h\in A_{jk}}{c}_{hj}^{+}|{\beta }_{h\to j}{|}^{\alpha }}{{\sum }_{h\in \mathrm{An}(j,G)}{c}_{hj}^{+}|{\beta }_{h\to j}{|}^{\alpha }}+\frac{1}{4}\frac{{\sum }_{h\in A_{jk}}{c}_{hj}^{-}|{\beta }_{h\to j}{|}^{\alpha }}{{\sum }_{h\in \mathrm{An}(j,G)}{c}_{hj}^{-}|{\beta }_{h\to j}{|}^{\alpha }},$$

$$c_{hj}^{+}=\{\begin{array}{l}{c}_h^{+},{\beta }_{h\to j}>0,\\ c_h^{-},{\beta }_{h\to j}<0,\end{array}{c}_{hj}^{-}=\{\begin{array}{l}{c}_h^{-},{\beta }_{h\to j}>0,\\ c_h^{+},{\beta }_{h\to j}<0.\end{array}$$

$$k_n\to \infty \text{and}{k}_n/n\to 0,n\to \infty .$$

$$\text{mbox}P\left(\hat{\pi }\notin {\Pi }_G\right)\to 0,\text{ as }n\to \infty .$$

$$\pi (i)<\pi (j)⟹j\notin \mathrm{an}(i,G),\text{for all }i,j\in V_O.$$

$$\text{mbox}P({\epsilon }_j>x)\sim {\ell }_j(x)x^{-{\alpha }_j},x\to \infty .$$

Full text

1 Introduction and background

Reasoning about the causal structure underlying a data generating process is a key scientific question in many disciplines. In recent years, much progress has been made in the formalisation of causal language [36, 50, 25]. In several situations, causal relationships manifest themselves only in extreme events. As stated by Cox and Wermuth [9, Sec. 8.7], large interventions (also named natural experiments) often carry information that is likely to be causal. In this light, one can view extreme observations as natural experiments that perturb a given system and facilitate causal analysis. On the one hand, existing causal methodology, focuses on moment related quantities of the distribution and is not tailored to estimate causal relationships from extreme events. On the other hand, the statistics of univariate extremes is relatively well understood [42] and there is a large set of tools for the analysis of heavy-tailed distributions. This work attempts to bring the fields of causality and extremes closer.

Let us first consider a bivariate random vector $(X_{1},X_{2})$ and assume that we are interested in the causal relationship between the two random variables, $X_{1}$ and $X_{2}$. We consider a linear structural causal model [36, Sec. 1.4] over variables including $(X_{1},X_{2})$ (without feedback mechanisms). We can then distinguish between the six scenarios of causal configurations shown in Figure 1 that include $X_{1}$, $X_{2}$ and possibly a third unobserved random variable $X_{0}$. This collection is complete in the following sense. Any structural causal model including $X_{1}$ and $X_{2}$ is interventionally equivalent [e.g., 38, Sec. 6.8] to one of the examples shown in Figure 1 when taking into account interventions on $X_{1}$ or $X_{2}$ only. The dashed edges can be interpreted as directed paths induced by a linear structural causal model (SCM); see Section 1.1 for a formal definition.

Here and in the sequel, we say that “$X_{1}$ is the cause of $X_{2}$” or “$X_{1}$ causes $X_{2}$” if there is a directed path from $X_{1}$ to $X_{2}$ in the SCM’s underlying directed acyclic graph (DAG). Assume that $X_{1}$ is the cause of $X_{2}$, and that both variables are heavy-tailed. Intuitively, if the causal relationship is monotonic, then an extremely large value of $X_{1}$ should cause an extreme value of $X_{2}$. The causal direction should, therefore, be strongly visible in the largest absolute values of the random variables. Here, “extreme” is to be seen in the respective scale of each variable, so it will make sense to consider the copula $\{F_{1}(X_{1}),F_{2}(X_{2})\}$, where $F_{j}$ is the marginal distribution of $X_{j}$, $j=1,2$. To exploit this intuition, we define the causal tail coefficient between variables $X_{1}$ and $X_{2}$ as

[TABLE]

if the limit exists. It reflects the causal relationship between $X_{1}$ and $X_{2}$ since, intuitively, if $X_{1}$ has a monotonically increasing causal influence on $X_{2}$, we expect $\Gamma_{12}$ to be close to one. Conversely, extremes of $X_{2}$ will not necessarily lead to extremes of $X_{1}$ and therefore, the coefficient $\Gamma_{21}$, where the roles of $X_{1}$ and $X_{2}$ in (1) are reversed, may be strictly smaller than one. This asymmetry will be made precise in Theorem 1 for linear structural causal models with heavy-tailed noise variables, and it forms the basis of our causal discovery algorithm.

A different perspective of our approach goes beyond the usual way of defining causality in the bulk of distribution. Namely, by looking at the signal in the tails, we might recover an extremal causal mechanism that is not necessarily present in the central part of the distribution. An example of this can be observed in financial markets. During calm periods, it is not clear whether any causal relationship exists among the financial variables. However, during turmoil, it is common to observe one specific stock or sector causing very negative (or positive) returns of other stocks (or sectors), displaying an extremal causal mechanism. In the finance literature, the concept of extremal causal mechanism is explained in terms of contagion, i.e., the spread of shocks across different markets [see 17]. For example, Rodriguez [44] uses a copula model with time-varying parameters to explain such contagion phenomena within countries in Latin America and Asia. On the other hand, there are also applications where the causal mechanism is present in the bulk of the distribution but absent in the tails. Seneviratne et al. [46] present this type of causal relationship between air temperature and the evapotranspiration of soil moisture. As the air temperature increases, the evapotranspiration process increases, too. This continues until the soil moisture resources are reduced to the point that a further increase in the temperature has no causal effect on the evapotranspiration.

Heavy-tailed distributions are an example of non-Gaussian models, which have received some attention in the causal literature. The LiNGAM algorithm [47] exploits non-Gaussianity through independent component analysis [8] to estimate an underlying causal structure. Misra and Kuruoglu [32] consider stable noise variables in a Bayesian network and develop a structure learning algorithm based on BIC. The work of Gissibl et al. [19] also studies causal questions related to extreme events. They consider max-linear models [18] where only the largest effect propagates to the descendants in a Bayesian network. The work by Naveau et al. [33] falls within the domain of attribution science. Namely, by studying extreme climate events, they try to answer counterfactual questions such as “what the Earth’s climate might have been without anthropogenic interventions”. Mhalla et al. [31] develop a method to estimate the causal relationships between bivariate extreme observations. It relies on the Kolmogorov complexity concept [see 28] adapted to high conditional quantiles. [13] reviews recent work on causality and sparsity in extremes.

The rest of the paper is organised as follows. Sections 1.1 and 1.2 briefly review structural causal models and some important concepts from extreme value theory. Section 2 contains a causal model for heavy-tailed distributions with positive coefficients. We prove that, given the underlying distribution, the causal tail coefficient allows us to distinguish between the causal scenarios shown in Figure 1. In Section 3, we introduce an algorithm named extremal ancestral search (EASE) that can be applied to a matrix of causal tail coefficients and that retrieves the causal order of the true graph, in the population case. We prove that our algorithm estimates a causal order even in the case where the causal tail coefficients are estimated empirically from data, as the sample size tends to infinity. In Section 4, we first generalise the results of the previous sections to the case of a structural causal model with real-valued coefficients. To do that, we introduce a more general causal tail coefficient that is sensitive to both the upper and lower tail of the variables. Second, we discuss the robustness properties of EASE in the presence of hidden confounders. Third, we analyse the properties of the causal tail coefficient when the noise variables have different tail indices. Section 5 contains experiments on simulated data and real-world applications.

The Appendix consists of six sections. Section A summarises important facts about regularly varying random variables. Section B contains the proofs of the results of the paper. Section C illustrates how the EASE algorithm retrieves a causal order of a DAG. Section D describes the settings used in the simulation study. Section E contains additional figures and tables. Section F presents further results for Section 5.2.

1.1 Structural causal models

A linear structural causal model, or linear SCM, (3; 36, Sec. 1.4) over variables $X_{1},\ldots,X_{p}$ is a collection of $p$ assignments

[TABLE]

where $\mathrm{pa}(j)\subseteq V=\{1,\ldots,p\}$ and $\beta_{jk}\in\mathbb{R}\setminus\{0\}$, together with a joint distribution over the noise variables $\varepsilon_{1},\ldots,\varepsilon_{p}$. Here, we assume that the noise variables are jointly independent and that the induced graph $G=(V,E)$, obtained by adding directed edges from the parents $\mathrm{pa}(j)$ to $j$, is a directed acyclic graph (DAG) with nodes $V$ and (directed) edges $E\subset V\times V$. To ease notation, we adopt the convention to sometimes identify a node with its corresponding random variable. To highlight the fact that $\mathrm{pa}(j)$ depends on a specific DAG $G$, we write $\mathrm{pa}(j,G)$, $j\in V$.

Structural causal models describe not only observational distributions but also interventional distributions. An intervention on $X_{j}$, for example, is defined as replacing the corresponding assignment (2) while leaving the other equations as they were. In practice, a causal model can be falsified via randomised experiments [see, e.g., 38, Sec. 6.8].

We define a directed path between node $j$ and $k$ as a sequence of distinct vertices such that successive pairs of vertices belong to the edge set $E$ of $G$. If there is a directed path from $j$ to $k$, we say that $j$ is an ancestor of $k$ in $G$. The set of ancestors of $j$ is denoted by $\mathrm{An}(j,G)$, and we write $\mathrm{an}(j,G)=\mathrm{An}(j,G)\setminus\{j\}$ when we consider the ancestors of $j$ except itself. A node $j$ that has no ancestors, i.e., $\mathrm{an}(j,G)=\varnothing$, is called a source node (or root node). Given two nodes $j,k\in V$, we say that $X_{j}$ causes $X_{k}$ if there is a directed path from $j$ to $k$ in $G$. Furthermore, given nodes $i,j,k\in V$, we say that $X_{i}$ is a confounder (or common cause) of $X_{j}$ and $X_{k}$ if there is a directed path from node $i$ to node $j$ and $k$ that does not include $k$ and $j$, respectively. Whenever a confounder is unobserved, we say it is a hidden confounder or hidden common cause. Finally, if $\mathrm{An}(j,G)\cap\mathrm{An}(k,G)=\varnothing$, then we say that there is no causal link between $X_{j}$ and $X_{k}$. A graph $G_{1}=(V_{1},E_{1})$ is called a subgraph of $G$ if $V_{1}\subseteq V$ and $E_{1}\subseteq(V_{1}\times V_{1})\cap E$. Recall that any subgraph of a DAG $G$ is also a DAG. For details on graphical models, we refer to Lauritzen [29].

1.2 Regularly varying functions and random variables

A positive, measurable function $f$ is said to be regularly varying with index $\alpha\in\mathbb{R}$, $f\in\mbox{RV}_{\alpha}$, if it is defined on some neighbourhood of infinity $[x_{0},\infty)$, $x_{0}>0$, and if for all $c>0$, $\lim_{x\to\infty}{f(cx)}/{f(x)}=c^{\alpha}$. If $\alpha=0$, $f$ is said to be slowly varying, $f\in\mbox{RV}_{0}$.

A random variable $X$ is said to be regularly varying with index $\alpha$ if

[TABLE]

for some $\ell\in\mbox{RV}_{0}$, where for any function $f$ and $g$, we write $f\sim g$ if $f(x)/g(x)\to 1$ as $x\to\infty$. If $X$ is regularly varying with index $\alpha$ then $cX$ is also regularly varying with the same index, for any $c>0$. For example, random variables with a Student’s-$t$, Pareto, Cauchy, or Fréchet distribution are regularly varying.

A characteristic property of regularly varying random variables is the max-sum-equivalence. The idea is that large sums of independent random variables tend to be driven by only one single large value. For this reason, the tail of the distribution of the maximum is equal to the tail of the distribution of the sum. For a rigorous formulation see Appendix A. We refer to Embrechts et al. [11, Sec. A3] for further details on regular variation and max-sum-equivalence.

2 The causal tail coefficient

To measure the causal effects in the extremes, we define the following parameter.

Definition 1**.**

Given two random variables $X_{1}$ and $X_{2}$, we define the causal tail coefficient

[TABLE]

if the limit exists, for $j,k=1,2$ and $j\neq k$.

The coefficient $\Gamma_{jk}$ lies between zero and one and is invariant under any marginal strictly increasing transformation since it depends on the rescaled margins $F_{j}(X_{j})$, for $j=1,2$. Below, we lay down the setup.

2.1 Setup

Consider a linear structural causal model (SCM) with an induced directed acyclic graph (DAG) $G$,

[TABLE]

where we assume that the coefficients $\beta_{jk}$ are strictly positive, $j,k\in V$. Let the independent noise variables $\varepsilon_{1},\ldots,\varepsilon_{p}$ be real-valued and regularly varying with comparable tails, i.e., there exists a tail-index $\alpha>0$ and $\ell\in\mbox{RV}_{0}$ such that for all $j\in V$, there exists $c_{j}>0$ that satisfies

[TABLE]

To simplify the notation, we rescale the variables $X_{j}$ such that $c_{j}=1$, $j\in V$. Furthermore, denote by $\beta_{k\to j}$ the sum of distinct weighted directed paths from node $k$ to node $j$, with $\beta_{j\to j}:=1$. Since $G$ is acyclic, we can express recursively each variable $X_{j}$, $j\in V$, as a weighted sum of the noise terms $\varepsilon_{1},\dots,\varepsilon_{k}$ that belong to the ancestors of $X_{j}$, that is,

[TABLE]

The noise terms in (5) are independent and regularly varying with comparable tails as in (4). Therefore, by using Lemma 3 of Appendix A, we can write

[TABLE]

Any probability distribution induced by an SCM is Markov with respect to the induced DAG $G$, and thus we can read off statistical independencies from it by $d$-separation (30; 36, Sec. 1.2.3). Conversely, to infer dependencies directly from the graph, one needs to assume that the distribution is faithful to the DAG $G$ [see 50, Sec. 2.3.3]. Most causal methods based on restricted SCMs, such as LiNGAM [see 47], RESIT [41], and Peters and Bühlmann [39], do not assume faithfulness. Similarly, in this work we require the milder assumption that $\beta_{j\to k}$ is non-zero if $j$ is an ancestor of $k$, i.e., $X_{j}$ causes $X_{k}$. This is automatically satisfied if the SCM has positive coefficients. In the sequel, we refer to this model as a heavy-tailed linear SCM.

2.2 Causal structure and the causal tail coefficient

In the setting of Section 2.1, the causal tail coefficient always exists and carries information about the underlying causal structure. In particular, it can be expressed in closed form.

Lemma 1**.**

Consider a heavy-tailed linear SCM over $p$ variables. Then, for $j,k\in V$ and $j\neq k$,

[TABLE]

where $A_{jk}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)$, and the sum over an empty index set equals zero.

For a proof see Appendix B.1. Lemma 1 provides a closed form expression for the causal tail coefficient $\Gamma_{jk}$, which can be written as a sum of two terms. The first term corresponds to the case when $X_{j}$ and $X_{k}$ are independent. The second term is non-negative and depends on the coefficients of the SCM and the tail index $\alpha>0$. By using matrix notation, it is possible to express the coefficient $\Gamma_{jk}$ more compactly. Consider the matrix of coefficients $\mathbf{B}$ of the DAG $G$, where $\mathbf{B}_{jk}=\beta_{jk}$, $j,k\in V$, and let $\mathbf{I}$ be the identity matrix. Furthermore, for any $\mathbf{M}\in\mathbb{R}^{p\times p}$, denote by $\mathbf{M}_{\alpha}$ the matrix where each entry of $\mathbf{M}$ is raised to the power $\alpha$. By applying the Neumann series, we obtain $\mathbf{H}=(\mathbf{I}-\mathbf{B})^{-1}$ where $\mathbf{H}_{jk}=\beta_{k\to j}$ for $j,k\in V$. Therefore, for $j,k\in V$ and $j\neq k$, we can write the causal tail coefficient as

[TABLE]

where $e_{j}\in\mathbb{R}^{p}$ is the $j$-th standard basis vector, and $e_{C}=\sum_{j\in C}e_{j}\in\mathbb{R}^{p}$ for any set $C\subseteq\{1,\dots,p\}$.

Example 1*.*

Consider the “diamond” graph $G=(V,E)$ in Figure 2, with $V=\{1,\dots,4\}$.

In this graph, for instance, it is easy to see that $\Gamma_{14}=1$. To compute $\Gamma_{41}$, we list the weighted directed paths from the ancestors of node 4 to node 4 itself, i.e.,

[TABLE]

Additionally, the set of common ancestors of node 1 and 4 is $A_{14}=\{1\}$. Putting everything together, by using Lemma 1, or formula (6), we obtain

[TABLE]

In this example, we see that it is possible to infer the causal relationship between $X_{1}$ and $X_{4}$ because $\Gamma_{14}>\Gamma_{41}$.

Consider now a general, heavy-tailed linear SCM over $p$ variables including $X_{1}$ and $X_{2}$, and inducing graph $G$. The following theorem shows that the causal tail coefficient, which is computable from the bivariate distribution of $X_{1}$ and $X_{2}$, see Equation (3), encodes the causal relationship between the two variables.

Theorem 1**.**

Consider a heavy-tailed linear SCM over $p$ variables including $X_{1}$ and $X_{2}$, as defined in Section 2.1. Then, knowledge of $\Gamma_{12}$ and $\Gamma_{21}$ allows us to distinguish the following cases: (a) $X_{1}$ causes $X_{2}$, (b) $X_{2}$ causes $X_{1}$, (c) there is no causal link between $X_{1}$ and $X_{2}$ (i.e., $\mathrm{An}(1,G)\cap\mathrm{An}(2,G)=\varnothing$), (d) there is a node $j\not\in\{1,2\}$, such that $X_{j}$ is a common cause of $X_{1}$ and $X_{2}$ and neither $X_{1}$ causes $X_{2}$ nor $X_{2}$ causes $X_{1}$. The corresponding values for $\Gamma_{12}$ and $\Gamma_{21}$ are depicted in Table 1.

For a proof see Appendix B.2. This result will also play a key role when estimating causal relationships from finitely many data. As a first remark, condition (a) and (b) might also include the presence of a common cause $X_{j}$. As a second remark, the empty entries in Table 1 cannot occur under the assumptions made in Section 2.1. For example, $\Gamma_{12}=\Gamma_{21}=1$ can only happen if some variables have different tail indices. One possibility is when the cause has a heavier tail than the effect. Another scenario is when a common cause $X_{j}$, for some $j\neq 1,2$, has heavier tails than the confounded variables $X_{1}$ and $X_{2}$. For further discussion on different tail indices see Section 4.3.

2.3 A non-parametric estimator

Consider a heavy-tailed linear SCM over $p$ variables including $X_{1}$ and $X_{2}$, with distributions $F_{1}$ and $F_{2}$, as described in Section 2.1. In order to construct a non-parametric estimator of $\Gamma_{12}$ and $\Gamma_{21}$ based on independent observations $(X_{i1},X_{i2})$, $i=1,\dots,n$, of $(X_{1},X_{2})$, we define the empirical distribution function of $X_{j}$ as

[TABLE]

for $j=1,2$. Denote by $g^{\leftarrow}$ the left continuous generalised inverse

[TABLE]

In addition, let the $(n-k)$-th order statistics be denoted by $X_{(n-k),j}=\widehat{F}_{j}^{\leftarrow}(1-k/n)$, for all $k=0,\dots,n-1$ and $j=1,2$, such that $X_{(1),j}\leq\dots\leq X_{(n),j}$. Replacing $F_{1}$ and $F_{2}$ in the definition of $\Gamma_{12}$ in (3) by the empirical counterparts, and the threshold $u$ by $u_{n}=1-k/n$, for some integer $0<k\leq n-1$, we define the estimator

[TABLE]

For this estimator to be consistent, a classical assumption in extreme value theory is that the number of upper order statistics $k=k_{n}$ depends on the sample size $n$ such that $k_{n}\to\infty$ and $k_{n}/n\to 0$ as $n\to\infty$. The first condition is needed to increase the effective sample size, whereas the second condition eliminates the approximation bias. The estimator $\widehat{\Gamma}_{21}=\widehat{\Gamma}_{21}^{(n)}$ is defined in an analogous way as (8).

Theorem 2**.**

Let $X_{i1}$ and $X_{i2}$, $i=1,\dots,n$, be independent copies of $X_{1}$ and $X_{2}$, respectively, where $X_{1}$ and $X_{2}$ are two of the $p$ variables of a heavy-tailed linear SCM.

• (A1)

Assume that the density functions $f_{j}=F^{\prime}_{j}$, $j=1,2$, exist and satisfy the von Mises’ condition

[TABLE]

• (A2)

Let $k_{n}\in\mathbb{N}$ be an intermediate sequence with

[TABLE]

Then the estimators $\widehat{\Gamma}_{12}$ and $\widehat{\Gamma}_{21}$ are consistent, as $n\to\infty$, i.e.,

[TABLE]

Remark 1**.**

The von Mises’ condition in (A1) is a very mild assumption that is satisfied by most univariate regularly varying distributions of interest. In our case $\gamma=1/\alpha$, where $\alpha$ is the common tail index of the noise variables.

For a proof of Theorem 2 see Appendix B.3. It uses several results from tail empirical process theory [e.g., 10, Sec. 2.2]. The main challenge comes from the fact that the variables $X_{1}$ and $X_{2}$ are tail dependent, and that the use of the empirical distribution function $\widehat{F}_{2}$ in (8) introduces dependence between the terms corresponding to different observations $i=1,\dots,n$. A related problem is studied in Cai et al. [5], where they derive asymptotic properties of the empirical estimator of the expected shortfall when another dependent variable is extreme. However, in contrast to Cai et al. [5], in the proof of Theorem 2, we work with a more explicit model and we consider the variables scaled to uniform margins, i.e., $\widehat{F}_{j}(X_{j})$ instead of $X_{j}$, $j=1,2$.

3 Causal discovery using extremes

We would like to recover the causal information from a dataset of $p$ variables under the model specification of Section 2.1. We develop an algorithm named extremal ancestral search (EASE) based on the causal tail coefficient defined in (3). We show that EASE can recover the causal order of the underlying graph in the population case (Section 3.1), and that it is consistent (Section 3.2).

3.1 Learning the causal order

Our goal is to recover the causal order of a heavy-tailed linear SCM over $p$ variables (as defined in Section 2.1) by observing $n$ i.i.d. copies of the random vector $X\in\mathbb{R}^{p}$. Given a DAG $G=(V,E)$, a permutation $\pi:\{1,\dots,p\}\rightarrow\{1,\dots,p\}$ is said to be a causal order (or topological order) of $G$ if $\pi(i)<\pi(j)$ for all $i$ and $j$ such that $i\in\mathrm{an}(j,G)$. We denote by $\Pi_{G}$ the set of all causal orders of $G$. For a permutation $\pi$ we sometimes use the notation $\pi=\left(\pi(1),\dots,\pi(p)\right)$.

A given causal order $\pi$ does not specify a unique DAG. As an example, the causal order $\pi=(1,2)$ comprises two DAGs: one where there is a directed edge between node 1 and node 2, and one where the two nodes are unconnected. On the other hand, there can be several causal orders for a given DAG. For example, a fully-disconnected DAG satisfies any causal order. However, even if the causal order does not identify a unique DAG, it still conveys important information. In particular, each causal order defines a class of DAGs that agree with respect to the non-ancestral relations. Therefore, once a causal order is available, one can estimate the complete DAG by using regularised regression methods. This idea has been exploited, e.g., by Shimizu et al. [48] and Bühlmann et al. [4]. In addition, Bühlmann et al. [4] and Peters and Bühlmann [37] argue that knowledge of a causal order is useful per se. In fact, given a correct causal order, one can construct a fully-connected DAG that describes interventional distribution across the variables.

For any heavy-tailed linear SCM and induced DAG $G=(V,E)$, we define the matrix $\Gamma\in\mathbb{R}^{p\times p}$ with entries $\Gamma_{ij}$, the causal tail coefficients between all pairs of variables $X_{i}$ and $X_{j}$, $i,j\in V$; see Definition 1. Theorem 1 tells us how the entries of $\Gamma$ encode the causal relationships between the random variables of the SCM. To recover the causal order of the DAG $G$, we propose the algorithm named extremal ancestral search (EASE).

EASE algorithm

Input: A matrix $\Gamma\in\mathbb{R}^{p\times p}$ of causal tail coefficients related to a DAG $G=(V,E)$ with $V=\{1,\dots,p\}$.

Returns: Permutation of the nodes $\pi:V\to\{1,\dots,p\}$.

Set $V_{1}=V$.

For $s\in\{1,\dots,p\}$

(a)

Let $M^{(s)}_{i}=\max_{j\in V_{s}\setminus\{i\}}\Gamma_{ji}$, for all $i\in V_{s}$.

(b)

Let $i_{s}\in\arg\min_{i\in V_{s}}M^{(s)}_{i}$.

(c)

Set $\pi(i_{s})=s$.

(d)

Set $V_{s+1}=V_{s}\setminus\{i_{s}\}$.

Return the permutation $\pi$.

Complexity: $O(p^{2})$.

The above is a greedy algorithm that identifies root nodes of the current subgraph at each step. In the first step, the algorithm finds a root node $i_{1}\in V$ as the one that minimises the score $M_{i}^{(1)}=\max_{j\neq i}\Gamma_{ji}$, $i\in V$. In fact, by Theorem 1, $M_{i}^{(1)}<1$ if and only if $i$ is a source node. Once the first node is selected, the algorithm searches for a second root node in the subgraph where $i_{1}$ is removed. The procedure continues until all nodes have been selected. In Appendix C, one example illustrates how EASE finds a causal order for a given DAG.

The next result states that, in the population case, the EASE algorithm yields a correct causal order of the underlying DAG.

Proposition 1**.**

Consider a heavy-tailed linear SCM over $p$ variables, as defined in Section 2.1, and let $G=(V,E)$ be the induced DAG. If the input $\Gamma$ is the matrix of causal tail coefficients associated with the SCM, then EASE returns a permutation $\pi$ that is a causal order of $G$.

For a proof see Appendix B.4.

3.2 Sample properties for the EASE algorithm

For finite samples, the EASE algorithm will take an estimate of the causal coefficient matrix $\Gamma$ as input. Based on the empirical non-parametric estimator $\widehat{\Gamma}$ and its asymptotic properties, we assess the performance of the algorithm. Let $\widehat{\Gamma}\in\mathbb{R}^{p\times p}$ denote the matrix where each entry $\widehat{\Gamma}_{ij}$ is defined as in (8) in Section 2.3, for $i,j\in V$. We say that a procedure makes a mistake when it returns a permutation $\pi\notin\Pi_{G}$. We derive an upper bound for the probability that EASE makes a mistake when the matrix $\Gamma$ is estimated by $\widehat{\Gamma}$.

Proposition 2**.**

Consider a heavy-tailed linear SCM over $p$ variables $X=(X_{1},\ldots,X_{p})$, as defined in Section 2.1, with induced DAG $G$. Let $\widehat{\Gamma}$ be the estimated causal coefficient matrix related to $G$. Let $\widehat{\pi}$ denote the permutation returned by EASE based on $\widehat{\Gamma}$. Then,

[TABLE]

where $\eta=\max_{u\notin\mathrm{An}(v,G)}\Gamma_{uv}<1$.

For a proof see Appendix B.5. The bound for the probability of making a mistake in the estimated causal order is expressed in terms of the distance between the true $\Gamma_{ij}$ and the estimated $\widehat{\Gamma}_{ij}$. This bound in combination with the consistency result of Theorem 2 yields the consistency of the EASE algorithm in the sample case.

Corollary 1**.**

Let $\widehat{\pi}$ be the permutation computed by EASE under the assumptions of Proposition 2. Let $k_{n}\in\mathbb{N}$ be an intermediate sequence with

[TABLE]

If the von Mises’ condition (9) holds, then the EASE algorithm is consistent, i.e.,

[TABLE]

The result above is for fixed dimension $p$. To prove consistency in a regime where $p$ scales with the sample size $n$ we would need to establish concentration inequalities for $\widehat{\Gamma}$ or asymptotic normality in Theorem 2. Both would require stronger assumptions on the tails of the noise variables and a second-order analysis in line with the proof of Theorem 2.

3.3 Computational complexity

The EASE algorithm is based on pairwise quantities and is therefore computationally efficient. To estimate the matrix $\widehat{\Gamma}$ of causal tail coefficients, which is the input for EASE, first, we need to rank the $n$ observations for each of the variables, with a computational complexity of $O(pn\log n)$. Then we compute the coefficients $\widehat{\Gamma}_{ij}$ for each pair $i,j\in V$, with a computational complexity of $O(k_{n}p^{2})$. The computational complexity of EASE grows with the square of the number of variables, i.e., $O(p^{2})$. The overall computational complexity of estimating the matrix $\widehat{\Gamma}$ and running the EASE algorithm is therefore $O\left(\max(pn\log n,k_{n}p^{2})\right)$.

4 Extensions

4.1 Real-valued coefficients

Until now, we have worked with a heavy-tailed linear SCM with positive coefficients (see Section 2.1 for a detailed explanation of the model). In the current section, we relax this assumption and let the coefficients of the SCM be real-valued, i.e., $\beta_{jk}\in\mathbb{R}$, $j,k\in V$. Additionally, we assume that $\beta_{j\to k}$ is non-zero, if $j$ is an ancestor of $k$. Given that the coefficients are real-valued, we want to consider both the upper and the lower tails of the variables. We assume that the noise variables $\varepsilon_{1},\dots,\varepsilon_{p}$, of the SCM have comparable upper and lower tails, that is, as $x\to\infty$

[TABLE]

where $c_{j}^{+},c_{j}^{-}>0$, $j\in V$ and $\ell\in\mbox{RV}_{0}$. Furthermore, we define a causal tail coefficient that is sensitive to both tails as

[TABLE]

if the limit exists, where $\sigma:x\mapsto|2x-1|$. Since $F_{j}(X_{j})\sim\text{Unif}[0,1]$, $j\in V$, we can rewrite (10) as

[TABLE]

where the first and second terms correspond to the cases where $X_{j}$ is extremely large and extremely small, respectively.

In the current setting, the coefficient defined in (10) always exists, and it has a closed form expression that encodes causal relationships between the variables.

Lemma 2**.**

Consider a heavy-tailed SCM over $p$ variables, where the coefficient $\beta_{jk}\in\mathbb{R}$, $j,k\in V$. Assume that $\beta_{j\to k}\neq 0$ if $j$ is an ancestor of $k$. Then, for $j,k\in V$ and $j\neq k$,

[TABLE]

where $A_{jk}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)$, and

[TABLE]

A proof is provided in Appendix B.6. The interpretation of the result is as follows. The baseline of the coefficient is $1/2$, which can be checked to be the value of $\Psi_{jk}$ when two variables are independent. The other two terms account for the equally weighted contribution from the lower and upper tail, respectively. The result stated in Lemma 2 allows us to extend Theorem 1 to the more general setting where the heavy-tailed SCM has coefficients $\beta_{jk}\in\mathbb{R}$, $j,k\in V$.

Theorem 3**.**

Consider a heavy-tailed linear SCM over $p$ variables including $X_{1}$ and $X_{2}$, and assume that $\beta_{jk}\in\mathbb{R}$, $j,k\in V$. In addition, assume that $\beta_{j\to k}\neq 0$ if $j$ is an ancestor of $k$. Then, knowledge of $\Psi_{12}$ and $\Psi_{21}$ allows us to distinguish the following cases: (a) $X_{1}$ causes $X_{2}$, (b) $X_{2}$ causes $X_{1}$, (c) there is no causal link between $X_{1}$ and $X_{2}$, (d) there is a node $j\not\in\{1,2\}$, such that $X_{j}$ is a common cause of $X_{1}$ and $X_{2}$ and neither $X_{1}$ causes $X_{2}$ nor $X_{2}$ causes $X_{1}$. The corresponding values for $\Psi_{12}$ and $\Psi_{21}$ are shown in Table 2.

The proof is identical to the proof of Theorem 1, replacing $\Gamma_{ij}$ with $\Psi_{ij}$ and by referring to Lemma 2 instead of Lemma 1. Moreover, as in Theorem 1, condition (a) and (b) can also include the presence of a common cause $X_{j}$. Theorem 3 implies that if we run the EASE algorithm based on the matrix $\Psi\in\mathbb{R}^{p\times p}$, containing the pairwise $\Psi_{ij}$, $i,j\in V$, then we retrieve a causal order of the underlying DAG. This is the analogue to Proposition 1 for heavy-tailed linear SCM with real-valued coefficients.

We define an empirical estimator $\widehat{\Psi}_{ij}$ of $\Psi_{ij}$ in a similar fashion as the estimator $\widehat{\Gamma}_{ij}$ in (8). The proof of Lemma 2 shows that the coefficient $\Psi_{ij}$ can be decomposed in the same way as $\Gamma_{ij}$ in the proof of Lemma 1. Therefore, following the lines of the proof of Theorem 2 with some minor modifications, we obtain the consistency $\widehat{\Psi}_{ij}\stackrel{{\scriptstyle P}}{{\longrightarrow}}\Psi_{ij}$ as $n\to\infty$ for any intermediate sequence $k_{n}\to\infty$ and $k_{n}/n\to 0$, and under the assumption of the von Mises’ condition for both the upper and the lower tail of $X_{j}$, $j\in V$.

We can then estimate a permutation $\widehat{\pi}$ by the EASE algorithm based on the matrix $\widehat{\Psi}\in\mathbb{R}^{p\times p}$ that contains the estimators $\widehat{\Psi}_{ij}$, $i,j\in V$, as entries. For this permutation we obtain the same bound for the probability that EASE makes a mistake as shown in Proposition 2 by replacing $\widehat{\Gamma}_{ij}$ and $\Gamma_{ij}$ by $\widehat{\Psi}_{ij}$ and $\Psi_{ij}$, respectively. This together with the consistency of $\widehat{\Psi}_{ij}$ yields the following result.

Corollary 2**.**

Assume the general setup of the heavy-tailed linear SCM with real-valued coefficients of this section. Let $\widehat{\pi}$ be the permutation computed by EASE based on the matrix $\widehat{\Psi}$. Assume the von Mises’ condition for the upper and the lower tail of $X_{j}$ and let $k_{n}\in\mathbb{N}$ be an intermediate sequence with

[TABLE]

Then, the EASE algorithm is consistent, i.e.,

[TABLE]

4.2 Presence of hidden confounders

A frequent assumption in causality is that one can observe all the relevant variables. However, in many real-world situations, it is hard, if not impossible, to do so. When some of the hidden variables are confounders (i.e., common causes), the causal inference process might be compromised. Therefore, an attractive property of a causal inference algorithm involves its robustness to hidden confounders. In this section, we show that EASE is capable of dealing with hidden common causes and, under certain assumptions, it recovers the causal order of the observed graph both in the population and in the asymptotic case.

Consider a heavy-tailed linear SCM with real-valued coefficients, as defined in Section 4.1, consisting of both observed and hidden variables. This SCM induces a DAG $G=(V,E)$, with $V=V_{O}\cup V_{H}$, $V_{O}\cap V_{H}=\varnothing$, where $V_{O}$ ($V_{H}$) denotes the set of nodes corresponding to the observed (hidden) variables. Our goal is to recover a causal order for the subset of the observed variables $X_{j}$, $j\in V_{O}$. In particular, we say that the EASE algorithm recovers a causal order $\pi$ over the observed variables if

[TABLE]

In fact, the results of the previous sections hold even in the presence of hidden confounders.

Regarding the population properties, Theorem 1 and 3 still apply: they state that the causal tail coefficients $\Gamma$ and $\Psi$ reflect the causal relationships between pairs of variables without taking into account other variables, e.g., by conditioning. In addition, the result of Proposition 1, and the corresponding extension in Section 4.1, are also valid. The proof of Proposition 1 depends only on the assumption that the input matrix contains the pairwise causal effects between the variables. Therefore, if we use matrix $\Gamma$ (or $\Psi$) as input for the EASE algorithm, we recover a causal order $\pi$ that satisfies (13).

Regarding the asymptotic properties, $\hat{\Gamma}$ and $\hat{\Psi}$ are consistent even in the presence of hidden common causes: in the proof of Theorem 2, the other variables do not appear. In addition, we can still find an upper bound for the probability that the EASE algorithm makes a mistake. To do so, one needs to adjust the proof of Proposition 2 by replacing the full DAG $G$ with the subgraph $G_{O}=(V_{O},E_{O})$ containing only the observed variables, where $E_{O}=E\cap(V_{O}\times V_{O})$. Combining the two previous arguments, it follows that Corollary 1 and 2 hold even in the presence of hidden confounders.

The ability to deal with hidden confounders is a property that, in general, is not shared by all methods in causality. For example, the PC algorithm [50, Sec. 5.4.2] might retrieve a Markov equivalence class that contains DAGs with a wrong causal order if some of the variables are not included in the analysis. Similarly, the standard version of the LiNGAM algorithm [47] might produce a wrong DAG in the presence of hidden common causes. Hoyer et al. [23], Entner and Hoyer [14], and Tashiro et al. [51] proposed extensions of LiNGAM that deal with hidden variables. While Entner and Hoyer [14], and Tashiro et al. [51] show good performance in practice, all three methods suffer from some drawbacks. For example, the LiNGAM version of Hoyer et al. [23] requires a priori the number of hidden variables in the SCM (or needs to estimate it from data). The main limitation of Entner and Hoyer [14] is that it recovers causal information only for subsets of variables that are not affected by hidden confounders. For some non-ancestral graphs, the method by Tashiro et al. [51] does not identify all ancestral relationships [52]. In addition, both the work of Hoyer et al. [23] and Tashiro et al. [51] are computationally intensive, with the latter showing a computational time that grows exponentially with the sample size and the number of observed variables. Among the constraint-based methods, Spirtes et al. [50, Sec. 6.7] proposed the FCI method, which is an extension to the PC algorithm that deals with arbitrarily many hidden confounders and produces a partial ancestral graph [see 53]. Due to the high number of independence tests, the FCI algorithm can be slow when the number of variables is large. For this reason, Claassen et al. [6] proposed the FCI+ algorithm, a faster version of FCI that is consistent in sparse high-dimensional settings with arbitrarily many hidden variables. In general, FCI type algorithms produce an equivalence class of graphs. They are not guaranteed to recover the causal order of the variables.

Compared to the methods mentioned above, our algorithm has the advantage of being computationally fast, and being able to produce a causal order without assumptions on the number of hidden variables and the sparsity of the true underlying DAG.

4.3 Noise variables with different tails

We have so far considered the case where the noise variables of a given SCM share the same tail coefficient $\alpha>0$ and the same slowly varying function $\ell$. Consider now a heavy-tailed SCM over $p$ variables, as defined in Section 2.1, with the difference that the noise variables have possibly different tail indices $\alpha_{1},\dots,\alpha_{p}>0$ and slowly varying functions $\ell_{1},\dots,\ell_{p}\in\mbox{RV}_{0}$, i.e., for $j=1,\dots,p$,

[TABLE]

We say that $\varepsilon_{j}$ has heavier (upper) tail than $\varepsilon_{k}$ if either $0<\alpha_{j}<\alpha_{k}$, or $\alpha_{j}=\alpha_{k}$ and $\ell_{j}(x)/\ell_{k}(x)\to\infty$ as $x\to\infty$. Denote by $G=(V,E)$ the DAG induced by the SCM. With similar arguments to the proof of Lemma 1, the causal tail coefficient for $j,k\in V$ can then be expressed as

[TABLE]

where $A_{jk}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)$, and the sum over an empty index set equals zero. From (14), we can study the different constellations of $X_{j}$ and $X_{k}$ and the corresponding values of $\Gamma_{jk}$. We obtain the following three statements.

If $X_{j}$ and $X_{k}$ are independent, the causal tail coefficient satisfies, as before, $\Gamma_{jk}=1/2$. 2. 2.

If $X_{j}$ is an ancestor of $X_{k}$ then, as before, $\Gamma_{jk}=1$. 3. 3.

In all other scenarios, it holds that $\Gamma_{jk}<1$ as long as the noise variables $\varepsilon_{h}$, $h\in A_{jk}$, of the common ancestors of $X_{j}$ and $X_{k}$ have tails that are lighter than (or as light as) the one of $\varepsilon_{j}$. On the other hand, if there is some common ancestor of $X_{j}$ and $X_{k}$ for which the noise variable’s tail is heavier than the one of $\varepsilon_{j}$, then $\Gamma_{jk}=1$.

These statements help to understand in which cases the values of $\Gamma_{jk}$ indicate a correct causal relation.

Example 2*.*

Suppose that $X_{j}$ is an ancestor of $X_{k}$, and there is possibly a common ancestor $X_{0}$ (which can also be a hidden confounder). Since $\Gamma_{jk}=1$, we will never mistakenly detect the existence of a causal effect from $X_{k}$ to $X_{j}$. If either $X_{j}$ or $X_{0}$ has a heavier tail than $X_{k}$, then $\Gamma_{kj}=1$ and we cannot detect the causal effect from $X_{j}$ to $X_{k}$.

Example 3*.*

Suppose neither $X_{j}$ causes $X_{k}$ nor $X_{k}$ causes $X_{j}$, and $X_{0}$ is a common ancestor of $X_{j}$ and $X_{k}$. If $X_{0}$ has a tail that is lighter than (or as light as) the one of $X_{k}$, but heavier than the one of $X_{j}$, then $\Gamma_{jk}=1>\Gamma_{kj}$. Therefore, the causal tail coefficient indicates a wrong causal effect from $X_{j}$ to $X_{k}$.

Whenever there exists a causal effect between two variables, we can, at worst, fail to detect it (that is, the causal tail coefficient does not indicate a causal effect in the wrong direction). When there is no causal connection between two variables, the causal tail coefficient might indicate a wrong causal effect. However, this does not affect the correctness of the EASE algorithm, on the population level. Indeed, if $\Gamma_{jk}=1>\Gamma_{kj}$, for $j,k\in V$, there are two possibilities. If $X_{j}$ is an ancestor of $X_{k}$, then the algorithm correctly chooses $j$ before $k$. If $X_{j}$ and $X_{k}$ share a common ancestor, but none of them is causing the other (see Example 3), then any permutation of $j$ and $k$ yields a valid causal order. Example 2 shows that EASE could make mistakes when $X_{j}$ is an ancestor of $X_{k}$ and $\Gamma_{jk}=\Gamma_{kj}=1$, since the causal tail coefficient does not indicate any causal effect. In this case, one could remove one variable at a time to obtain a subset $\tilde{V}\subset V$ that satisfies $\Gamma_{hm}<1$ or $\Gamma_{mh}<1$ for each $h,m\in\tilde{V}$. By applying the EASE algorithm to the subset of the remaining variables, one would recover a correct causal order on such subset.

For simplicity we only considered the causal tail coefficient $\Gamma$, but similar conclusions hold for $\Psi$.

5 Numerical results

5.1 Simulation study

We assess the performance of EASE in estimating a causal order of a graph induced by a heavy-tailed SCM. We simulate the SCMs with real-valued coefficients and different numbers of variables $p$ and samples $n$. The noise variables have Student’s $t$ distributions with different degrees of freedom $\alpha$ and we consider four different settings, including unobserved confounders and model misspecification; see Appendix D for details. Since the coefficients in the SCM are real-valued, we use the causal tail coefficient $\Psi$ defined in Section 4.1 for our EASE algorithm. Our code is available as an R package at https://github.com/nicolagnecco/causalXtreme. Scripts generating all our figures and results can be found at the same url.

5.1.1 Competing methods and evaluation metric

We compare our algorithm to three well-established methods in causality, the Rank PC algorithm [21], ICA-LiNGAM [47], and Pairwise LiNGAM [24].

The classic PC algorithm [50, Sec. 5.4.2] belongs to the class of constraint-based methods for causal discovery. It estimates the Markov equivalence class of a DAG, encoded as a completed partially directed acyclic graph (CPDAG). The PC algorithm retrieves a CPDAG by performing conditional independence tests on the variables. The Rank PC algorithm, proposed by Harris and Drton [21], is an extension of the PC algorithm and uses the rank-based Spearman correlation to perform the independence tests. This modification ensures that the method is more robust to non-Gaussian data.

The algorithms that fit our problem best are ICA-LiNGAM and Pairwise LiNGAM. ICA-LiNGAM, proposed by Shimizu et al. [47], leverages the results of independent component analysis (ICA) [8] to estimate the DAG of a linear SCM under the only assumption that the noise is non-Gaussian. Pairwise LiNGAM, proposed by Hyvärinen and Smith [24], is a likelihood-ratio-based method to identify the exogenous variables within the DirectLiNGAM framework. DirectLiNGAM, introduced by Shimizu et al. [48], is an algorithm based on two iterative steps, namely, finding an exogenous variable (i.e., a node in the DAG with no parents), and regressing this variable out of all the others. In this simulation study, we let ICA-LiNGAM and Pairwise LiNGAM return only a causal order (and not a complete DAG structure), in order to make a fair comparison with EASE.

The algorithms return different types of causal information. On the one hand, EASE, ICA-LiNGAM, and Pairwise LiNGAM estimate a causal order. On the other hand, the Rank PC algorithm computes a CPDAG that represents a Markov equivalence class of DAGs. Therefore, when it comes to evaluating the performance of the algorithms, it becomes crucial to use a measure that is meaningful for all of them. We choose the structural intervention distance (SID) proposed by Peters and Bühlmann [37]. The SID takes as input either a pair of DAGs or a DAG and a CPDAG and returns the number of falsely inferred interventional distributions [37, Definition 3]. We standardise the SID to lie between zero and one. For each method, we compute the distance between the simulated DAG, i.e., the ground truth, and the estimated DAG or CPDAG. An estimated causal order $\widehat{\pi}$ corresponds to a fully connected DAG $G=(V,E)$, where $(i,j)\in E$ if $\widehat{\pi}(i)<\widehat{\pi}(j)$. As a caveat, we slightly adapt the SID in the case of hidden confounders (see Setting 2 of our simulations), since it is not designed to work in such a situation.

5.1.2 Results

In this simulation experiment we use the implementation of the Rank PC, and ICA-LiNGAM algorithm developed by Kalisch et al. [27]. We implemented Pairwise LiNGAM in C++ and included it in our software package.

Regarding the hyperparameter settings, for the Rank PC algorithm, we perform a conditional independence test based on Spearman’s correlation coefficient, as proposed by Harris and Drton [21], and we set the level of the independence tests to 0.0005.

Concerning the choice of the number of exceedances $k_{n}$ in the EASE algorithm, we perform a small preliminary simulation. Figure 8 shows the SID of EASE for $k_{n}=\lfloor n^{\nu}\rfloor$ and different fractional exponents $\nu>0$. The best fractional exponent in Figure 8 seems to depend on the tail heaviness of the noise variables, and in particular it appears to be smaller for larger values of the degree of freedom $\alpha$ of the Student’s $t$ distribution. Our estimators $\widehat{\Gamma}_{ij}$ and $\widehat{\Psi}_{ij}$ are similar in construction to Hill’s estimator [22]. For the latter, the optimal number $k_{n}^{*}$ of exceedances depends on the tail index and an index related to a second-order condition; see Section 3.2 in [10] for details. For the Student’s $t$ distribution with $\alpha$ degrees of freedom it can be shown that $k_{n}^{*}\sim C_{\alpha}n^{1/(\alpha+1)}$, where $C_{\alpha}>0$ is a constant. This intuitive explanation coincides well with the optimal fractional exponents in Figure 8. In the sequel, we choose $k_{n}=\lfloor n^{0.4}\rfloor$ because it lies within the best range for the fractional exponent. This result also agrees with the assumptions of Theorem 2, where $k_{n}\to\infty$ and $k_{n}/n\to 0$, as $n\to\infty$.

Regarding the simulation settings, we let $n$ denote the number of observations, $p$ the number of variables, and $\alpha>0$ the tail index of the simulated distribution. For each combination of $n\in\{500,1000,10000\}$, $p\in\{4,7,10,15,20,30,50\}$ and $\alpha\in\{1.5,2.5,3.5\}$ we simulate 50 random SCMs under four different settings. The simulated data is independent of the data used to choose the best fractional exponent of $k_{n}$ (see Figure 8). The first setting corresponds to linear SCMs with real-valued coefficients described in Section 4.1. In the second setting, we introduce hidden confounders. The third setting corresponds to nonlinear SCMs. In the fourth setting, we first generate linear SCMs and then transform each variable to uniform margins. Further details on the generation of the SCMs are in Appendix D. For each simulation and setting we evaluate the performance of EASE, ICA-LiNGAM, Pairwise LiNGAM, and Rank PC algorithm with the SID. As a baseline, in each simulation, we also compute the SID of a randomly generated DAG, where we randomly choose the causal order, the sparsity and the edges of the graph.

Figure 3 displays the results of the simulations when the tail index $\alpha=1.5$. We can observe that EASE is quite robust across the four different settings. We explain this finding as follows. In the presence of hidden confounders (Setting 2), EASE can retrieve a correct causal order, asymptotically. Furthermore, the nonlinear setting used in this simulation (Setting 3) is such that the relationships between the variables are kept linear in the tails. Therefore, our algorithm is only moderately affected by this model misspecification. Finally, EASE is not affected by the transformation to uniform margins (Setting 4) because the causal tail coefficient $\Psi$ is invariant under any strictly monotone increasing transformation.

Compared to the other methods, we observe that EASE performs better than Rank PC across all settings, and better than ICA-LiNGAM in Setting 2 and 4. Pairwise LiNGAM is overall the best performing method, except in Setting 4. Also, Pairwise LiNGAM is less affected by misspecifications in the bulk of the data distribution (Setting 3), compared to ICA-LiNGAM. One reason is that Pairwise LiNGAM relies on ordinary least square regression that is sensitive to high-leverage points. In this particular setting, ICA-LiNGAM and Pairwise LiNGAM gain in robustness if we discard the data in the bulk of the distribution; see Table 3 in Appendix E. Furthermore, both ICA-LiNGAM and Pairwise LiNGAM are the algorithms with the best convergence for high dimension $p$, as $n$ increases. This result is not surprising because EASE uses only the $k_{n}<n$ upper order statistics to recover the causal structure. In addition, we notice that Pairwise LiNGAM outperforms ICA-LiNGAM, in agreement to the findings of Hyvärinen and Smith [24]. Regarding the Rank PC algorithm, we can see that it is quite stable under different settings, but it performs only marginally better than the random method. Moreover, in Figure 15 in Appendix E, we observe that the performance of the Rank PC algorithm is almost constant for significance levels of the independence test between $5\cdot 10^{-4}$ and $0.5$. The results do not change when we consider the standard PC algorithm, which is based on partial correlation as a conditional independence test. The results for tail indices $\alpha=2.5,3.5$ are in Figures 11 and 12 in Appendix E. Increasing values of $\alpha$ correspond to lighter tails, and we observe that it becomes more challenging for EASE to recover the correct causal order.

In addition to the competitive performance in the simulations, a further advantage of EASE is its computational efficiency. The algorithm performs computations only on the tails of the dataset and relies on simple non-parametric estimators of the causal tail coefficient; see Section 3.3. Figure 14 in Appendix E shows that EASE can be up to two orders of magnitude faster than the other methods.

5.2 Financial application

In general, one cannot easily reason about causality in financial markets. Several factors influence financial returns, and most of them are unobserved. In addition, the effect of these factors varies in time. However, under particular circumstances, it is possible to conjecture the existence of a specific causal relationship, with a reasonable degree of confidence. For example, in the Swiss financial market, one can argue that very large (both positive and negative) changes to the Euro Swiss franc exchange rate (EURCHF) induce changes in the Swiss Market Index (SMI), the main stock index in Switzerland. This is due to multiple reasons such as the multinational nature and the high export dependency of the Swiss economy. Consider, for instance, the decision of the Swiss National Bank (SNB) to discontinue the minimum exchange rate between Swiss franc and Euro, on January 15, 2015. This event can be deemed as a large intervention with a plausible causal interpretation (in the spirit of Cox and Wermuth [9, Sec. 8.7]). Following the SNB decision, the EURCHF plummeted more than 30 standard deviations, and all the stocks included in the SMI dropped in value on the same day.

For this reason, we consider the returns of the Euro Swiss franc exchange rate (EURCHF) and the three largest Swiss stocks in terms of market capitalisation, namely, Nestlé (NESN), Novartis (NOVN) and Roche (ROG). We choose to analyse three individual stocks instead of the SMI for three reasons. First, we deem it more appropriate to test our assumptions on more than two variables. Moreover, the three stocks make up 50% of the SMI composition, and thus they are representative of the index itself. Furthermore, all three companies are multinational corporations with a homogeneous exposure to foreign markets and a negligible fraction of revenues coming from the Swiss market [see 34, 35, 43]. The last point suggests that the effect (if any) of EURCHF on these stocks does not depend on the idiosyncrasies of each firm.

The dataset consists of daily returns spanning from January 2005 to September 2019 and includes $n=3832$ observations. The goal is to assess whether EASE can retrieve a correct causal order for the set of four variables. As ground truth, we conjecture that large changes in EURCHF (both positive and negative) will cause large changes in the stock returns, but not vice versa. Figure 4 shows the causal structure corresponding to our hypothesis.

Before running the EASE algorithm, we assess the tail behaviour of each variable by estimating the shape parameter $\xi$ of a generalised Pareto distribution [see 11, Sec. 3.4] on the threshold data. Recall that $\xi$ is the reciprocal of the tail index $\alpha=1/\xi$, if $\xi>0$. For each variable and each tail (upper and lower) we estimate the $\xi$ parameter using 200 observations, corresponding to the 95%-quantile, approximately. The estimated parameters and their standard errors for the upper tails are 0.31 (0.08) for EURCHF, 0.25 (0.08) for NESN, 0.16 (0.07) for NOVN, and 0.25 (0.09) for ROG. Regarding the lower tail, the estimated parameters and their standard errors are 0.27 (0.08) for EURCHF, 0.12 (0.08) for NESN, 0.17 (0.08) for NOVN, and 0.23 (0.10) for ROG. By adding and subtracting two standard errors from each point estimate of $\xi$, we observe that the lower shape parameter of Nestlé and Novartis is not significantly different from zero. Moreover, the shape parameters for the Euro Swiss franc exchange rate and for Roche are significantly different from zero in both tails. In addition, since the confidence intervals of all estimates are overlapping, the assumption of a common shape parameter seems reasonable. It seems however that the returns of Nestlé and Novartis have slightly lighter tails compared to the other two variables.

With the goal of recovering a causal order with EASE, first we estimate the $\widehat{\Psi}$ matrix from the full dataset, by setting the number of exceedances to $k=10$ (this corresponds to $\lfloor n^{0.3}\rfloor$, approximately). We run the EASE algorithm on the matrix $\widehat{\Psi}$, and we obtain the causal order $\widehat{\pi}^{-1}$ = (EURCHF, NOVN, ROG, NESN); this agrees with the proposed ground truth of Figure 4. As a comparison, also ICA-LiNGAM and Pairwise LiNGAM recover a causal order that agrees with our hypothesis.

Since our results are based on the $k=10$ upper order statistic, we assess the variability of the estimates $\widehat{\Psi}$ for different values of $k_{n}=\lfloor n^{\nu}\rfloor$, with $\nu\in[0.2,0.7]$. Figure 10 shows the estimated coefficients $\widehat{\Psi}$ for the pairs (EURCHF, NESN), (EURCHF, NOVN), and (EURCHF, ROGN), with the corresponding 90% bootstrap confidence intervals. In the three plots, the black (blue) line corresponds to the estimated coefficient $\widehat{\Psi}_{\text{EURCHF},i}$ ($\widehat{\Psi}_{i,\text{EURCHF}}$), with $i=$ NESN, NOVN, ROG. We can interpret the difference between the black and blue lines as a causal signal, since $\Psi_{ij}-\Psi_{ji}>0$ if variable $i$ causes variable $j$, for $i,j\in V$ (see Section 4.1). For the pairs (EURCHF, NESN) and (EURCHF, NOVN) the blue and the black lines overlap for all values of the upper order statistic $k$, and therefore any possible causal effect is not identified by the estimated coefficient $\widehat{\Psi}$. This result agrees with Example 2 of Section 4.3 which shows that the causal tail coefficients do not identify a causal signal when the ancestor has a heavier tail than its descendant — as is the case for the pairs (EURCHF, NESN) and (EURCHF, NOVN). In contrast, if we consider the pair (EURCHF, ROGN) we notice that the difference $\widehat{\Psi}_{\text{EURCHF},\text{ROGN}}-\widehat{\Psi}_{\text{ROGN},\text{EURCHF}}$ is positive for all fractional exponents $\nu\leq 0.4$. This can be explained by the fact that EURCHF and ROGN have comparable tail indices, and therefore it is easier for the coefficient $\Psi$ to detect a possible causal effect between the variables. In Appendix F, we show the dynamic evolution of the $\widehat{\Psi}$ coefficient across time.

Given the highly complex nature of financial markets, we do not take the conclusion of this experiment as a definite answer but rather consider it as an indication for a possible causal relationship in this data.

5.3 River data

We apply the EASE algorithm to the average daily discharges of the rivers located in the upper Danube basin. This dataset has been studied in Asadi et al. [1], Engelke and Hitz [12] and Mhalla et al. [31], and it is made available by the Bavarian Environmental Agency (http://www.gkd.bayern.de). We consider average daily discharges for 12 stations along the basin, representing the different tributaries and different sections of the Danube, and 11 of them are a subset of the 31 stations selected by Asadi et al. [1]. We exclude some of the 31 stations that are spatially very close since those are highly dependent and almost indistinguishable. For convenience, we name the stations with the same numbers used in Asadi et al. [1]. While Asadi et al. [1] decluster the data prior to their analysis in order to obtain independent samples, we use all observations despite the possible temporal dependence. In fact, for extreme value copulas, [54] show that the use of a larger but possibly dependent dataset can decrease the asymptotic estimation error. Moreover, Fawcett and Walshaw [15] argue that considering all exceedances over a high threshold reduces the bias of the maximum likelihood estimators compared to a declustered analysis. To account for the time dependence of the exceedances, they adjust the standard errors using methods presented by Smith [49]. In this experiment, we compute the standard errors according to the adjustment proposed by Fawcett and Walshaw [15].

The dataset spans from 1960 to 2009, where we consider only the summer months, i.e., June, July, and August. The rationale is that most of the extreme observations occur in summer due to heavy rainfall. The final dataset contains $n=4600$ observations. The rivers have an average volume that ranges between 20 $m^{3}/s$ (for the upstream rivers) and 1400 $m^{3}/s$ (for the downstream rivers). A map of the basin can be seen in Figure 13 in Appendix E.In order to implement our method, we first assess whether the equal tail index assumption is satisfied. To do so, we consider a regional model similar to the one presented by Asadi et al. [1]. We split the stations into four separate regions. Region 1 contains three stations in the southwest of the upper Danube basin and the catchment areas are located at mid-altitude; region 2 includes three stations in the Inn-Salzach basin whose tributaries are located in high-altitude alpine regions; region 3 contains four stations along the main Danube with large average water volume; region 4 comprises two stations in the north of the Danube. For each region, we fit a Poisson point process likelihood [7, Chap. 7] by considering exceedances over the 90% quantile and by constraining the shape parameter $\xi$ to be equal across the stations within the same region. To address the presence of temporal dependence in the exceedances, we adjust the standard errors as shown by Fawcett and Walshaw [15] and, based on these, we compute approximate confidence intervals. For each region, the estimated shape parameter and the corresponding confidence intervals are 0.167 (0.062, 0.273), 0.145 (0.047, 0.242), 0.133 (0.027, 0.239) and 0.229 (0.099, 0.358), respectively. The fact that the confidence intervals overlap suggests that our assumption of equal tail index across the variables is satisfied. Moreover, all confidence intervals do not include the zero value, and therefore the data can be deemed to be heavy-tailed.

We perform two separate analyses to identify causal structures both in space and time. Regarding the spatial structure, the goal is to recover the causal order of the network flow of the 12 stations on the rivers. This is a non-trivial task for two reasons. First, the water discharges at the stations can be confounded by rainfall that spreads out across the region of interest. Second, the large difference in water volume between the stations can further mask the causal structure of the water flow. The true DAG of the spatial disposition of the stations is shown in Figure 5. We run the EASE algorithm based on the $\Gamma$ coefficient, considering the upper tails, and setting the number of exceedances to $k=29=\lfloor n^{0.4}\rfloor$. The estimated causal order $\widehat{\pi}^{-1}=$ (23, 32, 26, 28, 19, 21, 11, 9, 7, 14, 13, 1) is correct, and the corresponding fully connected DAG has an SID equal to 0. We also run ICA-LiNGAM and Pairwise LiNGAM on the same dataset, and we obtain an SID of 0 and 0.053, respectively. Clearly, in this example, the causal structure in the bulk of the distribution is the same as in the extremes.

To assess the variability of our results, we compute the average SID of EASE over 50 bootstrap samples, for different values of the threshold parameter $k=\lfloor n^{\nu}\rfloor$, $\nu\in[0.2,0.7]$ — see Figure 9. From this figure, we observe that the fractional exponent $\nu\approx 0.4$ yields a good performance both in terms of SID and variability. The value for the optimal fractional exponent agrees with the empirical findings of Section 5.1.

Concerning the time-series analysis, we consider each station individually and try to recover the direction of time from the lagged data. Consider an autoregressive (AR) process of order $p\geq 1$,

[TABLE]

where the $\varepsilon_{j}$ are regularly varying with comparable tails. For a detailed discussion of such time series models we refer to Basrak and Segers [2] and Embrechts et al. [11, Chapter 7]. Peters et al. [40] prove that an AR($p$) process is time reversible, i.e., can be represented by an AR($p$) process in the reversed time direction, if and only if the noise is Gaussian. This means that for heavy-tailed random variables, one can in principle detect the direction of time from the data. For each station, we construct a dataset $D$ where the rows correspond to different days and the columns to different lags. We denote by $X_{0},X_{1},\dots,X_{6}$ the columns containing the current and lagged values of the station discharge. We then run EASE on the dataset $D$ and recover a causal order $\widehat{\pi}$ over the seven variables $X_{0},\dots,X_{6}$. We say that the direction of time is correctly inferred if the estimated causal order places the lags in the correct position, i.e., $\widehat{\pi}(j)<\widehat{\pi}(i)$ if $j<i$, for $i,j=0,\dots,6$. EASE successfully recovers the direction of time for 11 out of the 12 stations. As a comparison, ICA-LiNGAM, and Pairwise LiNGAM find the correct order for 9 and 11 stations, respectively.

So far, we have not considered a multivariate time series analysis of the dataset. On the one hand, in this particular application of the river data, the effects among the stations are almost instantaneous — due to the closeness of the water catchment areas, the fact that we have daily values and the river speed, which is about ten kilometres per hour. On the other hand, time information usually helps in estimating causal relations. For this reason, we apply multivariate Granger causality [20] to the river data, considering one day lag. For each pair of stations $(i,j)$, we say that $i$ Granger-causes $j$ if the corresponding $p$-value is significant at a 0.05 level, after the Bonferroni correction. We sort the significant $p$-values in ascending order, and we sequentially add the directed edge $(i,j)$ if nodes $i$ and $j$ are not connected. We continue until the skeleton of the resulting graph is connected, or all the $p$-values have been selected. The resulting directed tree achieves an SID of 0.083. Alternatively, if we sequentially add directed edges $(i,j)$ that do not create cycles (until all $p$-values have been selected) we obtain a DAG with an SID of 0.196.

6 Discussion and future work

In several real-world phenomena, the causal mechanisms between the variables appear more clearly during extreme events. Moreover, there are situations where the causal relationship in the bulk of the distribution differs from the structure in the tails. We have introduced an algorithm, named extremal ancestral search (EASE), that is shown to consistently recover the causal order of a DAG from extreme observations only. EASE has the advantage of relying on the pairwise causal tail coefficient, and therefore it is computationally efficient. In addition, our algorithm can deal with the presence of hidden confounders and performs well for small sample sizes and high dimensions. The EASE algorithm is robust to model misspecifications, such as nonlinear relationships in the bulk of the distribution, and strictly monotone increasing transformations applied marginally to each variable.

This work sheds light on a connection between causality and extremes, and thereby opens new directions of research. In particular, it might be interesting to study the properties of the causal tail coefficient under more general conditions. This includes high-dimensional settings where the dimension grows with the sample size, more general SCMs where the functional relations between the variables can be nonlinear, and settings where the noise variables have lighter tails. For example, in future research, one could combine EASE with regression techniques to obtain the complete DAG structure, compute the residuals and then test them for the assumption of common tail indices.

Another possible extension may consider multivariate time series data, where the temporal order of cause and effect could help to estimate causal relationships among variables, see, e.g., Granger [20]. Future work might study how to combine our approach with the Granger causality framework. For instance, one could first perform a Granger causality analysis, and then, apply EASE to the residuals of the vector autoregression model. For a careful study in this direction, it would also be necessary to investigate the statistical properties of the residuals.

Acknowledgements

We thank Cesare Miglioli, Stanislav Volgushev and Linbo Wang for helpful discussions. We are grateful to the editorial team and two anonymous referees for constructive comments that improved the paper. JP was supported by research grants from VILLUM FONDEN and the Carlsberg Foundation, and SE was supported by the Swiss National Science Foundation.

Appendix A Some facts about regular variation

In the sequel, for any two functions $f$, $g:\mathbb{R}\rightarrow\mathbb{R}$, we write $f\sim g$ if $f(x)/g(x)\to 1$ as $x\to\infty$. Also, we write $S_{p}:=Y_{1}+\dots+Y_{p}$, and $M_{p}:=\max(Y_{1},\dots,Y_{p})$.

Consider independent random variables $Y_{1},\dots,Y_{p}$ and assume that they have comparable upper tails, i.e., there exist $c_{j},\alpha>0$ and $\ell\in\mbox{RV}_{0}$ such that for all $j\in\{1,\dots,p\}$

[TABLE]

Lemma 3**.**

Let $Y_{1},\dots,Y_{p}$ be real-valued independent regularly varying random variables with comparable tails. Then,

[TABLE]

The proof for $p=2$ of Lemma 3 can be found in Feller [16, p. 278] and can be extended to a general $p$ using induction.

An important property of regularly varying random variable is the max-sum-equivalence presented in the following lemma [11, Sec. 1.3.1].

Lemma 4**.**

Let $Y_{1},\dots,Y_{p}$ be real-valued, independent regularly varying random variables with comparable tails. Then, as $x\to\infty$,

[TABLE]

Proof.

We can write, as $x\to\infty$,

[TABLE]

where $\Delta(x)$ contains terms of higher order interactions of the sets $\{Y_{j}>x\}$, $j=1,\dots,p$. Because of independence, the probability

[TABLE]

Similarly, this holds for the terms in $\Delta(x)$. Recalling, by Lemma 3, that $\mbox{P}(Y_{1}>x)+\dots+\mbox{P}(Y_{p}>x)\sim\mbox{P}(S_{p}>x)$, the result follows. ∎

Lemma 5**.**

Let $Y_{1},\dots,Y_{p}$ be real-valued independent regularly varying random variables with comparable tails. Then, for $j=1,\dots,p$,

[TABLE]

Proof.

For any $x>0$ and $\delta\in(0,1/(2p-2))$, we have, for $j=1,\dots,p$,

[TABLE]

Considering the upper bound, it holds, $j=1,\dots,p$,

[TABLE]

Regarding the lower bound, we get, as $x\to\infty$, $j=1,\dots,p$,

[TABLE]

Dividing everything by $\mbox{P}(Y_{j}>x)$, and letting first $x\to\infty$ and then $\delta\downarrow 0$, we get the desired result. ∎

Lemma 6**.**

Let $Y_{1},\dots,Y_{p}$ be real-valued independent regularly varying random variables with comparable tails. Then, as $x\to\infty$,

[TABLE]

Proof.

We first write $\mbox{P}\{S_{p}>x,M_{p}\leq x\}=\mbox{P}(S_{p}>x)-\mbox{P}\{S_{p}>x,M_{p}>x\}.$ Let $I:=\{1,\dots,p\}$. By definition of the maximum function and using the inclusion-exclusion principle, it follows that

[TABLE]

Regarding the summands in first term, it holds by Lemma 5, $h\in I$,

[TABLE]

The summands in the second term can be upper bounded by, $1\leq h<h^{\prime}\leq p$,

[TABLE]

The same holds for $\Delta(x)$ which contains terms of higher order interactions of the sets $\{Y_{j}>x\}$, $j\in I$. Putting everything together, we obtain

[TABLE]

where in the last equality we used Lemma 3. ∎

Appendix B Proofs

B.1 Proof of Lemma 1

Proof.

Let $j,k\in V$ and $j\neq k$. Recall that each variable $X_{h}$, $h\in V$, can be expressed as a weighted sum of the noise terms $\varepsilon_{1},\dots,\varepsilon_{p}$ belonging to the ancestors of $X_{h}$, as shown in (5). Therefore, we can write $X_{j}$ and $X_{k}$ as

[TABLE]

where $A_{jk}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)$, $A_{jk^{*}}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)^{c}$ and similarly for $A_{kj^{*}}$. We have

[TABLE]

The second summand can be bounded by

[TABLE]

by Lemma 6. For the first term, we use the inclusion-exclusion principle to write

[TABLE]

where $\Delta(x)$ contains terms of higher order interactions of the sets $\{\beta_{h\to j}\varepsilon_{h}>x\}$, $h\in\mathrm{An}(j,G)$. The probability

[TABLE]

The same holds for all finitely many terms in $\Delta(x)$. We further note that for all $h\in\mathrm{An}(j,G)$, by Lemma 5,

[TABLE]

Putting everything together, we can rewrite

[TABLE]

For $h\in A_{jk}$, let $c=\beta_{h\to j}/\beta_{h\to k}>0$, and note that for every $x>0$,

[TABLE]

Therefore, using Lemma 5 and that $F_{k}(cx)\to 1$ as $x\to\infty$, it follows that

[TABLE]

On the other hand, for $h\in A_{jk^{*}}$, we have that $X_{k}$ and $\varepsilon_{h}$ are independent, and therefore

[TABLE]

Consequently,

[TABLE]

where the second last equality follows from the fact that $\beta_{h\to j}\varepsilon_{h}$, $h\in A_{jk}$, are independent regularly varying random variables, see Lemma 3, and the last equality holds because we assume that the noise variables $\varepsilon_{j}$, $j\in V$, have comparable tails; see Section 2.1. ∎

B.2 Proof of Theorem 1

Proof.

Recall that $\mathrm{an}(j,G)=\mathrm{An}(j,G)\setminus\{j\}$ and define $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)$.

$(a).$

Suppose $X_{1}$ causes $X_{2}$, i.e., $1\in\mathrm{an}(2,G)$. This implies that $\mathrm{An}(1,G)\subset\mathrm{An}(2,G)$ and thus $A_{12}=\mathrm{An}(1,G)\subset\mathrm{An}(2,G)$. By applying Lemma 1 we obtain $\Gamma_{12}=1$ and $\Gamma_{21}\in(1/2,1)$.

Conversely, suppose that $\Gamma_{12}=1$ and $\Gamma_{21}\in(1/2,1)$. If $\Gamma_{12}=1$ then the numerator and denominator of the second term in Lemma 1 must be equal and strictly positive. This implies that $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)=\mathrm{An}(1,G)\neq\varnothing$. It follows that $\mathrm{An}(1,G)\subseteq\mathrm{An}(2,G)$, i.e., $1\in\mathrm{An}(2,G)$. At the same time, if $\Gamma_{21}\in(1/2,1)$, then the numerator of the second term in Lemma 1 must be positive and smaller than the denominator. This means that $A_{12}\neq\varnothing$ and $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)\subset\mathrm{An}(2,G)$. Thus, it follows that $\mathrm{An}(1,G)\subset\mathrm{An}(2,G)$. Therefore, $1\in\mathrm{an}(2,G)$, that is, $X_{1}$ causes $X_{2}$.

$(b).$

By symmetry, as case ($a$).

$(c).$

Suppose there is no causal link between $X_{1}$ and $X_{2}$, i.e., $\mathrm{An}(1,G)\cap\mathrm{An}(2,G)=\varnothing$. Then, $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)=\varnothing$ and by Lemma 1, we obtain $\Gamma_{12}=\Gamma_{21}=1/2$.

Suppose now that $\Gamma_{12}=\Gamma_{21}=1/2$. This means that the numerator of the second term in Lemma 1 must be equal to zero. This implies that $A_{12}=\varnothing$ and therefore $\mathrm{An}(1,G)\cap\mathrm{An}(2,G)=\varnothing$, that is, there is no causal link between $X_{1}$ and $X_{2}$.

$(d).$

Suppose there is a node $j\notin\{1,2\}$ such that $X_{j}$ is a common cause of $X_{1}$ and $X_{2}$, i.e., $j\in\mathrm{an}(1,G)$ and $j\in\mathrm{an}(2,G)$. Then $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)$ is non-empty. Since $\mathrm{An}(1,G)\neq\mathrm{An}(2,G)$, it follows that $A_{12}\subset\mathrm{An}(i,G)$, for $i=1,2$. Thus, according to Lemma 1 we have $\Gamma_{12},\Gamma_{21}\in(1/2,1)$.

Conversely, suppose that $\Gamma_{12},\Gamma_{21}\in(1/2,1)$. If $\Gamma_{12}\in(1/2,1)$, then the numerator of the second term in Lemma 1 must be positive and smaller than the denominator. This implies that $A_{12}\neq\varnothing$ and $A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)\subset\mathrm{An}(1,G)$. Similarly, if $\Gamma_{21}\in(1/2,1)$, it follows that $A_{21}=A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)\subset\mathrm{An}(2,G)$. This implies that $\mathrm{An}(1,G)\neq\mathrm{An}(2,G)$ and they are not disjoint. Therefore, there exists a node $j\notin\{1,2\}$ such that $j\in\mathrm{an}(1,G)$ and $j\in\mathrm{an}(2,G)$, i.e., $X_{j}$ is a common cause of $X_{1}$ and $X_{2}$.

∎

B.3 Proof of Theorem 2

Proof.

For simplicity we will write $k=k_{n}$ in the sequel. We only show the result for $\widehat{\Gamma}_{21}$, the proof for $\widehat{\Gamma}_{12}$ follows by symmetry. Recall that each variable $X_{h}$, $h\in V$, can be expressed as a weighted sum of the noise terms $\varepsilon_{1},\dots,\varepsilon_{p}$ belonging to the ancestors of $X_{h}$, as shown in (5). Therefore, we can write $X_{1}$ and $X_{2}$ as follows,

[TABLE]

where $A=A_{12}=\mathrm{An}(1,G)\cap\mathrm{An}(2,G)$ and $A_{j}=A_{jk^{*}}=\mathrm{An}(j,G)\setminus A$, for $j,k=1,2$. Thus, the estimator $\widehat{\Gamma}_{21}$ can be rewritten as

[TABLE]

Define the theoretical quantile function as

[TABLE]

Recall that $X_{(n-k),2}=\widehat{F}_{2}^{\leftarrow}(1-k/n)$ is the $(n-k)$-th order statistic of $X_{12},\dots,X_{n2}$, and as such an approximation to the theoretical quantile $U_{2}(n/k)$. Under the von Mises’ condition (9) this convergence of the intermediate order statistics can be made rigorous [10, Theorem 2.2.1], namely

[TABLE]

where $\alpha>0$ is the tail index of the variables in the SCM. This implies in particular that $X_{(n-k),2}\to\infty$, $n\to\infty$. For any $\delta_{1}>0$ define the event

[TABLE]

and note that by (18) it holds that $\mbox{P}(B_{n\delta_{1}})\to 1$ as $n\to\infty$.

Since the noise terms $\varepsilon_{h}$ are independent regularly varying random variables with comparable tails, then, by Lemma 3, we have

[TABLE]

Furthermore, from Resnick [42, Prop. 0.8] it holds, for all $h\in\mathrm{An}(2,G)$ and $x>0$,

[TABLE]

where $U_{2}$ is defined as in (17) for the distribution function $F_{2}$.

We treat the two terms in (16) separately. We can upper bound the absolute value of the second term, for any $\tau,\delta_{1}>0$, by

[TABLE]

where $\mbox{P}(B_{n\delta_{1}}^{c})\to 0$ as $n\to\infty$ by (18). The limit superior of the second term, as $n\to\infty$, can be bounded with Markov’s inequality by

[TABLE]

The last equality holds because $X_{2}$ is regularly varying with index $\alpha$ and because, by Lemma 6, we have

[TABLE]

Since $\delta_{1},\tau>0$ are arbitrary, it follows that $S_{2,n}=o_{P}(1)$.

For the first term, we use the inclusion-exclusion principle to write

[TABLE]

where $\Delta_{i}(X_{(n-k),2})$ contains terms of higher order interactions of the sets $\{\beta_{h\to 2}\varepsilon_{ih}>X_{(n-k),2}\}$, $h\in\mathrm{An}(2,G)$, $i=1,\dots,n$. First, we show that the terms $T_{2,n}$ and $T_{3,n}$ are $o_{P}(1)$. Considering $T_{2,n}$, for each $h,h^{\prime}\in\mathrm{An}(2,G)$, $h<h^{\prime}$, define

[TABLE]

We can upper bound its absolute value, for any $\tau,\delta_{1}>0$, by

[TABLE]

By Markov’s inequality,

[TABLE]

where in the last line we used property (20) and the fact that $k/n\to 0$ as $n\to\infty$, and $p_{h2},p_{2h^{\prime}}$ are defined as

[TABLE]

Since $\delta_{1},\tau$ are arbitrary, putting together the finitely many terms $h<h^{\prime}$ where $h,h^{\prime}\in\mathrm{An}(2,G)$, it follows that $T_{2,n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}0$. Using a similar argument as the one for $T_{2,n}$, one can show that $T_{3,n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}0$.

We want to show that $T_{1,n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}\Gamma_{21}$. Rewrite

[TABLE]

Using an argument similar to the one for $S_{2,n}$ on page B.3, one can show that $U_{3,n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}0$.

Regarding $U_{1,n}$, for each $h\in A$, define

[TABLE]

Note that, for $h\in A$, we have that both $\beta_{h\to 2}>0$ and $\beta_{h\to 1}>0$, therefore we can bound

[TABLE]

where $c=\beta_{h\to 1}/\beta_{h\to 2}>0$ and

[TABLE]

We show first that $V_{n}$ converges in probability to $p_{2h}$, defined in (21). This is motivated by the fact that $V_{n}$ is the empirical version of

[TABLE]

where the limit follows from property (20). We will study the asymptotic properties of $V_{n}$. For $x>0$ define

[TABLE]

which is a nonincreasing function of $x$. Observe that on the set $B_{n\delta_{1}}$ defined in (19), we may bound the random variable $V_{n}$ by

[TABLE]

We further compute the limits as $n\to\infty$

[TABLE]

An application of Chebyshev’s inequality yields

[TABLE]

For some $\tau>0$, choose $\delta_{1}>0$ such that $p_{2h}(1+\delta_{1})^{-\alpha}>p_{2h}-\tau$ and $p_{2h}(1-\delta_{1})^{-\alpha}<p_{2h}+\tau$. Then with (24), (25) and the fact that $\mbox{P}(B_{n\delta_{1}}^{c})\to 0$ for $B_{n\delta_{1}}$ in (19), we conclude

[TABLE]

that is, $V_{n}$ converges in probability to $p_{2h}$, as $n\to\infty$.

Furthermore, using (22) we will now show

[TABLE]

Indeed, for any $\tau,\delta_{1}>0$, as $n\to\infty$

[TABLE]

since $\widehat{F}_{1}(x)$ converges in probability to $F_{1}(x)$ for all $x\in\mathbb{R}$, and $U_{2}(n/k)\to\infty$, as $n\to\infty$. Moreover, with a similar argument as for $S_{2,n}$, one can show that $W_{n}$ defined in (23) converges in probability to 0, as $n\to\infty$.

Putting everything together, using (22), (26) and (27) we conclude that $U_{1,n}^{(h)}\stackrel{{\scriptstyle P}}{{\longrightarrow}}p_{2h}$, $h\in A$, and thus

[TABLE]

Considering the term $U_{2,n}$, for each $h\in A_{2}$, define

[TABLE]

On the event $B_{n\delta_{1}}$, we can bound $U_{2,n}^{(h)}$ by

[TABLE]

where we let, for all $x>0$,

[TABLE]

Since $X_{1}$ is independent of $\varepsilon_{h}$, for $h\in A_{2}$, the values in the sum can be seen as $M(x)=M_{n}(x)$ random samples out of $\{1/n,\dots,1\}$ without replacement, where $M(x)$ is Binomial with success probability

[TABLE]

Let $Z_{ni}$ be random samples out of $\{1/n,\dots,1\}$ without replacement, for all $i=1,\dots,n$, $n\in\mathbb{N}$. Then $u_{2,n}(x)$ has the same distribution as

[TABLE]

By a similar argument as in (25) the distribution of $M(x)$ satisfies for any fixed $x\in(0,\infty)$

[TABLE]

where $m(x)=m_{n}(x)=\lceil{kp_{2h}x^{-\alpha}}\rceil$. Thus, for any $\delta_{2}>0$ and any $x>0$, the probability of the event

[TABLE]

converges to 1 as $n\to\infty$. Consider the quantity

[TABLE]

Theorem 5.1 in [45] states that the limit in probability of this sum of samples without replacement is the same as the corresponding sum of samples with replacement. Therefore, for any $x\in(0,\infty)$, we have the convergence in probability

[TABLE]

For $\tau>0$, choose $\delta_{1},\delta_{2}>0$ small enough such that

[TABLE]

Then we can bound the probability

[TABLE]

since $\mbox{P}(B_{n\delta_{1}}^{c})$ and $\mbox{P}(C_{n\delta_{2},1-\delta_{1}}^{c})$ converge to 0 as $n\to\infty$, and the last term converges to [math] as a consequence of (31) and (32). Similarly, we can show that

[TABLE]

and since $\tau>0$ is arbitrary, $U_{2,n}^{(h)}\stackrel{{\scriptstyle P}}{{\longrightarrow}}p_{2h}/2$, for $h\in A_{2}$. Therefore,

[TABLE]

and $\widehat{\Gamma}_{21}\stackrel{{\scriptstyle P}}{{\longrightarrow}}\Gamma_{21}$.

∎

B.4 Proof of Proposition 1

Proof.

Let $s\in S:=\{1,\dots,p\}$ and denote by $i_{s}\in V$ the node chosen by the algorithm at step $s$. Denote by

[TABLE]

the set of nodes chosen by the algorithm before step $s$. Let $G_{s}=(V_{s},E_{s})$ be the subgraph of $G$ obtained by removing the nodes $H_{s}$ that are already chosen, i.e., $V_{s}=V\setminus H_{s}$ and $E_{s}=E\cap(V_{s}\times V_{s})$. Furthermore, define the score minimized by the algorithm to choose the node at step $s$,

[TABLE]

We want to show that EASE is a procedure that, for all $s\in S$, satisfies the statement

[TABLE]

We use strong induction. Namely, we prove that if $\Xi(s^{\prime})$ holds for all natural numbers $s^{\prime}<s$, then $\Xi(s)$ holds, too.

Fix $s\in S$ and suppose that for all $s^{\prime}\in S$, with $s^{\prime}<s$, $\Xi(s^{\prime})$ holds. Assume $\varphi:=i_{s}\in\arg\min_{i\in V_{s}}M_{i}^{(s)}$ and $\mathrm{an}(\varphi,G)\not\subseteq H_{s}$. Then, there exists a node $j\in V_{s}$ such that $j\in\mathrm{an}(\varphi,G)$, and by Theorem 1, $\Gamma_{j\varphi}=1$. It follows $M_{\varphi}^{(s)}=1$.

Also, since $G_{s}$ is a DAG, there exists a node $\ell\in V_{s}$ such that $\mathrm{an}(\ell,G_{s})=\varnothing$. If $\Xi(s^{\prime})$ holds for every natural number $s^{\prime}<s$, then $\mathrm{an}(\ell,G)\subseteq H_{s}$. Suppose not. Then, there exists a node $j\in V_{s}$ such that $j\in\mathrm{an}(\ell,G)$. Note that since $j\in\mathrm{an}(\ell,G)$ and $j\notin\mathrm{an}(\ell,G_{s})$, there exists a directed path from $j$ to $\ell$ in $G$ that is absent in $G_{s}$. Thus, there exists a node $h\in H_{s}$ that lies on such path, and it follows that $j\in\mathrm{an}(h,G)$, which is a contradiction. Since $\mathrm{an}(\ell,G)\subseteq H_{s}$, by Theorem 1 it holds $M_{\ell}^{(s)}<1=M_{\varphi}^{(s)}$, which is a contradiction.

Since $s\in S$ was arbitrary, we have proved that $\Xi(s)$ holds for all $s\in S$. Furthermore, note that $\Xi(1)$ holds as a special case of the argument above. Hence, we conclude that, for all $s^{\prime},s\in S$,

[TABLE]

and therefore $\pi(i_{s})=s$ is a causal order of $G$. ∎

B.5 Proof of Proposition 2

Proof.

If $\widehat{\pi}\not\in\Pi_{G}$, then there exists a node $i\in V$ that is chosen before one of its ancestors $u\in\mathrm{an}(i,G)$, i.e., $\widehat{\pi}(i)<\widehat{\pi}(u)$. Therefore, there exists a non-empty set $\tilde{V}\subseteq V$, with $i,u\in\tilde{V}$, such that

[TABLE]

Furthermore, since $G$ is a DAG, there exists a node $j\in\tilde{V}$ with no ancestors in $\tilde{V}$. Let $v\in\tilde{V}\setminus\{j\}$, and note that $v\notin\mathrm{An}(j,G)$. Thus, by (35), it follows that $\widehat{\Gamma}_{vj}-\widehat{\Gamma}_{ui}\geq 0$. Define $\widehat{\Delta}_{ij}:=|\widehat{\Gamma}_{ij}-\Gamma_{ij}|$ and note that the event $\widehat{\Gamma}_{vj}-\widehat{\Gamma}_{ui}\geq 0$ can be bounded by

[TABLE]

where $\Gamma_{ui}=1$ because $u\in\mathrm{an}(i,G)$, and $\Gamma_{vj}\leq\eta<1$. Therefore,

[TABLE]

This completes the proof of Proposition 2. ∎

B.6 Proof of Lemma 2

Proof.

Recall that each variable $X_{h}$, $h\in V$, can be expressed as a weighted sum of the noise terms $\varepsilon_{1},\dots,\varepsilon_{p}$ belonging to the ancestors of $X_{h}$, as shown in (5). Therefore, we can write $X_{j}$ and $X_{k}$ as follows,

[TABLE]

where $A_{jk}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)$, $A_{jk^{*}}=\mathrm{An}(j,G)\cap\mathrm{An}(k,G)^{c}$ and similarly for $A_{kj^{*}}$.

We treat the two terms of (11) separately. Consider the first term,

[TABLE]

Since $X_{j}$, $j\in V$, are regularly varying with index $\alpha>0$, using similar arguments as in Lemma 1, we can write

[TABLE]

Recall that $\sigma$ is defined as

[TABLE]

and that $0\leq\sigma(F_{k}(X_{k}))<1$, $k\in V$.

Consider the summands in (36) where $h\in A_{jk}$. We distinguish two cases. When the ratio $c=\beta_{h\to k}/\beta_{h\to j}>0$, we can bound the summand

[TABLE]

The last equality follows from Lemma 5 and since $\sigma(F_{k}(cx))\to 1$ as $x\to\infty$. When the ratio $c=\beta_{h\to k}/\beta_{h\to j}<0$, we obtain

[TABLE]

where the third inequality holds because, as $x\to\infty$, $F_{k}(X_{k})<F_{k}(cx)$ implies $\sigma(F_{k}(X_{k}))>\sigma(F_{k}(cx))$.

On the other hand, for all the summands in (36) where $h\in A_{jk^{*}}$ we have that $X_{k}$ and $\varepsilon_{h}$ are independent. Therefore,

[TABLE]

Consequently,

[TABLE]

Similarly, the second term can be shown to be

[TABLE]

Putting everything together yields the desired form of $\Psi_{jk}=\Psi_{jk}^{+}+\Psi_{jk}^{-}$. ∎

Appendix C Example of the EASE algorithm

Consider the DAG $G$ in Figure 6, with vertex set $V=\{1,2,3,4\}$. The set of causal orders of $G$ is $\Pi_{G}=\{(2,1,4,3),(2,1,3,4),(3,1,2,4)\}$. In Figure 7, we display the state-space tree of the EASE algorithm, i.e., the set of states that the algorithm can visit to find a causal order $\pi$ of $G$. Each state represents the status of the vector $\pi^{-1}$ during the algorithm evaluation. A state is red if all the paths below it lead to wrong causal order. The green states represent the causal orders of $G$.

Appendix D Experimental settings for the simulation study

The parameters of the simulation are the following.

• –

Distribution: Student’s $t$, with degrees of freedom $\alpha\in\{1.5,2.5,3.5\}$.

• –

Number of observations: $n\in\{500,1000,10000\}$.

• –

Number of variables: $p\in\{4,7,10,15,20,30,50\}$.

The settings that we consider are,

• –

Linear SCM,

• –

Linear SCM with hidden confounders,

• –

Nonlinear SCM,

• –

Linear SCM and uniform transformation of each variable.

For each combination of $n$, $p$, and $\alpha$ and each setting, we generate $n_{exp}=50$ random SCMs. Each SCM is built as follows.

Generate a random DAG.

1. (a)

Take a random permutation $\pi$ of the nodes $V=\{1,\dots,p\}$ that defines the causal order. 2. (b)

For each $i\in V$ such that $\pi(i)>1$, sample the number of parents

[TABLE]

from a binomial distribution. We set $q=\min\{5/(p-1),1/2\}$ so that, on average, there are 2.5 edges per node, when $p>10$. 3. (c)

Sample without replacement $n_{\mathrm{pa}}$ from $\{j\in V:\pi(j)<\pi(i)\}$ and name the resulting set $\mathrm{pa}(i)$. 2. 2.

Sample uniformly from $\{-0.9,-0.1\}\cup\{0.1,0.9\}$ the coefficients $\beta_{ij}$, where $i\in V$ and $j\in\mathrm{pa}(i,G)$. 3. 3.

In the case of hidden confounders,

1. (a)

Sample the number of confounding variables,

[TABLE]

from a binomial distribution. We set $q=2/(3p-3)$ so that, on average, there is one hidden confounder for every three nodes. 2. (b)

Sample without replacement $n_{\text{conf}}$ unordered pairs from $\left\{\{i,j\}:i,j\in V\right\}$ and name the resulting set $C$. 3. (c)

Update the parents of each node $i$ as $\mathrm{pa}(i):=\mathrm{pa}(i)\cup C_{i}$, where $C_{i}\subseteq C$ is the set of hidden confounders affecting node $i\in V$. Similarly, for each hidden confounder $c\in C$, set $\mathrm{pa}(c)=\varnothing$. 4. (d)

Sample uniformly from $\{-0.9,-0.1\}\cup\{0.1,0.9\}$ the coefficients $\beta_{ic}$, where $i\in V$ and $c\in C_{i}$. 4. 4.

Let $\widetilde{V}=V\cup C$111If there are no hidden confounders, then $C=\varnothing$.. For all $i\in\widetilde{V}$, sample $n$ i.i.d. copies of noise $\varepsilon_{i}\sim\text{Student's }t\ \text{with df}=\alpha$. 5. 5.

For each node $i\in\widetilde{V}$, generate

[TABLE]

1. (a)

in case of linear SCM, $f(X_{j})=X_{j}$, 2. (b)

in case of nonlinear SCM, $f(X_{j})=X_{j}\ \mathbf{1}\{\widehat{F}_{j}(X_{j})>0.95\}$, where $\widehat{F}_{j}$ is the empirical cdf of $X_{j}$. 6. 6.

In case of uniform margins, transform each variable by $X_{i}:=\widehat{F}_{i}(X_{i})$, where $\widehat{F}_{i}$ is the empirical cdf of $X_{i}$, $i\in\widetilde{V}$.

Appendix E Additional Figures and Tables

Appendix F Financial application

Starting from the financial application of Section 5.2, we study the dynamic evolution of the $\widehat{\Psi}$ coefficient across time. By looking at the time series of the EURCHF in Figure 16, we notice two periods of extremely high positive and negative returns, namely, August 2011 and January 2015. The second event was due to an unexpected decision of the Swiss National Bank (SNB) to remove the peg of 1.20 Swiss francs per Euro.

The idea is to estimate the $\widehat{\Psi}$ coefficients between the EURCHF and the three stocks on a rolling window of 1500 days. In this case, we use a threshold $k=10$ which corresponds approximately to a fractional exponent $\nu=0.3$, where we define $k=\lfloor n^{\nu}\rfloor$ and $n$ is the number of observations in the sample. In Figure 17, we notice that during turmoil periods, highlighted in red, the $\widehat{\Psi}$ coefficient is higher in the direction that goes from the EURCHF to the stocks. This is an example where the causal structure becomes clearer during extreme events. It is also interesting to observe that the $\widehat{\Psi}$ coefficient is quite stable during calmer market periods. For example, the stable black lines between 2012 and 2015 correspond to the currency peg maintained by the SNB. As a last note, the drop in the black lines in mid 2017 is due to the fact that the extreme event of 2011 ‘exits’ the rolling window.