<figure><img src="image_1.png" /></figure>

HAL
open science

# Causal Discovery from Time Series with Hybrids of Constraint-Based and Noise-Based Algorithms

Daria Bystrova, Charles Assaad, Julyan Arbel, Emilie Devijver, Éric Gaussier, Wilfried Thuiller

► To cite this version:

Daria Bystrova, Charles Assaad, Julyan Arbel, Emilie Devijver, Éric Gaussier, et al.. Causal Discovery from Time Series with Hybrids of Constraint-Based and Noise-Based Algorithms. Transactions on Machine Learning Research Journal, 2024. hal-04606168

HAL Id: hal-04606168

https://hal.science/hal-04606168

Submitted on 9 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Causal Discovery from Time Series with Hybrids of Constraint-Based and Noise-Based Algorithms

Daria Bystrova, Charles Assaad, Julyan Arbel, Emilie Devijver, Éric Gaussier, Wilfried Thuiller

► To cite this version:

Daria Bystrova, Charles Assaad, Julyan Arbel, Emilie Devijver, Éric Gaussier, et al.. Causal Discovery from Time Series with Hybrids of Constraint-Based and Noise-Based Algorithms. Transactions on Machine Learning Research Journal, 2024, pp.1-30. hal-04606168

HAL Id: hal-04606168

https://hal.science/hal-04606168v1

Submitted on 9 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CAUSAL DISCOVERY FROM TIME SERIES WITH HYBRIDS OF CONSTRAINT-BASED AND NOISE-BASED ALGORITHMS

A PREPRINT

<figure><img src="image_2.png" /></figure>

Daria Bystrova*
Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LECA, LJK, LIG, F38000, Grenoble, France

<figure><img src="image_3.png" /></figure>

Charles K. Assaad*
Sorbonne Université, INSERM, Institut Pierre Louis d'Epidémiologie et de Santé Publique, F75012, Paris, France

<figure><img src="image_4.png" /></figure>

Julyan Arbel
Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, F38000, Grenoble, France

<figure><img src="image_5.png" /></figure>

Emilie Devijver
Univ Grenoble Alpes, CNRS, Grenoble INP, LIG, F38000, Grenoble, France

<figure><img src="image_6.png" /></figure>

Eric Gaussier
Univ Grenoble Alpes, CNRS, Grenoble INP, LIG, F38000, Grenoble, France

<figure><img src="image_7.png" /></figure>

Wilfried Thuiller
Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, F38000, Grenoble, France

## ABSTRACT

Constraint-based methods and noise-based methods are two distinct families of methods proposed for uncovering causal graphs from observational data. However, both operate under strong assumptions that may be challenging to validate or could be violated in real-world scenarios. In response to these challenges, there is a growing interest in hybrid methods that amalgamate principles from both methods, showing robustness to assumption violations. This paper introduces a novel comprehensive framework for hybridizing constraint-based and noise-based methods designed to uncover causal graphs from observational time series. The framework is structured into two classes. The first class employs a noise-based strategy to identify a super graph, containing the true graph, followed by a constraint-based strategy to eliminate unnecessary edges. In the second class, a constraint-based strategy is applied to identify a skeleton, which is then oriented using a noise-based strategy. The paper provides theoretical guarantees for each class under the condition that all assumptions are satisfied, and it outlines some properties when assumptions are violated. To validate the efficacy of the framework, two algorithms from each class are experimentally tested on simulated data, realistic ecological data, and real datasets sourced from diverse applications. Notably, two novel datasets related to Information Technology monitoring are introduced within the set of considered real datasets. The experimental results underscore the robustness and effectiveness of the hybrid approaches across a broad spectrum of datasets.

**Keywords** Causal discovery · Time series · Noise-based · Constraint-based

## 1 Introduction

Recent technological advances allow collecting observational time series on complex dynamical systems in various fields, such as biodiversity monitoring in ecology [Dornelas et al., 2018], epidemiology [Meci et al., 2022, Arlegui

<footnote>*These authors contributed equally to this work.</footnote>
CAUSAL DISCOVERY FROM TIME SERIES WITH
HYBRIDS OF CONSTRAINT-BASED AND NOISE-BASED
ALGORITHMS

A PREPRINT

<figure><img src="image_2.png" /></figure>

Daria Bystrova*
Univ. Grenoble Alpes, CNRS, Inria,
Grenoble INP, LECA, LJK, LIG,
F38000, Grenoble, France

<figure><img src="image_3.png" /></figure>

Charles K. Assaad*
Sorbonne Université, INSERM,
Institut Pierre Louis d'Epidémiologie et de Santé Publique,
F75012, Paris, France

<figure><img src="image_4.png" /></figure>

Julyan Arbel
Univ. Grenoble Alpes, Inria, CNRS,
Grenoble INP, LJK,
F38000, Grenoble, France

<figure><img src="image_5.png" /></figure>

Emilie Devijver
Univ Grenoble Alpes, CNRS,
Grenoble INP, LIG,
F38000, Grenoble, France

<figure><img src="image_6.png" /></figure>

Eric Gaussier
Univ Grenoble Alpes, CNRS,
Grenoble INP, LIG,
F38000, Grenoble, France

<figure><img src="image_7.png" /></figure>

Wilfried Thuiller
Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS,
LECA,
F38000, Grenoble, France

ABSTRACT

Constraint-based methods and noise-based methods are two distinct families of methods proposed for uncovering causal graphs from observational data. However, both operate under strong assumptions that may be challenging to validate or could be violated in real-world scenarios. In response to these challenges, there is a growing interest in hybrid methods that amalgamate principles from both methods, showing robustness to assumption violations. This paper introduces a novel comprehensive framework for hybridizing constraint-based and noise-based methods designed to uncover causal graphs from observational time series. The framework is structured into two classes. The first class employs a noise-based strategy to identify a super graph, containing the true graph, followed by a constraint-based strategy to eliminate unnecessary edges. In the second class, a constraint-based strategy is applied to identify a skeleton, which is then oriented using a noise-based strategy. The paper provides theoretical guarantees for each class under the condition that all assumptions are satisfied, and it outlines some properties when assumptions are violated. To validate the efficacy of the framework, two algorithms from each class are experimentally tested on simulated data, realistic ecological data, and real datasets sourced from diverse applications. Notably, two novel datasets related to Information Technology monitoring are introduced within the set of considered real datasets. The experimental results underscore the robustness and effectiveness of the hybrid approaches across a broad spectrum of datasets.

**Keywords** Causal discovery · Time series · Noise-based · Constraint-based

# 1 Introduction

Recent technological advances allow collecting observational time series on complex dynamical systems in various fields, such as biodiversity monitoring in ecology [Dornelas et al., 2018], epidemiology [Meci et al., 2022, Arlegui

<footnote>*These authors contributed equally to this work.</footnote>
et al., 2023, Bales et al., 2023, Moreau et al., 2023], healthcare [Morid et al., 2023] and Information Technology (IT) monitoring systems [Tamburri et al., 2020, Assaad et al., 2023]. One of the key objectives in studying such dynamical systems is to understand the causal relationships between the system's components. To find these causal relations, experts can employ causal discovery methods for time series, which aim to build a causal graph from observational data. These methods can be categorized into several families, including Granger-based [Granger, 1969], constraint-based [Spirtes et al., 2000, Runge, 2020, Assaad et al., 2022c], score-based [Chickering, 2002b], continuous optimization-based [Zheng et al., 2018] and noise-based families [Hyvärinen et al., 2008, Peters et al., 2013] —for more details see Assaad et al. [2022a], Hasan et al. [2023], Gong et al. [2023]. Each family has its own set of assumptions, which may or may not be suitable for a specific dataset. Therefore, no single method stands out in all situations [Assaad et al., 2022a].

Hybrid frameworks combine several methods from different families to enhance graph inference [Hasan et al., 2023]. For non-temporal data, several authors propose to combine ideas from constraint-based and score-based methods to improve scalability [Tsamardinos et al., 2006] or robustness to small sample size [Ogarrio et al., 2016]. For temporal data, SVAR-GFCI, proposed by Malinsky and Spirtes [2018], is a time-series generalization of the hybrid method GFCI which is based on the score-based and constraint-based algorithms.

Another type of hybrid framework, which is our main focus in this work, is based on the combination of constraint-based and noise-based families. The advantage of constraint-based methods is that they are non-parametric (i.e., no assumption is made on the form of the underlying causal relationships), while the limitation is that they require strong non-testable assumptions and can only recover the causal graph up to its Markov equivalence class, i.e., orientation of some edges could be unknown in the inferred graph as several graphs represent the same conditional dependence structure. On the other hand, noise-based methods are capable of recovering true graphs. So, by combining methods from both families, although we require assumptions of both families of methods, some assumptions can be weakened and we can recover the true causal graph. There exist several methods of this type for non-temporal data, such as PClingam [Hoyer et al., 2008] or FRITL [Chen et al., 2021]. Assaad et al. [2021] introduced NBCB<sup>acyclic</sup> 2 for temporal data assuming that there are no cyclic causal relations between different time series.

This paper presents a hybrid framework for temporal data using noise-based and constraint-based algorithms. In this framework, we consider two different classes of methods, which we denote NBCB and CBNB. Both classes can infer different types of causal graphs that differentiate between instantaneous relations and lagged relations. To construct these types of causal graphs, NBCB and CBNB orient edges using a noise-based strategy and prune edges using a constraint-based strategy. The main difference between NBCB and CBNB is that the former starts by orienting the graph and then proceeds to pruning, while the latter starts by pruning and then proceeds to orientation. Most importantly, NBCB and CBNB combine the parts from corresponding methods in an efficient way, such that information in the first part improves the efficiency of the second part. At the core of CBNB lies the notion of an undirected cycle group, which we introduce to optimize the search for orientation. The advantage of the proposed hybrid framework is that, in practice, it has a trade-off performance between the constraint-based and noise-based algorithms. In comparison to the constrained-based family, the proposed methods do not require the so-called faithfulness assumption (but require a weaker assumption named adjacency faithfulness) and can recover the true causal graph instead of restricting only to the Markov equivalence class. In comparison to noise-based methods, the proposed methods provide a better pruning of edges when the sample size is small [Malinsky and Danks, 2018]. In this paper, we also show the distinct responses of NBCB and CBNB when assumptions are violated. This demonstration highlights the expected results across various scenarios and aids in determining the most suitable methodological approach for a given problem.

In summary, our main contributions are the following:

* We propose a hybrid framework for the causal discovery of time series that combines parts of noise-based and constraint-based algorithms. Within this framework, we derive two classes of algorithms, NBCB and CBNB, which we optimize to infer the causal graph from time series.

* We study theoretically to which extent each class of algorithms is robust against assumption violation.

* We provide extensive simulation studies and real data applications to illustrate the applicability of our approach and their enhanced capabilities against assumption violation compared to original methods.

* We introduce two novel datasets about IT monitoring within the set of considered real datasets.

The remainder of the paper is organized as follows: Section 2 describes the different types of causal graphs that can be used to represent causal relations between time series and the different assumptions related to those graphs. Section 3 discusses related work and particularly details the steps that compose noise-based and constraint-based algorithms.

<footnote n="2">²In Assaad et al. [2021], this method was called NBCB but in this work, we denote it as NBCB<sup>acyclic</sup> to explicitly point out that it assumes that the summary causal graph is acyclic. For more details about summary causal graphs see Section 2.</footnote>
et al., 2023, Bales et al., 2023, Moreau et al., 2023], healthcare [Morid et al., 2023] and Information Technology (IT) monitoring systems [Tamburri et al., 2020, Assaad et al., 2023]. One of the key objectives in studying such dynamical systems is to understand the causal relationships between the system's components. To find these causal relations, experts can employ causal discovery methods for time series, which aim to build a causal graph from observational data. These methods can be categorized into several families, including Granger-based [Granger, 1969], constraint-based [Spirtes et al., 2000, Runge, 2020, Assaad et al., 2022c], score-based [Chickering, 2002b], continuous optimization-based [Zheng et al., 2018] and noise-based families [Hyvärinen et al., 2008, Peters et al., 2013] —for more details see Assaad et al. [2022a], Hasan et al. [2023], Gong et al. [2023]. Each family has its own set of assumptions, which may or may not be suitable for a specific dataset. Therefore, no single method stands out in all situations [Assaad et al., 2022a].

Hybrid frameworks combine several methods from different families to enhance graph inference [Hasan et al., 2023]. For non-temporal data, several authors propose to combine ideas from constraint-based and score-based methods to improve scalability [Tsamardinos et al., 2006] or robustness to small sample size [Ogarrio et al., 2016]. For temporal data, SVAR-GFCI, proposed by Malinsky and Spirtes [2018], is a time-series generalization of the hybrid method GFCI which is based on the score-based and constraint-based algorithms.

Another type of hybrid framework, which is our main focus in this work, is based on the combination of constraint-based and noise-based families. The advantage of constraint-based methods is that they are non-parametric (i.e., no assumption is made on the form of the underlying causal relationships), while the limitation is that they require strong non-testable assumptions and can only recover the causal graph up to its Markov equivalence class, i.e., orientation of some edges could be unknown in the inferred graph as several graphs represent the same conditional dependence structure. On the other hand, noise-based methods are capable of recovering true graphs. So, by combining methods from both families, although we require assumptions of both families of methods, some assumptions can be weakened and we can recover the true causal graph. There exist several methods of this type for non-temporal data, such as PClingam [Hoyer et al., 2008] or FRITL [Chen et al., 2021]. Assaad et al. [2021] introduced NBCB<sup>acyclic</sup> 2 for temporal data assuming that there are no cyclic causal relations between different time series.

This paper presents a hybrid framework for temporal data using noise-based and constraint-based algorithms. In this framework, we consider two different classes of methods, which we denote NBCB and CBNB. Both classes can infer different types of causal graphs that differentiate between instantaneous relations and lagged relations. To construct these types of causal graphs, NBCB and CBNB orient edges using a noise-based strategy and prune edges using a constraint-based strategy. The main difference between NBCB and CBNB is that the former starts by orienting the graph and then proceeds to pruning, while the latter starts by pruning and then proceeds to orientation. Most, importantly, NBCB and CBNB combine the parts from corresponding methods in an efficient way, such that information in the first part improves the efficiency of the second part. At the core of CBNB lies the notion of an undirected cycle group, which we introduce to optimize the search for orientation. The advantage of the proposed hybrid framework is that, in practice, it has a trade-off performance between the constraint-based and noise-based algorithms. In comparison to the constrained-based family, the proposed methods do not require the so-called faithfulness assumption (but require a weaker assumption named adjacency faithfulness) and can recover the true causal graph instead of restricting only to the Markov equivalence class. In comparison to noise-based methods, the proposed methods provide a better pruning of edges when the sample size is small [Malinsky and Danks, 2018]. In this paper, we also show the distinct responses of NBCB and CBNB when assumptions are violated. This demonstration highlights the expected results across various scenarios and aids in determining the most suitable methodological approach for a given problem.

In summary, our main contributions are the following:

* We propose a hybrid framework for the causal discovery of time series that combines parts of noise-based and constraint-based algorithms. Within this framework, we derive two classes of algorithms, NBCB and CBNB, which we optimize to infer the causal graph from time series.

* We study theoretically to which extent each class of algorithms is robust against assumption violation.

* We provide extensive simulation studies and real data applications to illustrate the applicability of our approach and their enhanced capabilities against assumption violation compared to original methods.

* We introduce two novel datasets about IT monitoring within the set of considered real datasets.

The remainder of the paper is organized as follows: Section 2 describes the different types of causal graphs that can be used to represent causal relations between time series and the different assumptions related to those graphs. Section 3 discusses related work and particularly details the steps that compose noise-based and constraint-based algorithms.

<footnote n="2">²In Assaad et al. [2021], this method was called NBCB but in this work, we denote it as NBCB<sup>acyclic</sup> to explicitly point out that it assumes that the summary causal graph is acyclic. For more details about summary causal graphs see Section 2.</footnote>
$$
\begin{align}
X_t &:= a_x X_{t-1} + a_{yx} Y_t + a_{zx} Z_t + \xi_t^x \nonumber \\
Z_t &:= a_z Z_{t-1} + a_{xz} X_{t-1} + a_{yz} Y_t + a_{wz} W_{t-1} + \xi_t^z \nonumber \\
Y_t &:= a_y Y_{t-1} + a_{xy} X_{t-1} + a_{wy} W_{t-2} + \xi_t^y \tag{1} \\
W_t &:= a_w W_{t-1} + a_{yw} Y_t + a_{zw} Z_t + \xi_t^w \nonumber \\
U_t &:= a_u U_{t-1} + a_{wu} W_t + \xi_t^u. \nonumber
\end{align}
$$

<figure><img src="image_8.png" /><figcaption>Figure 1: Running example. Left: Dynamic structural causal model (dynamic SCM). Right: Associated full-time causal graph $G^f$ (FTCG).</figcaption></figure>

Section 4 introduces our main contribution, the hybrid framework, which consists of two classes NBCB and CBNB, each of which is detailed in dedicated subsections. In Section 5, NBCB and CBNB are compared to different causal discovery algorithms on simulated, realistic, and real datasets. Finally, Sections 6 and 7 discuss and conclude the paper.

## 2 Background

In this section, we first introduce some terminology, tools, and assumptions which are standard for the major part. We use upper case letters to denote observed random variables, lower case letters to represent deterministic constants, blackboard bold for sets, and Greek letter $\xi$ to denote noise. A directed graph is denoted as $\mathcal{G}$ and parents and descendants of $X$ in $\mathcal{G}$ are respectively denoted as $Pa_{\mathcal{G}}(X)$ and $Desc_{\mathcal{G}}(X)$. We denote $\langle X, \dots, Y \rangle$ a path in a graph which starts at node $X$ and ends at node $Y$ and we denote $\langle X - Z, \dots, W - Y \rangle$ a walk in a graph which starts at node $X$ and ends at node $Y$ (nodes are used for paths and edges are used for walks). The skeleton of a directed graph, is a graph that consists only of undirected edges that represent the same adjacencies as in the directed graph [Hasan et al., 2023].

Causal relations in a dynamical system can be represented by a dynamic structural causal model, an extension of structural causal model [SCM, Pearl, 2000] to time series. In such dynamic SCM, each point in a time series is defined by the function of its parents and some unobserved noise. Without loss of generality, we will use the linear dynamic SCM represented in Equation (1) as a running example.

Such a dynamic SCM can be represented graphically using a full-time causal graph (FTCG, Assaad et al. [2022a]), as represented in Figure 1 and denoted as $G^f = (\mathbb{E}^f, \mathbb{V}^f)$. The FTCG represents an infinite graph of the dynamical system through infinite nodes $\mathbb{V}^f$, representing observational random variables, and infinite edges $\mathbb{E}^f$, representing causal relationships. In this paper, we only consider FTCGs that are directed acyclic graphs (DAGs), i.e., all edges are directed and there exists no directed path in $G^f$ that starts and ends at the same node.

### 2.1 Assumptions

We assume that there are no hidden common causes, an assumption known as causal sufficiency.

**Assumption 1 (Causal sufficiency, Spirtes et al., 2000).** Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG. There exist no hidden common causes of any two observed nodes in the $\mathbb{V}^f$, i.e., the noise terms in the underlying dynamic SCM are jointly independent.

One of the most common assumptions is the causal Markov condition, which is assumed by most of the methods, and connects the causal graphs that correspond to the given SCM with the compatible probability distribution.

**Assumption 2 (Causal Markov condition, Spirtes et al., 2000).** Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG and $P$ be a probability distribution over the nodes in $\mathbb{V}^f$ generated by the causal structure represented by $G^f$. Every $X \in \mathbb{V}^f$ is independent of $\mathbb{V}^f \setminus \{Desc_{\mathcal{G}^f}(X) \cup Pa_{\mathcal{G}^f}(X)\}$ conditional on $Pa_{\mathcal{G}^f}(X)$.

In general, inferring causal graphs from data is possible under additional assumptions on the data-generating process. For the constraint-based family of methods, the necessary assumption for the correspondence between the graph
$$
\begin{align}
X_t &:= a_x X_{t-1} + a_{yx} Y_t + a_{zx} Z_t + \xi_t^x \nonumber \\
Z_t &:= a_z Z_{t-1} + a_{xz} X_{t-1} + a_{yz} Y_t + a_{wz} W_{t-1} + \xi_t^z \nonumber \\
Y_t &:= a_y Y_{t-1} + a_{xy} X_{t-1} + a_{wy} W_{t-2} + \xi_t^y \tag{1} \\
W_t &:= a_w W_{t-1} + a_{yw} Y_t + a_{zw} Z_t + \xi_t^w \nonumber \\
U_t &:= a_u U_{t-1} + a_{wu} W_t + \xi_t^u. \nonumber
\end{align}
$$

<figure><img src="image_8.png" /><figcaption>Figure 1: Running example. Left: Dynamic structural causal model (dynamic SCM). Right: Associated full-time causal graph $G^f$ (FTCG).</figcaption></figure>

Section 4 introduces our main contribution, the hybrid framework, which consists of two classes NBCB and CBNB, each of which is detailed in dedicated subsections. In Section 5, NBCB and CBNB are compared to different causal discovery algorithms on simulated, realistic, and real datasets. Finally, Sections 6 and 7 discuss and conclude the paper.

## 2 Background

In this section, we first introduce some terminology, tools, and assumptions which are standard for the major part. We use upper case letters to denote observed random variables, lower case letters to represent deterministic constants, blackboard bold for sets, and Greek letter $\xi$ to denote noise. A directed graph is denoted as $\mathcal{G}$ and parents and descendants of $X$ in $\mathcal{G}$ are respectively denoted as $Pa_{\mathcal{G}}(X)$ and $Desc_{\mathcal{G}}(X)$. We denote $\langle X, \dots, Y \rangle$ a path in a graph which starts at node $X$ and ends at node $Y$ and we denote $\langle X - Z, \dots, W - Y \rangle$ a walk in a graph which starts at node $X$ and ends at node $Y$ (nodes are used for paths and edges are used for walks). The skeleton of a directed graph, is a graph that consists only of undirected edges that represent the same adjacencies as in the directed graph [Hasan et al., 2023].

Causal relations in a dynamical system can be represented by a dynamic structural causal model, an extension of structural causal model [SCM, Pearl, 2000] to time series. In such dynamic SCM, each point in a time series is defined by the function of its parents and some unobserved noise. Without loss of generality, we will use the linear dynamic SCM represented in Equation (1) as a running example.

Such a dynamic SCM can be represented graphically using a full-time causal graph (FTCG, Assaad et al. [2022a]), as represented in Figure 1 and denoted as $G^f = (\mathbb{E}^f, \mathbb{V}^f)$. The FTCG represents an infinite graph of the dynamical system through infinite nodes $\mathbb{V}^f$, representing observational random variables, and infinite edges $\mathbb{E}^f$, representing causal relationships. In this paper, we only consider FTCGs that are directed acyclic graphs (DAGs), i.e., all edges are directed and there exists no directed path in $G^f$ that starts and ends at the same node.

### 2.1 Assumptions

We assume that there are no hidden common causes, an assumption known as causal sufficiency.

**Assumption 1 (Causal sufficiency, Spirtes et al., 2000).** Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG. There exist no hidden common causes of any two observed nodes in the $\mathbb{V}^f$, i.e., the noise terms in the underlying dynamic SCM are jointly independent.

One of the most common assumptions is the causal Markov condition, which is assumed by most of the methods, and connects the causal graphs that correspond to the given SCM with the compatible probability distribution.

**Assumption 2 (Causal Markov condition, Spirtes et al., 2000).** Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG and $P$ be a probability distribution over the nodes in $\mathbb{V}^f$ generated by the causal structure represented by $G^f$. Every $X \in \mathbb{V}^f$ is independent of $\mathbb{V}^f \setminus \{Desc_{\mathcal{G}^f}(X) \cup Pa_{\mathcal{G}^f}(X)\}$ conditional on $Pa_{\mathcal{G}^f}(X)$.

In general, inferring causal graphs from data is possible under additional assumptions on the data-generating process. For the constraint-based family of methods, the necessary assumption for the correspondence between the graph
<figure><img src="image_9.png" /><figcaption>Figure 2: Illustration of the violation of the faithfulness and the adjacency faithfulness assumption. (a) violates the faithfulness assumption but satisfies the adjacency faithfulness assumption and (b) violates the adjacency faithfulness which implies that it violates the faithfulness assumption.</figcaption></figure>

and the distribution is the faithfulness assumption [Spirtes et al., 2000], which states that all the conditional independence relations that are true in the probability distribution are entailed by the causal Markov condition applied to $G^f$. However, we consider in this work the following weaker version of the faithfulness assumption, called adjacency faithfulness.

**Assumption 3** (Adjacency Faithfulness, Ramsey et al., 2006). Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG. If two nodes X and Y in $\mathbb{V}^f$ are adjacent in $G^f$, then they are dependent conditionally on any subset of $\mathbb{V}^f\setminus\{X, Y\}$.

We provide an example of the violation of the adjacency faithfulness and faithfulness assumptions in Figure 2. Let us consider variables $W_t$ and $U_t$ and their past $(W_{t-1}, U_{t-1}, \dots)$ in the causal graph in Figure 2a. Let us assume that the corresponding linear SCM is composed of two equations:

$$ 
\begin{aligned}
W_t &= a_w W_{t-1} + \xi_t^w \\
U_t &= a_u U_{t-1} + a_{wu} W_t + \xi_t^u,
\end{aligned}
$$

which leads to

$$ U_t = a_u^2 W_{t-2} + (a_u a_{wu} + a_w w_{wu}) W_{t-1} + a_{wu} \xi_t^w + a_u \xi_{t-1}^u + \xi_t^u. $$

If the coefficients are such that $a_u = -a_w$, then the coefficient before $W_{t-1}$ is 0 and thus $W_{t-1} \perp U_t | U_{t-2}$ (the noise terms are jointly independent and independent of $U_{t-2}$ and $W_{t-1}$). If there exists such a combination of coefficients, we say that the path $\langle W_{t-1}, W_t, U_t \rangle$ is canceled by the path $\langle W_{t-1}, U_{t-1}, U_t \rangle$ (or vice versa). The independence $W_{t-1} \perp U_t | U_{t-2}$ is not entailed by the the causal Markov condition, so the faithfulness assumption is violated, while adjacency faithfulness is not violated as $W_{t-1}$ and $U_t$ are not adjacent. Let us now consider the graph with an extra edge as designed in Figure 2b. If the path $\langle W_{t-1}, U_t \rangle$ is canceled out by the two paths $\langle W_{t-1}, W_t, U_t \rangle$ and $\langle W_{t-1}, U_{t-1}, U_t \rangle$ then $W_{t-1} \perp U_t | \{U_{t-2}, W_{t-2}\}$ and the adjacency faithfulness assumption is violated. Notice that $W_{t-1} \perp U_t | \{U_{t-2}, W_{t-2}\}$ is not entailed by the the causal Markov condition, which means that the faithfulness assumption is also violated.

Note that adjacency faithfulness is proven to be equivalent to the minimality condition [Peters et al., 2017] introduced in Spirtes et al. [2000], which states that every proper subgraph of $G^f$ does not satisfy causal Markov condition.

In addition, the following assumption is required by the noise-based family to guarantee identifiability.

**Assumption 4** (Identifiable Functional Model, Peters et al., 2011). The SCM belongs to an identifiable functional model class, as defined in Peters et al. [2011].

For example, the linear dynamic SCM in Equation (1) satisfies Assumption 4 if the noise terms are non-Gaussian.

## 2.2 Causal graphs for time series

Using time series data has a particular advantage for discovering causal relationships between temporal variables since we can employ temporal priority, which states that the causal relationship between two variables is oriented such that the effect cannot happen before its cause. However, despite this advantage, working with full-time causal graphs (which was introduced before) is impractical due to their infinite dimension, which has led to the adoption of simpler causal graphs, assuming that causal relations between time series hold throughout time. This is formalized in the following assumption, in which $\mathbb{V}^f = (\mathbb{V}_{-\infty}, \dots, \mathbb{V}_t, \dots, \mathbb{V}_{\infty})$ where $\mathbb{V}$ denotes a vector of d time series.

**Assumption 5** (Consistency Throughout Time). Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be a full-time causal graph. There exists $\gamma$ in $\mathbb{N} \setminus \{0\}$ such that the causal structure of the graph consisting of the nodes $\{\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t\}$ is the same for every t.

We call the minimum value of $\gamma$ for which consistency throughout time holds the maximal temporal lag of the graph. For example, the maximal temporal lag of the graph of Figure 1 is $\gamma = 2$.
<figure><img src="image_9.png" /><figcaption>Figure 2: Illustration of the violation of the faithfulness and the adjacency faithfulness assumption. (a) violates the faithfulness assumption but satisfies the adjacency faithfulness assumption and (b) violates the adjacency faithfulness which implies that it violates the faithfulness assumption.</figcaption></figure>

and the distribution is the faithfulness assumption [Spirtes et al., 2000], which states that all the conditional independence relations that are true in the probability distribution are entailed by the causal Markov condition applied to $G^f$. However, we consider in this work the following weaker version of the faithfulness assumption, called adjacency faithfulness.

**Assumption 3** (Adjacency Faithfulness, Ramsey et al., 2006). Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be an FTCG. If two nodes X and Y in $\mathbb{V}^f$ are adjacent in $G^f$, then they are dependent conditionally on any subset of $\mathbb{V}^f\{X, Y\}$.

We provide an example of the violation of the adjacency faithfulness and faithfulness assumptions in Figure 2. Let us consider variables $W_t$ and $U_t$ and their past $(W_{t-1}, U_{t-1},...)$ in the causal graph in Figure 2a. Let us assume that the corresponding linear SCM is composed of two equations:

$$ 
\begin{aligned}
W_t &= a_w W_{t-1} + \xi_t^w \\
U_t &= a_u U_{t-1} + a_{wu} W_t + \xi_t^u,
\end{aligned}
$$

which leads to

$$ U_t = a_u^2 W_{t-2} + (a_u a_{wu} + a_w w_{wu}) W_{t-1} + a_{wu} \xi_t^w + a_u \xi_{t-1}^u + \xi_t^u. $$

If the coefficients are such that $a_u = -a_w$, then the coefficient before $W_{t-1}$ is 0 and thus $W_{t-1} \perp U_t | U_{t-2}$ (the noise terms are jointly independent and independent of $U_{t-2}$ and $W_{t-1}$). If there exists such a combination of coefficients, we say that the path $\langle W_{t-1}, W_t, U_t \rangle$ is canceled by the path $\langle W_{t-1}, U_{t-1}, U_t \rangle$ (or vice versa). The independence $W_{t-1} \perp U_t | U_{t-2}$ is not entailed by the the causal Markov condition, so the faithfulness assumption is violated, while adjacency faithfulness is not violated as $W_{t-1}$ and $U_t$ are not adjacent. Let us now consider the graph with an extra edge as designed in Figure 2b. If the path $\langle W_{t-1}, U_t \rangle$ is canceled out by the two paths $\langle W_{t-1}, W_t, U_t \rangle$ and $\langle W_{t-1}, U_{t-1}, U_t \rangle$ then $W_{t-1} \perp U_t | \{U_{t-2}, W_{t-2}\}$ and the adjacency faithfulness assumption is violated. Notice that $W_{t-1} \perp U_t | \{U_{t-2}, W_{t-2}\}$ is not entailed by the the causal Markov condition, which means that the faithfulness assumption is also violated.

Note that adjacency faithfulness is proven to be equivalent to the minimality condition [Peters et al., 2017] introduced in Spirtes et al. [2000], which states that every proper subgraph of $G^f$ does not satisfy causal Markov condition.

In addition, the following assumption is required by the noise-based family to guarantee identifiability.

**Assumption 4** (Identifiable Functional Model, Peters et al., 2011). The SCM belongs to an identifiable functional model class, as defined in Peters et al. [2011].

For example, the linear dynamic SCM in Equation (1) satisfies Assumption 4 if the noise terms are non-Gaussian.

## 2.2 Causal graphs for time series

Using time series data has a particular advantage for discovering causal relationships between temporal variables since we can employ temporal priority, which states that the causal relationship between two variables is oriented such that the effect cannot happen before its cause. However, despite this advantage, working with full-time causal graphs (which was introduced before) is impractical due to their infinite dimension, which has led to the adoption of simpler causal graphs, assuming that causal relations between time series hold throughout time. This is formalized in the following assumption, in which $\mathbb{V}^f = (\mathbb{V}_{-\infty}, \dots, \mathbb{V}_t, \dots, \mathbb{V}_{\infty})$ where $\mathbb{V}$ denotes a vector of d time series.

**Assumption 5** (Consistency Throughout Time). Let $G^f = (\mathbb{E}^f, \mathbb{V}^f)$ be a full-time causal graph. There exists $\gamma$ in $\mathbb{N} \setminus \{0\}$ such that the causal structure of the graph consisting of the nodes $\{\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t\}$ is the same for every t.

We call the minimum value of $\gamma$ for which consistency throughout time holds the maximal temporal lag of the graph. For example, the maximal temporal lag of the graph of Figure 1 is $\gamma = 2$.
<figure><img src="image_10.png" /><figcaption>Figure 3: Different causal graphs to represent the dynamic SCM in Equation (1): (a) window causal graph (WCG) with a maximal temporal lag equal to 2, (b) extended summary causal graph (ECG), and (c) summary causal graph (SCG).</figcaption></figure>

Under Assumption 5, the full-time causal graph is equivalent to the so-called window causal graph [Hyvärinen et al., 2008, Runge, 2020, Figure 3a].

**Definition 1** (Window Causal Graph, WCG). Let $\mathcal{G}^f = (\mathbb{E}^f, \mathbb{V}^f)$ be a full-time causal graph satisfying Assumption 5 with $\gamma$ the maximal temporal lag in $\mathcal{G}^f$. A window causal graph (WCG) $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ is the subgraph of $\mathcal{G}^f$ consisting of the nodes $\mathbb{V}^w = (\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t)$ and $\mathbb{E}^w$ contains all related edges.

Unfortunately, causal discovery methods still suffer in practice of the strong assumptions they rely on which are not always satisfied. Thus, in many applications, experts have to validate those graphs before using them. However, validating WCGs is challenging: if experts can usually identify causes and related effects, they do not know, in general, the exact temporal lags between them. Moreover, in some applications, experts are not interested in understanding exact causal relations with the temporal lags between them but rather opt for an abstract representation of the causal relationships between variables. An extended summary causal graph [Assaad et al., 2022c, Figure 3b] is one type of abstraction that can differentiate between instantaneous relations and lagged relations without giving precise information about the lag.

**Definition 2** (Extended Summary Causal Graph, ECG). Let $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ be a WCG with maximal temporal lag $\gamma$ and nodes $(\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t)$. The underlying extended summary causal graph (ECG) $\mathcal{G}^e = (\mathbb{E}^e, \mathbb{V}^e)$ consists of the nodes $\mathbb{V}^e = (\mathbb{V}_{t-}, \mathbb{V}_t)$ and the set of directed edges $\mathbb{E}^e$ is defined as follows: for all $X_t, Y_t \in \mathbb{V}_t$, $X_t \neq Y_t$, there exists a directed edge from $X_t$ to $Y_t$ (denoted as $X_t \to Y_t$) if and only if the same directed edge exists in $\mathcal{G}^w$; for all $X_{t-} \in \mathbb{V}_{t-}$, $Y_t \in \mathbb{V}_t$, there exists a directed edge from $X_{t-}$ to $Y_t$ (denoted as $X_{t-} \to Y_t$) if and only if there exists at least one temporal lag $\ell > 0$ such that there exists a directed edge between $X_{t-\ell}$ and $Y_t$ in $\mathcal{G}^w$.

The ECG presents several interesting characteristics: similarly to the WCG, it inherits the acyclicity of the FTCG; it is possible to infer it from data without passing a WCG without any additional assumption (it only needs the assumptions that are needed to infer a WCG); the ECG and the WCG are equivalent when the maximal temporal lag $\gamma$ is 1.

Another more abstract graphical representation is the summary causal graph [Arnold et al., 2007, Peters et al., 2013, Assaad et al., 2022a, Figure 3c], which provides an overview of the causal relationships between time series without any information about the temporal lag (it does not differentiate between instantaneous relations and lagged relations).

**Definition 3** (Summary Causal Graph, SCG). Let $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ be a WCG with maximal temporal lag $\gamma$. The underlying summary causal graph (SCG) of $\mathcal{G}^w$ is $\mathcal{G}^s = (\mathbb{E}^s, \mathbb{V}^s)$ where $\mathbb{V}^s$ contains one node for each time series and the set of directed edges $\mathbb{E}^s$ are defined as follows: for all $X, Y \in \mathbb{V}^s$, $X \neq Y$, there exists a directed edge from $X$ to $Y$ (denoted as $X \to Y$) if and only if there exists at least one temporal lag $\ell \ge 0$ such that there exists a directed edge between $X_{t-\ell}$ and $Y_t$ in $\mathcal{G}^w$.

Note that, unlike the FTCG, the WCG and the ECG, the SCG may contain cycles and in particular self-loops, and two edges in opposite directions, as illustrated in the running example in Figure 3. Unlike the ECG, it is not possible to infer the SCG directly from data without additional assumptions. For example, NBCBacyclic [Assaad et al., 2021] can discover the summary causal graph assuming it is acyclic, which limits the range of potential WCGs and ECGs.

In this paper, we focus on methods that infer the WCG or the ECG from observational data.
<figure><img src="image_10.png" /><figcaption>Figure 3: Different causal graphs to represent the dynamic SCM in Equation (1): (a) window causal graph (WCG) with a maximal temporal lag equal to 2, (b) extended summary causal graph (ECG), and (c) summary causal graph (SCG).</figcaption></figure>

Under Assumption 5, the full-time causal graph is equivalent to the so-called window causal graph [Hyvärinen et al., 2008, Runge, 2020, Figure 3a].

**Definition 1** (Window Causal Graph, WCG). Let $\mathcal{G}^f = (\mathbb{E}^f, \mathbb{V}^f)$ be a full-time causal graph satisfying Assumption 5 with $\gamma$ the maximal temporal lag in $\mathcal{G}^f$. A window causal graph (WCG) $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ is the subgraph of $\mathcal{G}^f$ consisting of the nodes $\mathbb{V}^w = (\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t)$ and $\mathbb{E}^w$ contains all related edges.

Unfortunately, causal discovery methods still suffer in practice of the strong assumptions they rely on which are not always satisfied. Thus, in many applications, experts have to validate those graphs before using them. However, validating WCGs is challenging: if experts can usually identify causes and related effects, they do not know, in general, the exact temporal lags between them. Moreover, in some applications, experts are not interested in understanding exact causal relations with the temporal lags between them but rather opt for an abstract representation of the causal relationships between variables. An extended summary causal graph [Assaad et al., 2022c, Figure 3b] is one type of abstraction that can differentiate between instantaneous relations and lagged relations without giving precise information about the lag.

**Definition 2** (Extended Summary Causal Graph, ECG). Let $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ be a WCG with maximal temporal lag $\gamma$ and nodes $(\mathbb{V}_{t-\gamma}, \dots, \mathbb{V}_t)$. The underlying extended summary causal graph (ECG) $\mathcal{G}^e = (\mathbb{E}^e, \mathbb{V}^e)$ consists of the nodes $\mathbb{V}^e = (\mathbb{V}_{t-}, \mathbb{V}_t)$ and the set of directed edges $\mathbb{E}^e$ is defined as follows: for all $X_t, Y_t \in \mathbb{V}_t$, $X_t \neq Y_t$, there exists a directed edge from $X_t$ to $Y_t$ (denoted as $X_t \to Y_t$) if and only if the same directed edge exists in $\mathcal{G}^w$; for all $X_{t-} \in \mathbb{V}_{t-}$, $Y_t \in \mathbb{V}_t$, there exists a directed edge from $X_{t-}$ to $Y_t$ (denoted as $X_{t-} \to Y_t$) if and only if there exists at least one temporal lag $\ell > 0$ such that there exists a directed edge between $X_{t-\ell}$ and $Y_t$ in $\mathcal{G}^w$.

The ECG presents several interesting characteristics: similarly to the WCG, it inherits the acyclicity of the FTCG; it is possible to infer it from data without passing a WCG without any additional assumption (it only needs the assumptions that are needed to infer a WCG); the ECG and the WCG are equivalent when the maximal temporal lag $\gamma$ is 1.

Another more abstract graphical representation is the summary causal graph [Arnold et al., 2007, Peters et al., 2013, Assaad et al., 2022a, Figure 3c], which provides an overview of the causal relationships between time series without any information about the temporal lag (it does not differentiate between instantaneous relations and lagged relations).

**Definition 3** (Summary Causal Graph, SCG). Let $\mathcal{G}^w = (\mathbb{E}^w, \mathbb{V}^w)$ be a WCG with maximal temporal lag $\gamma$. The underlying summary causal graph (SCG) of $\mathcal{G}^w$ is $\mathcal{G}^s = (\mathbb{E}^s, \mathbb{V}^s)$ where $\mathbb{V}^s$ contains one node for each time series and the set of directed edges $\mathbb{E}^s$ are defined as follows: for all $X, Y \in \mathbb{V}^s$, $X \neq Y$, there exists a directed edge from $X$ to $Y$ (denoted as $X \to Y$) if and only if there exists at least one temporal lag $\ell \ge 0$ such that there exists a directed edge between $X_{t-\ell}$ and $Y_t$ in $\mathcal{G}^w$.

Note that, unlike the FTCG, the WCG and the ECG, the SCG may contain cycles and in particular self-loops, and two edges in opposite directions, as illustrated in the running example in Figure 3. Unlike the ECG, it is not possible to infer the SCG directly from data without additional assumptions. For example, NBCBacyclic [Assaad et al., 2021] can discover the summary causal graph assuming it is acyclic, which limits the range of potential WCGs and ECGs.

In this paper, we focus on methods that infer the WCG or the ECG from observational data.
We use throughout for clarity the notation $G^* = (V^*, E^*)$ to refer to either $G^e = (E^e, V^e)$ or $G^w = (E^w, V^w)$. Moreover, we use the notation $t^*$ as a wildcard that represents time step $t-i$, where $i \in \{0, \dots, \gamma\}$ in the context of WCG, and that represents either $t$ or $t$ in the context of ECG.

Another important notion that we will use in this paper is the notion of a causal order [Peters et al., 2017] which we define, in the following, only for instantaneous nodes.

**Definition 4** (Causal Order of Instantaneous Nodes). Given a causal graph $G^*$, we call a bijective mapping $\pi: V_t \mapsto \{1, \dots, d\}$, a causal order of instantaneous nodes if it satisfies

$$\pi(X_t) < \pi(Y_t) \quad \text{if } Y_t \in \text{Desc}_{G^*}(X_t), \forall X_t, Y_t \in V_t.$$

Because of the acyclicity assumption of the WCG (resp. ECG), there is always a causal order between instantaneous nodes (and a fortiori between all nodes) but it is not necessarily unique [Peters et al., 2017]. For example, in Figure 3a, there exists a causal order $\pi_1$ such that $\pi_1(X_t) = 3$ and $\pi_1(W_t) = 4$ and there exists another causal order $\pi_2$ such that $\pi_2(W_t) = 3$ and $\pi_2(X_t) = 4$, while $\pi_1(Y_t) = \pi_2(Y_t) = 1$, $\pi_1(Z_t) = \pi_2(Z_t) = 2$ and $\pi_1(U_t) = \pi_2(U_t) = 5$. Note that the number of possible causal orders is not limited to two, for example, we can obtain a new causal order from $\pi_2$, by inverting the positions of $U_t$ and $X_t$.

# 3 Related works

In this section, we start by giving a general literature review and then we give additional details on the steps that compose noise-based and constraint-based algorithms.

## 3.1 Literature review

In this section, we describe several methods from the main causal discovery families related to our work. More details can be found in thorough reviews such as Assaad et al. [2022a], Hasan et al. [2023], Gong et al. [2023]. All those methods rely on consistency throughout time (Assumption 5) to reduce the infinite dimension of the full-time causal graph. Moreover, to consider meaningful graphs, causal sufficiency and causal Markov condition (Assumptions 1 and 2) are assumed.

*Granger Causality*, introduced by Granger in 1969 and improved in subsequent works [Granger, 2004, Arnold et al., 2007], stands as one of the earliest methods designed for identifying causal relationships among time series. This approach primarily considers linear relationships and temporal priorities, operating under the assumption that the past of a cause is both necessary and sufficient for optimally forecasting its effect. It is important to note that Granger causality, by relying on temporal priority, is limited in its ability to infer instantaneous causal relations. Furthermore, Granger causality methods typically utilize the past of one time series, up to a defined maximal temporal lag, to predict the present value of another time series without explicitly distinguishing between the importance of different lags. Consequently, these methods can be employed in constructing an ECG if we assume that there are no instantaneous relations.

*Constraint-based approaches* [Spirtes et al., 2000] are certainly the most popular approaches for discovering causal graphs. These methods are based on conditional independence tests, do not depend on any specific distribution form, and require the faithfulness assumption. In the case of causal sufficiency (Assumption 1), they are usually based on the PC-algorithm [Spirtes et al., 2000], initially introduced for non-temporal data. In theory, a constraint-based algorithm can only infer a representative [known as a CPDAG, Chickering, 2002a] of the Markov equivalence class [Verma and Pearl, 1990]. Fast Causal Inference (FCI) algorithm can be used when Assumption 1 is not satisfied and infers a Partially Ancestral Graph (PAG), which is a representative of a class of equivalent Maximal ancestral graphs (MAG) and allows representing the existence of hidden confounders in the causal graph. For time series, several algorithms infer WCGs, such as PCMCI [Runge et al., 2019] and PCMCI$^+$ [Runge, 2020], extensions of the PC-algorithm, and tsFCI [Entner and Hoyer, 2010] or SVAR-FCI [Malinsky and Spirtes, 2018], extensions of the FCI algorithm. PCGCE proposed by Assaad et al. [2022c], is based on the adaptation of the PC-algorithm and infers directly an ECG.

*Score-based approaches* [Chickering, 2002a] search over the space of possible graphs, trying to maximize a score that reflects how well the graph fits the data. Greedy Equivalence Search [GES, Chickering, 2002a] is one of the first score-based methods that, under the faithfulness assumption, finds the CPDAG similarly to the PC-algorithm. There have been several recent modifications of GES, such as the more efficient Fast Greedy Search [FGS, Ramsey, 2015] and Selective Greedy Equivalence Search [Chickering and Meek, 2015].

*Continuous optimization-based approaches* Unlike score-based methods, this category of algorithms avoids the combinatorial greedy search over DAGs by employing gradient-based optimization. Notears [Zheng et al., 2018] is considered as the first approach to reformulate the greedy search as a continuous optimization problem within this category.
We use throughout for clarity the notation $G^* = (\mathbb{V}^*, \mathbb{E}^*)$ to refer to either $G^e = (\mathbb{E}^e, \mathbb{V}^e)$ or $G^w = (\mathbb{E}^w, \mathbb{V}^w)$. Moreover, we use the notation $t^*$ as a wildcard that represents time step $t-i$, where $i \in \{0, \dots, \gamma\}$ in the context of WCG, and that represents either $t$ or $t$ in the context of ECG.

Another important notion that we will use in this paper is the notion of a causal order [Peters et al., 2017] which we define, in the following, only for instantaneous nodes.

**Definition 4 (Causal Order of Instantaneous Nodes).** Given a causal graph $G^*$, we call a bijective mapping $\pi: \mathbb{V}_t \mapsto \{1, \dots, d\}$, a causal order of instantaneous nodes if it satisfies

$$ \pi(X_t) < \pi(Y_t) \quad \text{if } Y_t \in \text{Desc}_{G^*}(X_t), \forall X_t, Y_t \in \mathbb{V}_t. $$

Because of the acyclicity assumption of the WCG (resp. ECG), there is always a causal order between instantaneous nodes (and a fortiori between all nodes) but it is not necessarily unique [Peters et al., 2017]. For example, in Figure 3a, there exists a causal order $\pi_1$ such that $\pi_1(X_t) = 3$ and $\pi_1(W_t) = 4$ and there exists another causal order $\pi_2$ such that $\pi_2(W_t) = 3$ and $\pi_2(X_t) = 4$, while $\pi_1(Y_t) = \pi_2(Y_t) = 1$, $\pi_1(Z_t) = \pi_2(Z_t) = 2$ and $\pi_1(U_t) = \pi_2(U_t) = 5$. Note that the number of possible causal orders is not limited to two, for example, we can obtain a new causal order from $\pi_2$, by inverting the positions of $U_t$ and $X_t$.

# 3 Related works

In this section, we start by giving a general literature review and then we give additional details on the steps that compose noise-based and constraint-based algorithms.

## 3.1 Literature review

In this section, we describe several methods from the main causal discovery families related to our work. More details can be found in thorough reviews such as Assaad et al. [2022a], Hasan et al. [2023], Gong et al. [2023]. All those methods rely on consistency throughout time (Assumption 5) to reduce the infinite dimension of the full-time causal graph. Moreover, to consider meaningful graphs, causal sufficiency and causal Markov condition (Assumptions 1 and 2) are assumed.

**Granger Causality**, introduced by Granger in 1969 and improved in subsequent works [Granger, 2004, Arnold et al., 2007], stands as one of the earliest methods designed for identifying causal relationships among time series. This approach primarily considers linear relationships and temporal priorities, operating under the assumption that the past of a cause is both necessary and sufficient for optimally forecasting its effect. It is important to note that Granger causality, by relying on temporal priority, is limited in its ability to infer instantaneous causal relations. Furthermore, Granger causality methods typically utilize the past of one time series, up to a defined maximal temporal lag, to predict the present value of another time series without explicitly distinguishing between the importance of different lags. Consequently, these methods can be employed in constructing an ECG if we assume that there are no instantaneous relations.

**Constraint-based approaches** [Spirtes et al., 2000] are certainly the most popular approaches for discovering causal graphs. These methods are based on conditional independence tests, do not depend on any specific distribution form, and require the faithfulness assumption. In the case of causal sufficiency (Assumption 1), they are usually based on the PC-algorithm [Spirtes et al., 2000], initially introduced for non-temporal data. In theory, a constraint-based algorithm can only infer a representative [known as a CPDAG, Chickering, 2002a] of the Markov equivalence class [Verma and Pearl, 1990]. Fast Causal Inference (FCI) algorithm can be used when Assumption 1 is not satisfied and infers a Partially Ancestral Graph (PAG), which is a representative of a class of equivalent Maximal ancestral graphs (MAG) and allows representing the existence of hidden confounders in the causal graph. For time series, several algorithms infer WCGs, such as PCMCI [Runge et al., 2019] and PCMCI$^+$ [Runge, 2020], extensions of the PC-algorithm, and tsFCI [Entner and Hoyer, 2010] or SVAR-FCI [Malinsky and Spirtes, 2018], extensions of the FCI algorithm. PCGCE proposed by Assaad et al. [2022c], is based on the adaptation of the PC-algorithm and infers directly an ECG.

**Score-based approaches** [Chickering, 2002a] search over the space of possible graphs, trying to maximize a score that reflects how well the graph fits the data. Greedy Equivalence Search [GES, Chickering, 2002a] is one of the first score-based methods that, under the faithfulness assumption, finds the CPDAG similarly to the PC-algorithm. There have been several recent modifications of GES, such as the more efficient Fast Greedy Search [FGS, Ramsey, 2015] and Selective Greedy Equivalence Search [Chickering and Meek, 2015].

**Continuous optimization-based approaches** Unlike score-based methods, this category of algorithms avoids the combinatorial greedy search over DAGs by employing gradient-based optimization. Notears [Zheng et al., 2018] is considered as the first approach to reformulate the greedy search as a continuous optimization problem within this category.
This method assumes linearity and considers only DAGs whose topological order follows increasing marginal variance (rescaling data can change or reverse their inferences) [Kaiser and Sipos, 2021]. Recently, a new extension of Notears called Dynotears [Pamfil et al., 2020] was presented to infer a WCG from time series.

*Noise-based approaches* discover causal relations using footprints produced by the causal asymmetry in the data, namely the noise. They assume an identifiable functional model, as introduced in Assumption 4. For time series, one of the most popular algorithms in this family is VarLiNGAM [Hyvärinen et al., 2008], which can discover the WCG assuming linear autoregressive models with non-Gaussian noise. TiMINo [Peters et al., 2013] discovers the SCG assuming a nonlinear additive noise model. Note that TiMINo only assumes non-Gaussian noise when causal relations are linear. These approaches also require the minimality condition (Assumption 3).

*Hybrid frameworks* integrate methods from different families to improve the inference of the graph by mitigating the limitations of one algorithm through its combination with another algorithm. For nontemporal data, one group of methods combines constraint-based and score-based methods. Greedy Fast Causal inference [GFCI, Ogarrio et al., 2016] is a combination of the constraint-based FCI algorithm and the GES algorithm, which allows addressing FCI limited applicability on small sample size data and correcting the graph inferred by GES in the presence of latent confounders (removing causal sufficiency, stated in Assumption 1). HCM hybrid [Li et al., 2022] aims to discover the causal structure from mixed nontemporal data and combines part of the constraint-based method and greedy search adapted for mixed data. SVAR-GFCI introduced by Malinsky and Spirtes [2018] extends GFCI for temporal data. Another type of hybrid method is a combination of constraint-based and noise-based methods. Hoyer et al. [2008] presented a hybrid approach for nontemporal called PClingam that starts with a constraint-based procedure to find the pattern of the graph and then uses a noise-based procedure to orient some of the non-oriented edges (depending on the noise distribution), which allows obtaining a more informative causal graph. Another algorithm called FRITL [Chen et al., 2021] aims to infer more informative causal graphs with the presence of hidden confounders using the noise-based LiNGAM methods to refine the output of the FCI algorithm. For time series, a combination of noise-based and constraint-based methods was proposed by Assaad et al. [2021]—to infer an SCG assuming there are no cycles between different time series—which starts with inferring causal order using an additive noise model and pruning unnecessary edges using conditional independence tests. The combination of methods from different families generally relies on the assumptions required by employed methods (as for example in the case of GFCI or FRTIL), while sometimes hybridization allows for relaxation of some assumptions, e.g., faithfulness in the case of NBCB<sup>acyclic</sup> or non-Gaussian noise in the case of PClingam.

## 3.2 Details on noise-based and constraint-based methods

Let us first describe the basic steps of noise-based methods. A noise-based algorithm that infers a causal graph relies on the fact that, under an identifiable functional model (Assumption 4), a prediction model of a target node Y where the predictors are the true causes of Y should yield residuals (that represent the noise) that are independent of the causes. The procedure of such an algorithm can be divided into two main steps:

NB1. Find the causal order between instantaneous nodes $V_t$ by recursively performing regression and independence tests between the predictors and residuals (noise). But note that even if only instantaneous nodes $V_t$ are accounted for, lagged nodes are used in the regression as a means to account for confounder bias. For example, assuming linearity, we use the following regression model to compute the residuals of each $Y_t \in V_t$:

$$Y_t = \sum_{X_t \in V_t \setminus \{Y_t\}} a_{xy}X_t + \sum_{Z_{t-\ell} \in V^* \setminus V_t} a_{zy\ell}Z_{t-\ell} + \xi_t^y. \quad (2)$$

NB2. Find which set of predictors is not needed to keep the independence between other predictors and residuals. The former set of predictors is considered as not causally related to the node that is predicted and the latter set of predictors is considered as the causes of this node.

In general, step NB1 requires only Assumptions 1, 2, 4 and 5. Note that it has been shown that the causal order can be identified even if Assumption 3 is violated [Peters et al., 2014]<sup>3</sup>. However, Peters et al. [2014] only considered the identifiability of the causal order in case of violation of Assumption 3 up to the additive noise model case. However, Assumption 3 is required for step NB2. The output of step NB1 is a causal order $\hat{\pi}$. For step NB1 one can use, for example, TiMINo [Peters et al., 2013] or VarLiNGAM [Hyvärinen et al., 2008] algorithms.

<footnote n="3"><sup>3</sup>Note that Peters et al. [2014] considered the case where the minimality assumption is violated which is equivalent to the case where Assumption 3 is violated [Peters et al., 2017].</footnote>
This method assumes linearity and considers only DAGs whose topological order follows increasing marginal variance (rescaling data can change or reverse their inferences) [Kaiser and Sipos, 2021]. Recently, a new extension of Notears called Dynotears [Pamfil et al., 2020] was presented to infer a WCG from time series.

*Noise-based approaches* discover causal relations using footprints produced by the causal asymmetry in the data, namely the noise. They assume an identifiable functional model, as introduced in Assumption 4. For time series, one of the most popular algorithms in this family is VarLiNGAM [Hyvärinen et al., 2008], which can discover the WCG assuming linear autoregressive models with non-Gaussian noise. TiMINo [Peters et al., 2013] discovers the SCG assuming a nonlinear additive noise model. Note that TiMINo only assumes non-Gaussian noise when causal relations are linear. These approaches also require the minimality condition (Assumption 3).

*Hybrid frameworks* integrate methods from different families to improve the inference of the graph by mitigating the limitations of one algorithm through its combination with another algorithm. For nontemporal data, one group of methods combines constraint-based and score-based methods. Greedy Fast Causal inference [GFCI, Ogarrio et al., 2016] is a combination of the constraint-based FCI algorithm and the GES algorithm, which allows addressing FCI limited applicability on small sample size data and correcting the graph inferred by GES in the presence of latent confounders (removing causal sufficiency, stated in Assumption 1). HCM hybrid [Li et al., 2022] aims to discover the causal structure from mixed nontemporal data and combines part of the constraint-based method and greedy search adapted for mixed data. SVAR-GFCI introduced by Malinsky and Spirtes [2018] extends GFCI for temporal data. Another type of hybrid method is a combination of constraint-based and noise-based methods. Hoyer et al. [2008] presented a hybrid approach for nontemporal called PClingam that starts with a constraint-based procedure to find the pattern of the graph and then uses a noise-based procedure to orient some of the non-oriented edges (depending on the noise distribution), which allows obtaining a more informative causal graph. Another algorithm called FRITL [Chen et al., 2021] aims to infer more informative causal graphs with the presence of hidden confounders using the noise-based LiNGAM methods to refine the output of the FCI algorithm. For time series, a combination of noise-based and constraint-based methods was proposed by Assaad et al. [2021]—to infer an SCG assuming there are no cycles between different time series—which starts with inferring causal order using an additive noise model and pruning unnecessary edges using conditional independence tests. The combination of methods from different families generally relies on the assumptions required by employed methods (as for example in the case of GFCI or FRTIL), while sometimes hybridization allows for relaxation of some assumptions, e.g., faithfulness in the case of NBCB^acyclic or non-Gaussian noise in the case of PClingam.

## 3.2 Details on noise-based and constraint-based methods

Let us first describe the basic steps of noise-based methods. A noise-based algorithm that infers a causal graph relies on the fact that, under an identifiable functional model (Assumption 4), a prediction model of a target node Y where the predictors are the true causes of Y should yield residuals (that represent the noise) that are independent of the causes. The procedure of such an algorithm can be divided into two main steps:

NB1. Find the causal order between instantaneous nodes $V_t$ by recursively performing regression and independence tests between the predictors and residuals (noise). But note that even if only instantaneous nodes $V_t$ are accounted for, lagged nodes are used in the regression as a means to account for confounder bias. For example, assuming linearity, we use the following regression model to compute the residuals of each $Y_t \in V_t$:

$$Y_t = \sum_{X_t \in V_t \setminus \{Y_t\}} a_{xy} X_t + \sum_{Z_{t-\ell} \in V^* \setminus V_t} a_{zy\ell} Z_{t-\ell} + \xi_t^y. \quad (2)$$

NB2. Find which set of predictors is not needed to keep the independence between other predictors and residuals. The former set of predictors is considered as not causally related to the node that is predicted and the latter set of predictors is considered as the causes of this node.

In general, step NB1 requires only Assumptions 1, 2, 4 and 5. Note that it has been shown that the causal order can be identified even if Assumption 3 is violated [Peters et al., 2014]^3. However, Peters et al. [2014] only considered the identifiability of the causal order in case of violation of Assumption 3 up to the additive noise model case. However, Assumption 3 is required for step NB2. The output of step NB1 is a causal order $\hat{\pi}$. For step NB1 one can use, for example, TiMINo [Peters et al., 2013] or VarLiNGAM [Hyvärinen et al., 2008] algorithms.

<footnote>^3 Note that Peters et al. [2014] considered the case where the minimality assumption is violated which is equivalent to the case where Assumption 3 is violated [Peters et al., 2017].</footnote>
Turning to constraint-based methods, we focus in this paper on the methods that are based on the PC-algorithm, which requires causal sufficiency (Assumption 1). It follows that most constraint-based algorithms for WCGs or ECGs can be divided into two main steps:

CB1. Initialize a fully connected graph such that lagged relations are oriented using temporal priority (if $X_{t*}$ is adjacent to $Y_t$ and $t^* < t$ then $X_{t*} \rightarrow Y_t$) and instantaneous relations are unoriented, and then prune unnecessary edges using the following procedure:

(a) Eliminate edges between nodes $(X_{t*}, Y_t)$ if $X_{t*} \perp Y_t$.

(b) For each pair of nodes $(X_{t*}, Y_t)$ having an edge between them, and for each subset of nodes $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n=1$ such that $\forall Z_{t'} \in S, Z_{t'}$ is adjacent to $Y_t$, eliminate the edge between $X_{t*}$ and $Y_t$ if $X_{t*} \perp Y_t | S$.

(c) Iteratively repeat step (b) while increasing the size of the conditioning set $n$ by 1 until there are no more adjacent pairs $(X_{t*}, Y_t)$, such that there is a subset $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n$ that was not tested.

CB2. Orient instantaneous relations using some rules (in case of causal sufficiency, see Spirites et al., 2000, Meek, 1995).

Note that for step CB1, we only need Assumptions 1, 2, 3 and 5, while steps CB1 and CB2 together require the faithfulness assumption, which is stronger than Assumptions 3 [Ramsey et al., 2006]. The output of the CB1 step is a partially oriented graph $\hat{G}^*$ that shares the same skeleton as the true graph $G^*$ and where all lagged relations are oriented and all instantaneous relations are non-oriented.

In our work, we consider two algorithms from this family. The first one, PCMCI$^+$ [Runge, 2020], is an adaptation of the PC-algorithm to time series that can discover WCGs, and the second one, PCGCE [Assaad et al., 2022c], is an adaptation of the PC-algorithm that can discover ECGs.

# 4 Hybrids of constraint-based and noise-based algorithms

Here, we present a hybrid framework for discovering causal graphs from observational temporal data, which contains two classes of methods: noise-based-then-constraint-based (NBCB) and constraint-based-then-noise-based (CBNB). Both classes of methods operate under identical assumptions. They both necessitate the standard Assumptions 1 and 2 in addition to Assumption 5 for stationary time series. Moreover, both classes require Assumption 4 which is essential for noise-based methods. As both classes use only part of the algorithms from the constraint-based family, they do not require the faithfulness assumption but only the weaker adjacency faithfulness, stated in Assumption 3. In the rest of the paper, we assume that methods from CBNB or NBCB classes are built using the algorithms from constraint-based and noise-based families that are correct under the given assumptions.

Note that we do not need an acyclicity assumption of the true SCG, as in NBCB$^{acyclic}$, since we primarily focus on inferring a WCG or an ECG. If needed, we can then deduce the SCG from the inferred WCG or the inferred ECG.

In the following sections, we provide a detailed description of the NBCB and CBNB classes of methods.

## 4.1 NBCB class of algorithms

Each of the methods in the NBCB class has two major steps. The first part of NBCB constructs a fully connected graph such that lagged relations are oriented using temporal priority, on which the NB1 step is applied in order to find the causal order between all instantaneous nodes.

In the second part, we aim to use CB1 to prune the edges. This step can be optimized by using additional information on the causal order. We modify CB1 step to CB1', by taking into account the causal order of the nodes in the following way:

CB1'. Start with the initialization of the fully connected oriented graph as instantaneous relations can be oriented using the causal order $\pi$.

(a) Eliminate edges between nodes$^4$ $(X_{t*}, Y_t)$ if $X_{t*} \perp Y_t$.

(b) For each pair of nodes $(X_{t*}, Y_t)$ having an edge between them, and for each subset of nodes $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n=1$ such that $\forall Z_{t'} \in S, Z_{t'}$ is parent of $Y_t$, eliminate the edge between $X_{t*}$ and $Y_t$ if $X_{t*} \perp Y_t | S$.

<footnote>$^4$We recall that $t^*$ is a wildcard that represent time step $t-i$, where $i \in \{0, \dots, \gamma\}$ in the context of WCG and represent either $t-$ or $t$ in context of ECG.</footnote>
Turning to constraint-based methods, we focus in this paper on the methods that are based on the PC-algorithm, which requires causal sufficiency (Assumption 1). It follows that most constraint-based algorithms for WCGs or ECGs can be divided into two main steps:

CB1. Initialize a fully connected graph such that lagged relations are oriented using temporal priority (if $X_{t*}$ is adjacent to $Y_t$ and $t^* < t$ then $X_{t*} \rightarrow Y_t$) and instantaneous relations are unoriented, and then prune unnecessary edges using the following procedure:

(a) Eliminate edges between nodes $(X_{t*}, Y_t)$ if $X_{t*} \perp Y_t$.

(b) For each pair of nodes $(X_{t*}, Y_t)$ having an edge between them, and for each subset of nodes $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n=1$ such that $\forall Z_{t'} \in S, Z_{t'}$ is adjacent to $Y_t$, eliminate the edge between $X_{t*}$ and $Y_t$ if $X_{t*} \perp Y_t | S$.

(c) Iteratively repeat step (b) while increasing the size of the conditioning set $n$ by 1 until there are no more adjacent pairs $(X_{t*}, Y_t)$, such that there is a subset $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n$ that was not tested.

CB2. Orient instantaneous relations using some rules (in case of causal sufficiency, see Spirites et al., 2000, Meek, 1995).

Note that for step CB1, we only need Assumptions 1, 2, 3 and 5, while steps CB1 and CB2 together require the faithfulness assumption, which is stronger than Assumptions 3 [Ramsey et al., 2006]. The output of the CB1 step is a partially oriented graph $\hat{G}^*$ that shares the same skeleton as the true graph $G^*$ and where all lagged relations are oriented and all instantaneous relations are non-oriented.

In our work, we consider two algorithms from this family. The first one, PCMCI$^+$ [Runge, 2020], is an adaptation of the PC-algorithm to time series that can discover WCGs, and the second one, PCGCE [Assaad et al., 2022c], is an adaptation of the PC-algorithm that can discover ECGs.

# 4 Hybrids of constraint-based and noise-based algorithms

Here, we present a hybrid framework for discovering causal graphs from observational temporal data, which contains two classes of methods: noise-based-then-constraint-based (NBCB) and constraint-based-then-noise-based (CBNB). Both classes of methods operate under identical assumptions. They both necessitate the standard Assumptions 1 and 2 in addition to Assumption 5 for stationary time series. Moreover, both classes require Assumption 4 which is essential for noise-based methods. As both classes use only part of the algorithms from the constraint-based family, they do not require the faithfulness assumption but only the weaker adjacency faithfulness, stated in Assumption 3. In the rest of the paper, we assume that methods from CBNB or NBCB classes are built using the algorithms from constraint-based and noise-based families that are correct under the given assumptions.

Note that we do not need an acyclicity assumption of the true SCG, as in NBCB$^{acyclic}$, since we primarily focus on inferring a WCG or an ECG. If needed, we can then deduce the SCG from the inferred WCG or the inferred ECG.

In the following sections, we provide a detailed description of the NBCB and CBNB classes of methods.

## 4.1 NBCB class of algorithms

Each of the methods in the NBCB class has two major steps. The first part of NBCB constructs a fully connected graph such that lagged relations are oriented using temporal priority, on which the NB1 step is applied in order to find the causal order between all instantaneous nodes.

In the second part, we aim to use CB1 to prune the edges. This step can be optimized by using additional information on the causal order. We modify CB1 step to CB1', by taking into account the causal order of the nodes in the following way:

CB1'. Start with the initialization of the fully connected oriented graph as instantaneous relations can be oriented using the causal order $\pi$.

(a) Eliminate edges between nodes$^4$ $(X_{t*}, Y_t)$ if $X_{t*} \perp Y_t$.

(b) For each pair of nodes $(X_{t*}, Y_t)$ having an edge between them, and for each subset of nodes $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n=1$ such that $\forall Z_{t'} \in S, Z_{t'}$ is parent of $Y_t$, eliminate the edge between $X_{t*}$ and $Y_t$ if $X_{t*} \perp Y_t | S$.

<footnote>$^4$We recall that $t^*$ is a wildcard that represent time step $t-i$, where $i \in \{0, \dots, \gamma\}$ in the context of WCG and represent either $t$- or $t$ in context of ECG.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 4: Illustration of the NBCB algorithm for the running example. (a) Output of the NB1 step (b) Output of the CB1'.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 5: Illustration of items (2) and (3) in Proposition 1.</figcaption></figure>

(c) Iteratively repeat step (b) while increasing the size of the conditioning set $n$ by 1 until there are no more adjacent pairs $(X_{t*}, Y_t)$, such that there is a subset $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n$ that was not tested.

Note that in CB1' step (a) and (c) are the same as in step CB1. The NBCB class of methods can be presented as a combination of the NB1 and CB1', for which we provide a schematic illustration in Figure 4.

**Theorem 1.** Let $\mathcal{G}^f = (\mathbb{V}^f, \mathbb{E}^f)$ be an FTCG. Under Assumptions 1, 2, 3, 4, 5 and given perfect conditional independence information about all pairs of variables in $\mathbb{V}^f$, any algorithm in the NBCB class returns the correct WCG or the correct ECG compatible with $\mathcal{G}^f$.

The proof of Theorem 1 is available in Appendix A. It is important to note that after having obtained the WCG or the ECG, we can deduce the correct SCG using Definition 3.

The NBCB class is robust to the violation of Assumption 3. This means that even if Assumption 3 is violated, NBCB is still capable of providing valuable information on causal relationships, as demonstrated in the following proposition. We restrict ourselves to a class of functional identifiable models known as restrictive additive noise models class [Peters et al., 2014], which is stronger than Assumption 4.

**Proposition 1 (Violation of Assumption 3).** Under Assumptions 1, 2, 5, given that the SCM is a restrictive additive noise model [Peters et al., 2014], and given a correct causal order between instantaneous nodes, the NBCB class would give a WCG or an ECG such that, for each pair of nodes $X_{t*}$ and $Y_t$, one of the following possibilities holds true:

(1) The causal relationship between $X_{t*}$ and $Y_t$ is correctly identified.

(2) $X_{t*}$ and $Y_t$ are not adjacent in the inferred graph, but they are adjacent in the true graph.

(3) $X_{t*}$ and $Y_t$ are adjacent in the inferred graph, but they are not adjacent in the true graph.

The proof of Proposition 1 is presented in Appendix A. Proposition 1 states that if there is an oriented edge in a graph $\hat{\mathcal{G}}^*$ inferred by NBCB, then it can not have an opposite orientation in the true graph. In the following, we illustrate cases (2) and (3) of the Proposition 1. Let us consider that there exists linear dynamic SCM with the corresponding
<figure><img src="image_11.png" /><figcaption>Figure 4: Illustration of the NBCB algorithm for the running example. (a) Output of the NB1 step (b) Output of the CB1'.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 5: Illustration of items (2) and (3) in Proposition 1.</figcaption></figure>

(c) Iteratively repeat step (b) while increasing the size of the conditioning set $n$ by 1 until there are no more adjacent pairs $(X_{t*}, Y_t)$, such that there is a subset $S \subseteq \mathbb{V}^* \setminus \{X_{t*}, Y_t\}$ of size $n$ that was not tested.

Note that in CB1' step (a) and (c) are the same as in step CB1. The NBCB class of methods can be presented as a combination of the NB1 and CB1', for which we provide a schematic illustration in Figure 4.

**Theorem 1.** Let $\mathcal{G}^f = (\mathbb{V}^f, \mathbb{E}^f)$ be an FTCG. Under Assumptions 1, 2, 3, 4, 5 and given perfect conditional independence information about all pairs of variables in $\mathbb{V}^f$, any algorithm in the NBCB class returns the correct WCG or the correct ECG compatible with $\mathcal{G}^f$.

The proof of Theorem 1 is available in Appendix A. It is important to note that after having obtained the WCG or the ECG, we can deduce the correct SCG using Definition 3.

The NBCB class is robust to the violation of Assumption 3. This means that even if Assumption 3 is violated, NBCB is still capable of providing valuable information on causal relationships, as demonstrated in the following proposition. We restrict ourselves to a class of functional identifiable models known as restrictive additive noise models class [Peters et al., 2014], which is stronger than Assumption 4.

**Proposition 1 (Violation of Assumption 3).** Under Assumptions 1, 2, 5, given that the SCM is a restrictive additive noise model [Peters et al., 2014], and given a correct causal order between instantaneous nodes, the NBCB class would give a WCG or an ECG such that, for each pair of nodes $X_{t*}$ and $Y_t$, one of the following possibilities holds true:

(1) The causal relationship between $X_{t*}$ and $Y_t$ is correctly identified.

(2) $X_{t*}$ and $Y_t$ are not adjacent in the inferred graph, but they are adjacent in the true graph.

(3) $X_{t*}$ and $Y_t$ are adjacent in the inferred graph, but they are not adjacent in the true graph.

The proof of Proposition 1 is presented in Appendix A. Proposition 1 states that if there is an oriented edge in a graph $\hat{\mathcal{G}}^*$ inferred by NBCB, then it can not have an opposite orientation in the true graph. In the following, we illustrate cases (2) and (3) of the Proposition 1. Let us consider that there exists linear dynamic SCM with the corresponding
true WCG $\mathcal{G}^w$ presented in Figure 5a. As long as adjacency faithfulness is violated, we can assume that there exists a path $\langle Z_{t-2}, Z_{t-1}, Z_t, X_t \rangle$ which is canceled by the path $\langle Z_{t-2}, X_t \rangle$ (see example of faithfulness violation in Section 2.1), such that in the observed data $Z_{t-2} \perp X_t | X_{t-1}$. Thus, the edge $Z_{t-2} - X_t$ would be removed, so we obtain case (2) in Proposition 1. A similar argument works for the edge $Z_{t-2} - Y_t$. Further, these two errors would propagate in the next step for the pair of nodes $Y_t$ and $X_t$. As in the inferred graph $\hat{\mathcal{G}}^w$, edges $Z_{t-2} - X_t$ and $Z_{t-2} - Y_t$ are absent (see Figure 5b), then $X_t \perp Y_t | Pa_{\hat{\mathcal{G}}^w}(X_t) \cup Pa_{\hat{\mathcal{G}}^w}(Y_t)$, as $Z_{t-2}$ is the common confounder in the true graph, thus the edge $X_t - Y_t$ would not be removed and we obtain the case (3) in Proposition 1.

We want to highlight that Proposition 1 holds when the CB1' step of NBCB is implemented using a PC-style algorithm. Other constraint-based algorithms could be envisioned, such as a greedy approach where the iterative process of finding conditional independence while progressively expanding the conditioning set (steps (b) and (c) of CB1') is replaced by a single step of conditioning on all parents. With such a greedy algorithm, the concerns outlined in case (3) of the proposition could be alleviated.

Note that NBCB would perform poorly when the SCM is not an identifiable functional model (Assumption 4 is violated). In this case, NB1 would give an incorrect causal order, and the errors would propagate in the second step.

**Algorithm 1:** Noise-Based-then-Constraint-Based (NBCB)

**Input:** A multivariate time series, a maximal temporal lag $\gamma$ and a significance threshold $\alpha$, NB1, CB1', an independence measure $I()$, and a conditional independence test $CI()$

**Result:** $\hat{\mathcal{G}}^*$ (WCG or ECG) and $\hat{\mathcal{G}}^s$ (SCG)

Find the causal order $\hat{\pi}$ between all instantaneous nodes using NB1 which takes $\gamma$ and $I()$ as hyper-parameters;

Discover $\hat{\mathcal{G}}^*$ using CB1' which takes $\gamma$, $\alpha$, $CI()$, and $\hat{\pi}$ as hyper-parameters;

Deduce the SCG $\hat{\mathcal{G}}^s$ from $\hat{\mathcal{G}}^*$ using Definition 3.

The pseudo-code of NBCB is given in Algorithm 1, which involves the abstract steps NB1 and CB1'. Step NB1 can be directly obtained from any existing noise-based algorithm. In the experimental section of this paper, step NB1 is directly based on the VarLiNGAM algorithm, while step CB1' is either based on the modification of the PCTMI+ algorithm or on the modification of the PCGCE algorithm. More details on specific versions of NB1 and CB1' is given in Appendix B.

## 4.2 CBNB class of algorithms

As far as we know, there exists no previous work that investigated the case where the constraint-based part of a hybrid method is executed before the noise-based part in causal discovery from time series (for non-temporal data, see PClingam, Hoyer et al., 2008). In theory, there is no clear argument as to why one part should be executed before the other. So here, we present a new class of methods called CBNB, where the constraint-based part is before the noise-based part.

In the first step of CBNB, CB1 is used to infer $\hat{\mathcal{G}}^*$. Then, we are going to orient edges between instantaneous nodes $V_t$ by finding the causal order between them.

A vanilla approach would be to apply NB1 over all instantaneous nodes (while taking into account lagged common confounders), as it was done in NBCB, to get the causal order and then orient the instantaneous relations accordingly. However, we argue that this vanilla approach, despite being correct, is not optimal since it does not take into account the knowledge that has already been acquired through the construction of the output graph in step CB1. Therefore, we modify the NB1 step to NB1' by finding the causal order within different groups of instantaneous nodes separately in the following way:

NB1'. Find the causal order within instantaneous groups of nodes $I_t \subseteq V_t$ by recursively performing regression and independence tests between the predictors and residuals (noise). But note that even if only the instantaneous group of nodes $I_t$ are accounted for, the lagged parents $Pa_{\hat{\mathcal{G}}^*}(I_t) \setminus V_t$ are used in the regression as a means to account for confounder bias.

Not any groups of instantaneous nodes $I_t$ would yield a correct causal order since, in general, for any $X_t, Y_t \in I_t$, $I_t \cup Pa_{\hat{\mathcal{G}}^*}(I_t)$ does not necessarily contain any subset that remove confounding bias between $X_t$ and $Y_t$. Therefore, these groups should be selected carefully. To find these groups, we first need to define an undirected cycle walk.

**Definition 5** (Undirected Cycle Walk). Let $\hat{\mathcal{G}}^*$ be the output of the CB1 step. A sequence of edges $u = \langle e_1, \dots, e_n \rangle$, for $n \ge 2$, is an undirected cycle walk iff
true WCG $\mathcal{G}^w$ presented in Figure 5a. As long as adjacency faithfulness is violated, we can assume that there exists a path $\langle Z_{t-2}, Z_{t-1}, Z_t, X_t \rangle$ which is canceled by the path $\langle Z_{t-2}, X_t \rangle$ (see example of faithfulness violation in Section 2.1), such that in the observed data $Z_{t-2} \perp X_t | X_{t-1}$. Thus, the edge $Z_{t-2} - X_t$ would be removed, so we obtain case (2) in Proposition 1. A similar argument works for the edge $Z_{t-2} - Y_t$. Further, these two errors would propagate in the next step for the pair of nodes $Y_t$ and $X_t$. As in the inferred graph $\hat{\mathcal{G}}^w$, edges $Z_{t-2} - X_t$ and $Z_{t-2} - Y_t$ are absent (see Figure 5b), then $X_t \perp Y_t | Pa_{\hat{\mathcal{G}}^w}(X_t) \cup Pa_{\hat{\mathcal{G}}^w}(Y_t)$, as $Z_{t-2}$ is the common confounder in the true graph, thus the edge $X_t - Y_t$ would not be removed and we obtain the case (3) in Proposition 1.

We want to highlight that Proposition 1 holds when the CB1' step of NBCB is implemented using a PC-style algorithm. Other constraint-based algorithms could be envisioned, such as a greedy approach where the iterative process of finding conditional independence while progressively expanding the conditioning set (steps (b) and (c) of CB1') is replaced by a single step of conditioning on all parents. With such a greedy algorithm, the concerns outlined in case (3) of the proposition could be alleviated.

Note that NBCB would perform poorly when the SCM is not an identifiable functional model (Assumption 4 is violated). In this case, NB1 would give an incorrect causal order, and the errors would propagate in the second step.

**Algorithm 1:** Noise-Based-then-Constraint-Based (NBCB)

**Input:** A multivariate time series, a maximal temporal lag $\gamma$ and a significance threshold $\alpha$, NB1, CB1', an independence measure $I()$, and a conditional independence test $CI()$

**Result:** $\hat{\mathcal{G}}^*$ (WCG or ECG) and $\hat{\mathcal{G}}^s$ (SCG)

Find the causal order $\hat{\pi}$ between all instantaneous nodes using NB1 which takes $\gamma$ and $I()$ as hyper-parameters;

Discover $\hat{\mathcal{G}}^*$ using CB1' which takes $\gamma$, $\alpha$, $CI()$, and $\hat{\pi}$ as hyper-parameters;

Deduce the SCG $\hat{\mathcal{G}}^s$ from $\hat{\mathcal{G}}^*$ using Definition 3.

The pseudo-code of NBCB is given in Algorithm 1, which involves the abstract steps NB1 and CB1'. Step NB1 can be directly obtained from any existing noise-based algorithm. In the experimental section of this paper, step NB1 is directly based on the VarLiNGAM algorithm, while step CB1' is either based on the modification of the PCTMI+ algorithm or on the modification of the PCGCE algorithm. More details on specific versions of NB1 and CB1' is given in Appendix B.

## 4.2 CBNB class of algorithms

As far as we know, there exists no previous work that investigated the case where the constraint-based part of a hybrid method is executed before the noise-based part in causal discovery from time series (for non-temporal data, see PClingam, Hoyer et al., 2008). In theory, there is no clear argument as to why one part should be executed before the other. So here, we present a new class of methods called CBNB, where the constraint-based part is before the noise-based part.

In the first step of CBNB, CB1 is used to infer $\hat{\mathcal{G}}^*$. Then, we are going to orient edges between instantaneous nodes $\mathbb{V}_t$ by finding the causal order between them.

A vanilla approach would be to apply NB1 over all instantaneous nodes (while taking into account lagged common confounders), as it was done in NBCB, to get the causal order and then orient the instantaneous relations accordingly. However, we argue that this vanilla approach, despite being correct, is not optimal since it does not take into account the knowledge that has already been acquired through the construction of the output graph in step CB1. Therefore, we modify the NB1 step to NB1' by finding the causal order within different groups of instantaneous nodes separately in the following way:

NB1'. Find the causal order within instantaneous groups of nodes $\mathbb{I}_t \subseteq \mathbb{V}_t$ by recursively performing regression and independence tests between the predictors and residuals (noise). But note that even if only the instantaneous group of nodes $\mathbb{I}_t$ are accounted for, the lagged parents $Pa_{\hat{\mathcal{G}}^*}(\mathbb{I}_t) \setminus \mathbb{V}_t$ are used in the regression as a means to account for confounder bias.

Not any groups of instantaneous nodes $\mathbb{I}_t$ would yield a correct causal order since, in general, for any $X_t, Y_t \in \mathbb{I}_t$, $\mathbb{I}_t \cup Pa_{\hat{\mathcal{G}}^*}(\mathbb{I}_t)$ does not necessarily contain any subset that remove confounding bias between $X_t$ and $Y_t$. Therefore, these groups should be selected carefully. To find these groups, we first need to define an undirected cycle walk.

**Definition 5** (Undirected Cycle Walk). Let $\hat{\mathcal{G}}^*$ be the output of the CB1 step. A sequence of edges $\mathbf{u} = \langle e_1, \dots, e_n \rangle$, for $n \ge 2$, is an undirected cycle walk iff
<figure><img src="image_13.png" /><figcaption>Figure 6: Illustration of the CBNB algorithm for the running example. (a) Output of the CB1 step (b) Detection of the undirected cycle group, where red and blue colours denote two different cycle groups of instantaneous nodes (c) Output of the NB1' on each of the cycle groups.</figcaption></figure>

• $\forall e_i \in \mathbf{u}, e_i$ is an undirected edge in $\hat{\mathcal{G}}^*$;

• $\forall e_i, e_j \in \mathbf{u}$, such that $i < j$, no node in $e_i$ coincides with any node in $e_j$, except if $j = i + 1$ where the right-hand side node in $e_i$ always coincide with the left-hand side node in $e_{i+1}$, or if $i = 1, j = n$, where the left-hand side node in $e_1$ and the right-hand side node in $e_n$ always coincide.

Note that in our definition an undirected edge can form an undirected cycle walk. We provide a schematic illustration in Figure 6 for the running example. In Figure 6b, $\langle X_t - Y_t, Y_t - Z_t, Z_t - X_t \rangle$ and $\langle Z_t - Y_t, Y_t - W_t, W_t - Z_t \rangle$ are two undirected cycle walks. Notice that $Z_t - Y_t$ is in both undirected cycle walks and would be considered twice, which is not only a computational problem but might also induce bias in practice. Thus, defining a undirected cycle walk alone is not sufficient. So, we combine all undirected cycle walks that share at least one edge in a undirected cycle group<fnref n="5" /> which is defined as follows:

**Definition 6 (Undirected Cycle Group).** Let $\hat{\mathcal{G}}^*$ be the output of the CB1 step. $\mathcal{C}$ is an undirected cycle group of $\hat{\mathcal{G}}^*$ iff $\mathcal{C}$ is a set of undirected cycle walks and $\forall \mathbf{u}_1, \mathbf{u}_2 \in \mathcal{C}, \mathbf{u}_1 \cap \mathbf{u}_2 \neq \emptyset$.

We assume that two undirected cycle walks $\mathbf{u}_1$ and $\mathbf{u}_2$ do not intersect (i.e., $\mathbf{u}_1 \cap \mathbf{u}_2 = \emptyset$) if there is no common edge between these two walks. For simplicity, we consider two undirected cycle walks that are reversed versions of one another as equivalent and so only one of these two undirected cycle walks will be accounted for. In Figure 6b, for the WCG, there are two *undirected cycle groups*, $\mathcal{C}_1 = \{\langle X_t - Y_t, Y_t - Z_t, Z_t - X_t \rangle, \langle Z_t - Y_t, Y_t - W_t, W_t - Z_t \rangle, \dots \}$ and $\mathcal{C}_2 = \{\langle U_t - W_t, W_t - U_t \rangle\}$. The same example applies to the corresponding ECG. Given an *undirected cycle group* $\mathcal{C}$, we are interested in the set of its nodes, meaning the set of nodes that belong to undirected cycle walks described by the undirected cycle group. Similarly, we say that an edge $X_t - Y_{t'}$ belongs to $\mathcal{C}$ iff $\exists \mathbf{u} \in \mathcal{C}$ such that $\langle X_t, Y_{t'} \rangle \in \mathbf{u}$.

Having brought out the concept of *undirected cycle group*, we can now describe the second (noise-based) part of CBNB. Given the output of the CB1 step, CBNB searches for all undirected cycle groups. Note that each edge is a member of exactly one undirected cycle group. Then for each undirected cycle group $\mathcal{C}$, CBNB uses NB1' to find the causal order $\pi$ between all the nodes belonging to undirected cycle group $\mathcal{C}$. Finally, using this causal order, CBNB orient all edges belonging to $\mathcal{C}$.

Note that the CBNB step where NB1' is applied for different *undirected cycle groups* can be parallelized, which can significantly improve the computational time for high dimensional problems.

**Theorem 2.** Let $\mathcal{G}^f = (\mathbb{V}^f, \mathbb{E}^f)$ be an FTCG. Under Assumptions 1, 2, 3, 4, 5 and given perfect conditional independence information about all pairs of variables in $\mathbb{V}^f$, any algorithm in the CBNB class returns the correct WCG or the correct ECG compatible with $\mathcal{G}^f$.

The proof of Theorem 2 is available in Appendix A. It is also important to note that as in the case of NBCB, after having obtained the WCG or the ECG, we can deduce the correct SCG using Definition 3.

<footnote n="5">⁵Note that our definition of cycle group is different than the one introduced in Spirtes [1995].</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 6: Illustration of the CBNB algorithm for the running example. (a) Output of the CB1 step (b) Detection of the undirected cycle group, where red and blue colours denote two different cycle groups of instantaneous nodes (c) Output of the NB1' on each of the cycle groups.</figcaption></figure>

• $\forall e_i \in \mathbf{u}, e_i$ is an undirected edge in $\hat{\mathcal{G}}^*$;

• $\forall e_i, e_j \in \mathbf{u}$, such that $i < j$, no node in $e_i$ coincides with any node in $e_j$, except if $j = i + 1$ where the right-hand side node in $e_i$ always coincide with the left-hand side node in $e_{i+1}$, or if $i = 1, j = n$, where the left-hand side node in $e_1$ and the right-hand side node in $e_n$ always coincide.

Note that in our definition an undirected edge can form an undirected cycle walk. We provide a schematic illustration in Figure 6 for the running example. In Figure 6b, $\langle X_t - Y_t, Y_t - Z_t, Z_t - X_t \rangle$ and $\langle Z_t - Y_t, Y_t - W_t, W_t - Z_t \rangle$ are two undirected cycle walks. Notice that $Z_t - Y_t$ is in both undirected cycle walks and would be considered twice, which is not only a computational problem but might also induce bias in practice. Thus, defining a undirected cycle walk alone is not sufficient. So, we combine all undirected cycle walks that share at least one edge in a undirected cycle group<fnref n="5" /> which is defined as follows:

**Definition 6 (Undirected Cycle Group).** Let $\hat{\mathcal{G}}^*$ be the output of the CB1 step. $\mathcal{C}$ is an undirected cycle group of $\hat{\mathcal{G}}^*$ iff $\mathcal{C}$ is a set of undirected cycle walks and $\forall \mathbf{u}_1, \mathbf{u}_2 \in \mathcal{C}, \mathbf{u}_1 \cap \mathbf{u}_2 \neq \emptyset$.

We assume that two undirected cycle walks $\mathbf{u}_1$ and $\mathbf{u}_2$ do not intersect (i.e., $\mathbf{u}_1 \cap \mathbf{u}_2 = \emptyset$) if there is no common edge between these two walks. For simplicity, we consider two undirected cycle walks that are reversed versions of one another as equivalent and so only one of these two undirected cycle walks will be accounted for. In Figure 6b, for the WCG, there are two *undirected cycle groups*, $\mathcal{C}_1 = \{\langle X_t - Y_t, Y_t - Z_t, Z_t - X_t \rangle, \langle Z_t - Y_t, Y_t - W_t, W_t - Z_t \rangle, \dots \}$ and $\mathcal{C}_2 = \{\langle U_t - W_t, W_t - U_t \rangle\}$. The same example applies to the corresponding ECG. Given an *undirected cycle group* $\mathcal{C}$, we are interested in the set of its nodes, meaning the set of nodes that belong to undirected cycle walks described by the undirected cycle group. Similarly, we say that an edge $X_t - Y_t'$ belongs to $\mathcal{C}$ iff $\exists \mathbf{u} \in \mathcal{C}$ such that $\langle X_t, Y_t' \rangle \in \mathbf{u}$.

Having brought out the concept of *undirected cycle group*, we can now describe the second (noise-based) part of CBNB. Given the output of the CB1 step, CBNB searches for all undirected cycle groups. Note that each edge is a member of exactly one undirected cycle group. Then for each undirected cycle group $\mathcal{C}$, CBNB uses NB1' to find the causal order $\pi$ between all the nodes belonging to undirected cycle group $\mathcal{C}$. Finally, using this causal order, CBNB orient all edges belonging to $\mathcal{C}$.

Note that the CBNB step where NB1' is applied for different *undirected cycle groups* can be parallelized, which can significantly improve the computational time for high dimensional problems.

**Theorem 2.** Let $\mathcal{G}^f = (\mathbb{V}^f, \mathbb{E}^f)$ be an FTCG. Under Assumptions 1, 2, 3, 4, 5 and given perfect conditional independence information about all pairs of variables in $\mathbb{V}^f$, any algorithm in the CBNB class returns the correct WCG or the correct ECG compatible with $\mathcal{G}^f$.

The proof of Theorem 2 is available in Appendix A. It is also important to note that as in the case of NBCB, after having obtained the WCG or the ECG, we can deduce the correct SCG using Definition 3.

<footnote n="5">⁵Note that our definition of cycle group is different than the one introduced in Spirtes [1995].</footnote>
The CBNB class is robust to the violation of Assumption 4. This means that even if Assumption 4 is violated, CBNB is still capable of providing valuable information on causal relationships, as demonstrated in the following proposition.

**Proposition 2 (Violation of Assumption 4).** *Under Assumptions 1, 2, 3, 5 and given perfect conditional independence information about all pairs of variables, CBNB is guaranteed to find the correct skeleton of the WCG or the ECG.*

The proof of Proposition 2 is available in Appendix A. Note that since CBNB gives the correct skeleton of WCGs and ECGs, the correct skeleton of SCG can also be deduced. If Assumption 3 is violated, then the skeleton obtained in the CB1 step of CBNB is not reliable, thereby affecting the reliability of the CBNB algorithm’s results.

The pseudo-code of CBNB is given in Algorithm 2, which involves the abstract steps CB1 and NB1'. In the experimental section of this paper, step CB1 is directly based either on the PCMCI+ algorithm or on the PCGCE algorithm, while step NB1' is based on the modification of the VarLiNGAM algorithm. More details on specific versions of CB1 and NB1' is given in Appendix B.

**Algorithm 2: Constraint-Based-then-Noise-Based (CBNB)**

**Input:** A multivariate time series, a maximal temporal lag $\gamma$ and a significance threshold $\alpha$, CB1, NB1', an independence measure $I()$, and a conditional independence test $CI()$

**Result:** $\hat{G}^*$ (WCG or ECG) and $\hat{G}^s$ (SCG)

Initialize $\hat{G}^*$ as the output of CB1 which takes $\gamma$, $\alpha$, and $CI()$ as hyper-parameters;

for each undirected cycle group $C$ in $\hat{G}^*$ do

$I_t$: set of nodes that belong to $C$;

Find the causal order $\pi$ between nodes in $I_t$ using NB1' which takes $\gamma$, $I()$, and $\hat{G}^*$ as hyper-parameters;

Orient instantaneous edges between $I_t$ in $\hat{G}^*$ using the order $\pi$;

Deduce the SCG $\hat{G}^s$ from $\hat{G}^*$ using Definition 3

## 4.3 Complexity analysis for NBCB and CBNB classes of methods

When considering $d$ time series, finding the causal order is of complexity $d^2 f(n, d)$ where $f(n, d)$ is the complexity of the user-specific regression method with time series length $n$. Then, pruning the graph by conditional independence tests is of complexity $\frac{(dq)^2(dq-1)^{k-1}}{(k-1)!}$, where $k$ represents the maximal degree of any node, $q$ is equal to $\tau$ when considering a WCG and $q$ is equal to 2 when considering an ECG, and each operation consists in conducting significance test to a conditional independence measure. Both CBNB and NBCB have the same complexity in the worst case, $d^2 f(n, d) + \frac{(dq)^2(dq-1)^{k-1}}{(k-1)!}$ which corresponds for CBNB to the case where there is only one undirected cycle group.

# 5 Experiments

In this section, we propose first an extensive analysis<fnref n="6" /> both on simulated data, generated from basic causal structures, and on simulated but realistic benchmarks. We then perform an analysis on different real datasets.

## 5.1 Experimental setup

**Baselines and hyper-parameters.** We compare NBCB and CBNB with five state-of-the-art methods:

* the multivariate version of Granger Causality denoted GCMVL [Arnold et al., 2007];

* the constraint-based methods PCMCI$^+$ [Runge, 2020] for which we use the Python code available at https://github.com/jakobrunge/tigramite and PCGCE [Assaad et al., 2022c] for which the main Python code available at https://github.com/ckassad/PCGCE. For PCGCE, as the authors suggested, we reduce the dimensionality of $V_{t-}$ in the ECG to 1 using PCA;

* the continuous optimization-based method Dynotears [Pamfil et al., 2020], for which the Python code is available at https://github.com/quantumblacklabs/causalnex. For this method, we set the hyperparameters to their recommended values ($\lambda_W = \lambda_A = 0.05$ and $\alpha_W = \alpha_A = 0.01$);

<footnote n="6">⁶A Python code of all our methods and of every experimentation is available in https://github.com/ckassad/Hybrids_of_CB_and_NB_for_Time_Series.</footnote>
The CBNB class is robust to the violation of Assumption 4. This means that even if Assumption 4 is violated, CBNB is still capable of providing valuable information on causal relationships, as demonstrated in the following proposition.

**Proposition 2** (Violation of Assumption 4). *Under Assumptions 1, 2, 3, 5 and given perfect conditional independence information about all pairs of variables, CBNB is guaranteed to find the correct skeleton of the WCG or the ECG.*

The proof of Proposition 2 is available in Appendix A. Note that since CBNB gives the correct skeleton of WCGs and ECGs, the correct skeleton of SCG can also be deduced. If Assumption 3 is violated, then the skeleton obtained in the CB1 step of CBNB is not reliable, thereby affecting the reliability of the CBNB algorithm’s results.

The pseudo-code of CBNB is given in Algorithm 2, which involves the abstract steps CB1 and NB1'. In the experimental section of this paper, step CB1 is directly based either on the PCMCI+ algorithm or on the PCGCE algorithm, while step NB1' is based on the modification of the VarLiNGAM algorithm. More details on specific versions of CB1 and NB1' is given in Appendix B.

**Algorithm 2:** Constraint-Based-then-Noise-Based (CBNB)

**Input:** A multivariate time series, a maximal temporal lag $\gamma$ and a significance threshold $\alpha$, CB1, NB1', an independence measure $I()$, and a conditional independence test $CI()$

**Result:** $\hat{G}^*$ (WCG or ECG) and $\hat{G}^s$ (SCG)

Initialize $\hat{G}^*$ as the output of CB1 which takes $\gamma$, $\alpha$, and $CI()$ as hyper-parameters;

for each undirected cycle group $C$ in $\hat{G}^*$ do

$\mathbb{I}_t$: set of nodes that belong to $C$;

Find the causal order $\pi$ between nodes in $\mathbb{I}_t$ using NB1' which takes $\gamma$, $I()$, and $\hat{G}^*$ as hyper-parameters;

Orient instantaneous edges between $\mathbb{I}_t$ in $\hat{G}^*$ using the order $\pi$;

Deduce the SCG $\hat{G}^s$ from $\hat{G}^*$ using Definition 3

## 4.3 Complexity analysis for NBCB and CBNB classes of methods

When considering $d$ time series, finding the causal order is of complexity $d^2 f(n, d)$ where $f(n, d)$ is the complexity of the user-specific regression method with time series length $n$. Then, pruning the graph by conditional independence tests is of complexity $\frac{(dq)^2(dq-1)^{k-1}}{(k-1)!}$, where $k$ represents the maximal degree of any node, $q$ is equal to $\tau$ when considering a WCG and $q$ is equal to 2 when considering an ECG, and each operation consists in conducting significance test to a conditional independence measure. Both CBNB and NBCB have the same complexity in the worst case, $d^2 f(n, d) + \frac{(dq)^2(dq-1)^{k-1}}{(k-1)!}$ which corresponds for CBNB to the case where there is only one undirected cycle group.

# 5 Experiments

In this section, we propose first an extensive analysis<fnref n="6" /> both on simulated data, generated from basic causal structures, and on simulated but realistic benchmarks. We then perform an analysis on different real datasets.

## 5.1 Experimental setup

**Baselines and hyper-parameters.** We compare NBCB and CBNB with five state-of-the-art methods:

* the multivariate version of Granger Causality denoted GCMVL [Arnold et al., 2007];

* the constraint-based methods PCMCI$^+$ [Runge, 2020] for which we use the Python code available at https://github.com/jakobrunge/tigramite and PCGCE [Assaad et al., 2022c] for which the main Python code available at https://github.com/ckassad/PCGCE. For PCGCE, as the authors suggested, we reduce the dimensionality of $V_{t-}$ in the ECG to 1 using PCA;

* the continuous optimization-based method Dynotears [Pamfil et al., 2020], for which the Python code is available at https://github.com/quantumblacklabs/causalnex. For this method, we set the hyperparameters to their recommended values ($\lambda_W = \lambda_A = 0.05$ and $\alpha_W = \alpha_A = 0.01$);

<footnote n="6">⁶A Python code of all our methods and of every experimentation is available in https://github.com/ckassad/Hybrids_of_CB_and_NB_for_Time_Series.</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 7: Summary causal graphs corresponding to the data simulated in Section 5.2. The last two models correspond to an unfaithful distribution.</figcaption></figure>

* the noise-based method VarLiNGAM [Hyvärinen et al., 2008] for which we use the Python code available at https://github.com/cdt15/lingam where the regularization parameter in the adaptive Lasso is selected using BIC.

For all the methods, the maximal temporal lag is set to $\gamma = 5$ and the significant threshold for hypothesis testing to $\alpha = 0.05$. We test two versions of each of our classes which we denote NBCB-w, NBCB-e, CBNB-w, and CBNB-e. Methods with "-w" suffix are based on the algorithms that infer a WCG, while methods with "-e" are based on the algorithms that infer an ECG. In NBCB-w and NBCB-e, the NB1 step is based on the VarLiNGAM algorithm and the CB1' step is respectively based on the PCMCI+ algorithm and the PCGCE algorithm. In CBNB-w and CBNB-e, the CB1 step is respectively based on the PCMCI+ algorithm and the PCGCE algorithm and the NB1' step is based on the VarLiNGAM algorithm. The pseudo-code of each version of the NB1, CB1, NB1', CB1' steps are given in Appendix B. To find undirected cycle walks in CBNB-w and CBNB-e, we use an adapted version of Paton's algorithm [Paton, 1969]. For all methods that require a conditional independence test, we use a test based on partial correlation, which assumes Gaussian distributions but which has been successfully used on non-iid data [Peters et al., 2013].

For real data, we also use nonlinear versions of our hybrid algorithms as well of the PCMCI and the PCGCE algorithms. We use kernel-based conditional independence test which combines Gaussian process regression with a distance correlation test on the residuals [Runge et al., 2019], and a Gaussian process regression in the noise-based part. All nonlinear versions are denoted by the suffix "-nl".

**Evaluation.** In the different experimental settings, we compared the results concerning the F1 score of the orientations in the SCG obtained without considering self causes, as it is treated differently depending on the methods. When there are three datasets or more, we report the mean and the variance for the F1 score.

## 5.2 Simulated data with Gaussian and non-Gaussian noise

The simulated datasets correspond to four causally sufficient SCGs presented in Figure 7, extracted from WCGs $G^w$, among which three are acyclic, two correspond to an unfaithful distribution, and one is cyclic. The generating process of all datasets is the following: for all $Y$, for all $t > 0$,

$$Y_t = a_y Y_{t-1} + \sum_{X_{t-\ell} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-\ell} + 0.1 \xi_t^y,$$

where $a_y, a_{xy} \in U([-1, -0.1] \cup [0.1, 1])$ and for each parent $X$ we randomly choose if $X$ causes $Y$ instantaneously or with a lag of 1, i.e., $\ell \in \{0, 1\}$. Regarding the noise, we consider two different settings, in the first the noise is drawn from a uniform distribution, i.e., $\xi^y \sim U([-1, 1])$, and in the second the noise is drawn from a Gaussian distribution, i.e., $\xi^y \sim N(0, 1)$. For each setting and each structure in Figure 7 we generate 100 datasets of 1000 timestamps.

For the unfaithful diamond structure in Figure 7, following Zhalama et al. [2016], we set $a_x, a_y, a_z, a_w$ to zero, and $a_{zw} = -a_{xy}a_{yw}/a_{xz}$. All relationships are considered instantaneous. Consequently, in the distribution consistent with this model, we observe $X \perp W$, which violates the faithfulness assumption as this independence is not entailed by the
<figure><img src="image_14.png" /><figcaption>Figure 7: Summary causal graphs corresponding to the data simulated in Section 5.2. The last two models correspond to an unfaithful distribution.</figcaption></figure>

* the noise-based method VarLiNGAM [Hyvärinen et al., 2008] for which we use the Python code available at https://github.com/cdt15/lingam where the regularization parameter in the adaptive Lasso is selected using BIC.

For all the methods, the maximal temporal lag is set to $\gamma = 5$ and the significant threshold for hypothesis testing to $\alpha = 0.05$. We test two versions of each of our classes which we denote NBCB-w, NBCB-e, CBNB-w, and CBNB-e. Methods with "-w" suffix are based on the algorithms that infer a WCG, while methods with "-e" are based on the algorithms that infer an ECG. In NBCB-w and NBCB-e, the NB1 step is based on the VarLiNGAM algorithm and the CB1' step is respectively based on the PCMCI+ algorithm and the PCGCE algorithm. In CBNB-w and CBNB-e, the CB1 step is respectively based on the PCMCI+ algorithm and the PCGCE algorithm and the NB1' step is based on the VarLiNGAM algorithm. The pseudo-code of each version of the NB1, CB1, NB1', CB1' steps are given in Appendix B. To find undirected cycle walks in CBNB-w and CBNB-e, we use an adapted version of Paton's algorithm [Paton, 1969]. For all methods that require a conditional independence test, we use a test based on partial correlation, which assumes Gaussian distributions but which has been successfully used on non-iid data [Peters et al., 2013].

For real data, we also use nonlinear versions of our hybrid algorithms as well of the PCMCI and the PCGCE algorithms. We use kernel-based conditional independence test which combines Gaussian process regression with a distance correlation test on the residuals [Runge et al., 2019], and a Gaussian process regression in the noise-based part. All nonlinear versions are denoted by the suffix "-nl".

**Evaluation.** In the different experimental settings, we compared the results concerning the F1 score of the orientations in the SCG obtained without considering self causes, as it is treated differently depending on the methods. When there are three datasets or more, we report the mean and the variance for the F1 score.

## 5.2 Simulated data with Gaussian and non-Gaussian noise

The simulated datasets correspond to four causally sufficient SCGs presented in Figure 7, extracted from WCGs $G^w$, among which three are acyclic, two correspond to an unfaithful distribution, and one is cyclic. The generating process of all datasets is the following: for all $Y$, for all $t > 0$,

$$Y_t = a_y Y_{t-1} + \sum_{X_{t-\ell} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-\ell} + 0.1 \xi_t^y,$$

where $a_y, a_{xy} \in U([-1, -0.1] \cup [0.1, 1])$ and for each parent $X$ we randomly choose if $X$ causes $Y$ instantaneously or with a lag of 1, i.e., $\ell \in \{0, 1\}$. Regarding the noise, we consider two different settings, in the first the noise is drawn from a uniform distribution, i.e., $\xi^y \sim U([-1, 1])$, and in the second the noise is drawn from a Gaussian distribution, i.e., $\xi^y \sim N(0, 1)$. For each setting and each structure in Figure 7 we generate 100 datasets of 1000 timestamps.

For the unfaithful diamond structure in Figure 7, following Zhalama et al. [2016], we set $a_x, a_y, a_z, a_w$ to zero, and $a_{zw} = -a_{xy}a_{yw}/a_{xz}$. All relationships are considered instantaneous. Consequently, in the distribution consistent with this model, we observe $X \perp W$, which violates the faithfulness assumption as this independence is not entailed by the
<table><caption>Table 1: Results obtained on the simulated data of Section 5.2 for the different structures with 1000 observations with non-Gaussian noise (top panel) and with Gaussian noise (bottom panel). We report the mean and the variance of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Diamond</th><th>Cyclic Diamond</th><th>Unf. Diamond</th><th>Adj. unf. Diamond</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">Non-Gaussian noise</td></tr><tr><td>NBCB-w</td><td><strong>0.94</strong> ± 0.01</td><td><strong>0.81</strong> ± 0.01</td><td>0.86 ± 0.01</td><td>0.8 ± 0.01</td></tr><tr><td>CBNB-w</td><td><strong>0.94</strong> ± 0.01</td><td><strong>0.8</strong> ± 0.01</td><td>0.9 ± 0.01</td><td>0.84 ± 0.01</td></tr><tr><td>NBCB-e</td><td>0.74 ± 0.02</td><td>0.72 ± 0.01</td><td>0.95 ± 0.01</td><td><strong>0.86</strong> ± 0.01</td></tr><tr><td>CBNB-e</td><td>0.74 ± 0.02</td><td>0.7 ± 0.02</td><td><strong>0.96</strong> ± 0.01</td><td><strong>0.86</strong> ± 0.01</td></tr><tr><td>GCMVL</td><td>0.86 ± 0.01</td><td>0.68 ± 0.01</td><td>0.04 ± 0.01</td><td>0.04 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td>0.92 ± 0.01</td><td>0.75 ± 0.01</td><td>0.47 ± 0.04</td><td>0.44 ± 0.03</td></tr><tr><td>PCGCE</td><td>0.69 ± 0.02</td><td>0.66 ± 0.01</td><td>0.5 ± 0.01</td><td>0.45 ± 0.01</td></tr><tr><td>Dynotears</td><td>0.03 ± 0.01</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td></tr><tr><td>VarLiNGAM</td><td><strong>0.99</strong> ± 0.01</td><td>0.79 ± 0.01</td><td><strong>0.98</strong> ± 0.01</td><td><strong>0.87</strong> ± 0.01</td></tr><tr><td colspan="5" style="text-align: center;">Gaussian noise</td></tr><tr><td>NBCB-w</td><td>0.78 ± 0.03</td><td><strong>0.77</strong> ± 0.01</td><td><strong>0.52</strong> ± 0.05</td><td><strong>0.48</strong> ± 0.01</td></tr><tr><td>CBNB-w</td><td>0.8 ± 0.04</td><td><strong>0.75</strong> ± 0.01</td><td><strong>0.52</strong> ± 0.05</td><td><strong>0.48</strong> ± 0.06</td></tr><tr><td>NBCB-e</td><td>0.64 ± 0.02</td><td>0.67 ± 0.02</td><td><strong>0.52</strong> ± 0.07</td><td><strong>0.44</strong> ± 0.06</td></tr><tr><td>CBNB-e</td><td>0.72 ± 0.03</td><td>0.65 ± 0.02</td><td><strong>0.53</strong> ± 0.07</td><td><strong>0.44</strong> ± 0.05</td></tr><tr><td>GCMVL</td><td><strong>0.87</strong> ± 0.01</td><td>0.7 ± 0.01</td><td>0.03 ± 0.01</td><td>0.01 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td><strong>0.93</strong> ± 0.01</td><td><strong>0.75</strong> ± 0.01</td><td>0.42 ± 0.05</td><td>0.4 ± 0.04</td></tr><tr><td>PCGCE</td><td>0.69 ± 0.02</td><td>0.65 ± 0.01</td><td>0.5 ± 0.02</td><td><strong>0.44</strong> ± 0.01</td></tr><tr><td>Dynotears</td><td>0.06 ± 0.02</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td></tr><tr><td>VarLiNGAM</td><td>0.78 ± 0.03</td><td>0.74 ± 0.01</td><td>0.5 ± 0.07</td><td>0.42 ± 0.06</td></tr></tbody></table>

causal Markov condition. However, it does not violate the adjacency faithfulness assumption since $X$ is not adjacent to $W$ in the graph. Similarly, for the adjacency unfaithful diamond structure, we set $a_{xw} = -a_{xy}a_{yw} - a_{xz}a_{zw}$, which violates the adjacency faithfulness assumption since in addition to $X$ being adjacent to $W$ in the graph, in the distribution, we observe $X \perp W$, and this independence is not entailed by the causal Markov condition.

In Table 1, we report the results for the setting with non-Gaussian noise (top panel), and with Gaussian noise (bottom panel).

Let us start with the case of non-Gaussian noise, which is easier to handle for most methods. For the first two structures, diamond and cyclic diamond, all assumptions required by both families of methods are satisfied, as well as by competitors. In this case, VarLiNGAM has the best performance for diamond, followed by CBNB-w and NBCB-w, while CBNB-w and NBCB-w have the best performance for cyclic diamond, followed by VarLiNGAM. PCMCI$^+$ comes closest to these methods. As expected from these results, we can see that CBNB-w and NBCB-w perform better than both original noise-based and constraint-based methods or as a trade-off between them. We can note that the same conclusion can be made for NBCB-e and CBNB-e. For the unfaithful structures, constraint-based methods PCMCI$^+$ and PCGCE have a drop in performance due to faithfulness violation. VarLiNGAM has the best results, which makes sense as it does not rely on faithfulness, followed by CBNB-e, with close results by NBCB-e. We can also note that our methods experience a noticeable drop in performance only for adjacency unfaithful structures, confirming that they do not require a full faithfulness assumption but are surprisingly still competitive for adjacency unfaithful structure, illustrating some robustness.

Considering Gaussian noise is more challenging, as the SCM is not an identifiable functional model (Assumption 4 is violated), which is needed for VarLiNGAM (and its use within the proposed methods). For the first diamond structure, as expected, PCMCI$^+$ performs best, closely followed by GCMVL and CBNB-w, NBCB-w, and VarLiNGAM. NBCB-e and CBNB-e demonstrate lower performance than VarLiNGAM, due to lower results of PCCGE. For the cyclic diamond structure, NBCB-w has the best results, followed by CBNB-w and PCMCI$^+$. For unfaithful structures with Gaussian noise, the constraint-based methods PCMCI$^+$ and PCGCE again experience a drop in performance due to faithfulness violation. Our methods yield the best results. Specifically, CBNB-e performs best for the unfaithful diamond, with the rest of our methods being the second-best. NBCB-w and CBNB-w work best for the adjacency unfaithful diamond, followed by NBCB-e, CBNB-e algorithms, and PCGCE.
<table><caption>Table 1: Results obtained on the simulated data of Section 5.2 for the different structures with 1000 observations with non-Gaussian noise (top panel) and with Gaussian noise (bottom panel). We report the mean and the variance of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Diamond</th><th>Cyclic Diamond</th><th>Unf. Diamond</th><th>Adj. unf. Diamond</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Non-Gaussian noise</td></tr><tr><td>NBCB-w</td><td><strong>0.94</strong> ± 0.01</td><td><strong>0.81</strong> ± 0.01</td><td>0.86 ± 0.01</td><td>0.8 ± 0.01</td></tr><tr><td>CBNB-w</td><td><strong>0.94</strong> ± 0.01</td><td><strong>0.8</strong> ± 0.01</td><td>0.9 ± 0.01</td><td>0.84 ± 0.01</td></tr><tr><td>NBCB-e</td><td>0.74 ± 0.02</td><td>0.72 ± 0.01</td><td>0.95 ± 0.01</td><td><strong>0.86</strong> ± 0.01</td></tr><tr><td>CBNB-e</td><td>0.74 ± 0.02</td><td>0.7 ± 0.02</td><td><strong>0.96</strong> ± 0.01</td><td><strong>0.86</strong> ± 0.01</td></tr><tr><td>GCMVL</td><td>0.86 ± 0.01</td><td>0.68 ± 0.01</td><td>0.04 ± 0.01</td><td>0.04 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td>0.92 ± 0.01</td><td>0.75 ± 0.01</td><td>0.47 ± 0.04</td><td>0.44 ± 0.03</td></tr><tr><td>PCGCE</td><td>0.69 ± 0.02</td><td>0.66 ± 0.01</td><td>0.5 ± 0.01</td><td>0.45 ± 0.01</td></tr><tr><td>Dynotears</td><td>0.03 ± 0.01</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td></tr><tr><td>VarLiNGAM</td><td><strong>0.99</strong> ± 0.01</td><td>0.79 ± 0.01</td><td><strong>0.98</strong> ± 0.01</td><td><strong>0.87</strong> ± 0.01</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Gaussian noise</td></tr><tr><td>NBCB-w</td><td>0.78 ± 0.03</td><td><strong>0.77</strong> ± 0.01</td><td><strong>0.52</strong> ± 0.05</td><td><strong>0.48</strong> ± 0.01</td></tr><tr><td>CBNB-w</td><td>0.8 ± 0.04</td><td><strong>0.75</strong> ± 0.01</td><td><strong>0.52</strong> ± 0.05</td><td><strong>0.48</strong> ± 0.06</td></tr><tr><td>NBCB-e</td><td>0.64 ± 0.02</td><td>0.67 ± 0.02</td><td><strong>0.52</strong> ± 0.07</td><td><strong>0.44</strong> ± 0.06</td></tr><tr><td>CBNB-e</td><td>0.72 ± 0.03</td><td>0.65 ± 0.02</td><td><strong>0.53</strong> ± 0.07</td><td><strong>0.44</strong> ± 0.05</td></tr><tr><td>GCMVL</td><td><strong>0.87</strong> ± 0.01</td><td>0.7 ± 0.01</td><td>0.03 ± 0.01</td><td>0.01 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td><strong>0.93</strong> ± 0.01</td><td><strong>0.75</strong> ± 0.01</td><td>0.42 ± 0.05</td><td>0.4 ± 0.04</td></tr><tr><td>PCGCE</td><td>0.69 ± 0.02</td><td>0.65 ± 0.01</td><td>0.5 ± 0.02</td><td><strong>0.44</strong> ± 0.01</td></tr><tr><td>Dynotears</td><td>0.06 ± 0.02</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td><td>0.0 ± 0.0</td></tr><tr><td>VarLiNGAM</td><td>0.78 ± 0.03</td><td>0.74 ± 0.01</td><td>0.5 ± 0.07</td><td>0.42 ± 0.06</td></tr></tbody></table>

causal Markov condition. However, it does not violate the adjacency faithfulness assumption since $X$ is not adjacent to $W$ in the graph. Similarly, for the adjacency unfaithful diamond structure, we set $a_{xw} = -a_{xy}a_{yw} - a_{xz}a_{zw}$, which violates the adjacency faithfulness assumption since in addition to $X$ being adjacent to $W$ in the graph, in the distribution, we observe $X \perp W$, and this independence is not entailed by the causal Markov condition.

In Table 1, we report the results for the setting with non-Gaussian noise (top panel), and with Gaussian noise (bottom panel).

Let us start with the case of non-Gaussian noise, which is easier to handle for most methods. For the first two structures, diamond and cyclic diamond, all assumptions required by both families of methods are satisfied, as well as by competitors. In this case, VarLiNGAM has the best performance for diamond, followed by CBNB-w and NBCB-w, while CBNB-w and NBCB-w have the best performance for cyclic diamond, followed by VarLiNGAM. PCMCI$^+$ comes closest to these methods. As expected from these results, we can see that CBNB-w and NBCB-w perform better than both original noise-based and constraint-based methods or as a trade-off between them. We can note that the same conclusion can be made for NBCB-e and CBNB-e. For the unfaithful structures, constraint-based methods PCMCI$^+$ and PCGCE have a drop in performance due to faithfulness violation. VarLiNGAM has the best results, which makes sense as it does not rely on faithfulness, followed by CBNB-e, with close results by NBCB-e. We can also note that our methods experience a noticeable drop in performance only for adjacency unfaithful structures, confirming that they do not require a full faithfulness assumption but are surprisingly still competitive for adjacency unfaithful structure, illustrating some robustness.

Considering Gaussian noise is more challenging, as the SCM is not an identifiable functional model (Assumption 4 is violated), which is needed for VarLiNGAM (and its use within the proposed methods). For the first diamond structure, as expected, PCMCI$^+$ performs best, closely followed by GCMVL and CBNB-w, NBCB-w, and VarLiNGAM. NBCB-e and CBNB-e demonstrate lower performance than VarLiNGAM, due to lower results of PCCGE. For the cyclic diamond structure, NBCB-w has the best results, followed by CBNB-w and PCMCI$^+$. For unfaithful structures with Gaussian noise, the constraint-based methods PCMCI$^+$ and PCGCE again experience a drop in performance due to faithfulness violation. Our methods yield the best results. Specifically, CBNB-e performs best for the unfaithful diamond, with the rest of our methods being the second-best. NBCB-w and CBNB-w work best for the adjacency unfaithful diamond, followed by NBCB-e, CBNB-e algorithms, and PCGCE.
<table><caption>Table 2: Results for realistic datasets of Section 5.3 generated using the Lotka–Volterra model with five species (left column) and ten species (right column). We report the mean and the variance of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in bold green.</caption><thead><tr><th></th><th>Lotka-Volterra(5)</th><th>Lotka-Volterra(10)</th></tr></thead><tbody><tr><td>NBCB-w</td><td>0.41 ± 0.03</td><td><b>0.28</b> ± 0.01</td></tr><tr><td>CBNB-w</td><td>0.38 ± 0.03</td><td><b>0.24</b> ± 0.01</td></tr><tr><td>NBCB-e</td><td><b>0.47</b> ± 0.02</td><td><b>0.24</b> ± 0.01</td></tr><tr><td>CBNB-e</td><td><b>0.44</b> ± 0.02</td><td>0.23 ± 0.01</td></tr><tr><td>GCMVL</td><td>0.19 ± 0.03</td><td>0.11 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td>0.36 ± 0.03</td><td>0.22 ± 0.01</td></tr><tr><td>PCGCE</td><td><b>0.44</b> ± 0.02</td><td>0.22 ± 0.01</td></tr><tr><td>Dynotears</td><td>0.18 ± 0.05</td><td>0.15 ± 0.01</td></tr><tr><td>VarLiNGAM</td><td>0.43 ± 0.06</td><td>0.23 ± 0.02</td></tr></tbody></table>

We highlight the consistently poor performance of Dynotears in Table 1, where this method has the lowest result in all scenarios. This can be attributed to the fact that in our simulated data, the variances do not increase in accordance with the topological order of the WCG.

Comparing PCMCI<sup>+</sup> and PCGCE, we can see that PCGCE has lower performance in general, except for the unfaithful cases for both types of noise distribution. This empirical observation suggests that PCGCE is more robust to assumption violation. This behavior is also inherited by CBNB-e and NBCB-e methods.

We can conclude that when necessary assumptions are satisfied, CBNB-w and NBCB-w are either trade-offs between the PCMCI<sup>+</sup> and VarLiNGAM or perform better, and NBCB-e and CBNB-e, are trade-offs between PCGCE and VarLiNGAM. For the unfaithful structures under assumption violation, all our methods are more robust compared to constraint and noise-based families.

Table 1 presents the F1 score of the orientations in the SCG, which is not suitable to illustrate the Proposition 2. Thus, in Table 6 in Appendix E we present the F1 score on the adjacencies, which illustrate the robustness of CBNB class to Assumption 4. More precisely, CBNB-w has the same performance as PCMCI<sup>+</sup> and better than the results of NBCB-w, in the case of Gaussian noise, for all structures, except the adjacency unfaithful diamond. We also see similar results for CBNB-e. Proposition 1 is difficult to see from Table 1, as we do not evaluate separately orientations and adjacencies, due to the specific structure of the inferred graph.

## 5.3 Realistic ecological data from the Lotka–Volterra model

We consider the simulation with multi-species generalization of the Ricker model introduced by Poggiato et al. [2022], which is analogous to the generalized Lotka–Volterra model with abiotic control presented in the same paper and is commonly used in ecological studies. Ricker model with abiotic control in discrete time for abundance of species Y at time step *t* has the following form:

$$ Y_t = \begin{cases} Y_{t-1} \exp \left( \Delta t \left( \sum_{X_{t-1} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-1} + \bar{Y}(-a_y) \exp \left( -\frac{(o_y - x)^2}{2\sigma_y^2} \right) \right) + \xi_t^y \right), & \text{(preys)} \\ Y_{t-1} \exp \left( \Delta t \left( \sum_{X_{t-1} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-1} - \mu \right) + \xi_t^y \right), & \text{(predators)} \end{cases} $$

where $Y_t$ is the abundance of species Y at time *t*, the upper equation is related to preys and the lower to predator species, $\bar{Y}$ is the abundance of species Y in the stationary state, $a_{xy}$ is the strength of the effect of species X on species Y and $a_y$ is strength of the effect of species Y on itself, $\xi_t^y$ is an i.i.d Gaussian random variable with variance $\sigma_r$, $o_y$ is a niche optimum for species Y, $x$ is the environmental variable, $\mu$ is the extinction rate of the predator. We run this simulation for the number of species $S = \{5, 10\}$ with the following fixed parameters: fixed environment $x = 0.5$, number of time steps $T = 1000$, $\mu = 0.05$, $\sigma_r = 0.2$, for each species $o_y$ randomly sampled from $U([0.05, 0, 95])$, the interaction matrix related to coefficients $a_{xy}$ and $a_y$ is obtained through a randomly generated WCG $G^w$ which is compatible with an SCG $G^s$. The SCG is constrained to contain only bi-directed edges and to encompass precisely 3 trophic levels, representing the hierarchical positions of species in the food chain. Specifically, these levels include basal species or prey ($L_1$), their predators ($L_2$), and the predators of predators ($L_3$). The second constraint ensures that $\forall X \in L_i$ and $\forall Y \in L_j$, if $X \rightleftharpoons Y$ in $G^s$ then $|i - j| = 1$ [for more details see Poggiato et al., 2022]. The interaction strength is randomly sampled for all interactions. We generate 100 graphs and thus we obtain 100 datasets.

For the Lotka-Volterra(5) datasets, NBCB-e performs better than the other methods followed by CBNB-e and PCGCE as shown in Table 2. Close to them perform VarLiNGAM, NBCB-w and CBNB-w. Dynotears and GCMVL have the
<table><caption>Table 2: Results for realistic datasets of Section 5.3 generated using the Lotka–Volterra model with five species (left column) and ten species (right column). We report the mean and the variance of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in bold green.</caption><thead><tr><th></th><th>Lotka–Volterra(5)</th><th>Lotka–Volterra(10)</th></tr></thead><tbody><tr><td>NBCB-w</td><td>0.41 ± 0.03</td><td><b>0.28</b> ± 0.01</td></tr><tr><td>CBNB-w</td><td>0.38 ± 0.03</td><td><b>0.24</b> ± 0.01</td></tr><tr><td>NBCB-e</td><td><b>0.47</b> ± 0.02</td><td><b>0.24</b> ± 0.01</td></tr><tr><td>CBNB-e</td><td><b>0.44</b> ± 0.02</td><td>0.23 ± 0.01</td></tr><tr><td>GCMVL</td><td>0.19 ± 0.03</td><td>0.11 ± 0.01</td></tr><tr><td>PCMCI<sup>+</sup></td><td>0.36 ± 0.03</td><td>0.22 ± 0.01</td></tr><tr><td>PCGCE</td><td><b>0.44</b> ± 0.02</td><td>0.22 ± 0.01</td></tr><tr><td>Dynotears</td><td>0.18 ± 0.05</td><td>0.15 ± 0.01</td></tr><tr><td>VarLiNGAM</td><td>0.43 ± 0.06</td><td>0.23 ± 0.02</td></tr></tbody></table>

We highlight the consistently poor performance of Dynotears in Table 1, where this method has the lowest result in all scenarios. This can be attributed to the fact that in our simulated data, the variances do not increase in accordance with the topological order of the WCG.

Comparing PCMCI<sup>+</sup> and PCGCE, we can see that PCGCE has lower performance in general, except for the unfaithful cases for both types of noise distribution. This empirical observation suggests that PCGCE is more robust to assumption violation. This behavior is also inherited by CBNB-e and NBCB-e methods.

We can conclude that when necessary assumptions are satisfied, CBNB-w and NBCB-w are either trade-offs between the PCMCI<sup>+</sup> and VarLiNGAM or perform better, and NBCB-e and CBNB-e, are trade-offs between PCGCE and VarLiNGAM. For the unfaithful structures under assumption violation, all our methods are more robust compared to constraint and noise-based families.

Table 1 presents the F1 score of the orientations in the SCG, which is not suitable to illustrate the Proposition 2. Thus, in Table 6 in Appendix E we present the F1 score on the adjacencies, which illustrate the robustness of CBNB class to Assumption 4. More precisely, CBNB-w has the same performance as PCMCI<sup>+</sup> and better than the results of NBCB-w, in the case of Gaussian noise, for all structures, except the adjacency unfaithful diamond. We also see similar results for CBNB-e. Proposition 1 is difficult to see from Table 1, as we do not evaluate separately orientations and adjacencies, due to the specific structure of the inferred graph.

## 5.3 Realistic ecological data from the Lotka–Volterra model

We consider the simulation with multi-species generalization of the Ricker model introduced by Poggiato et al. [2022], which is analogous to the generalized Lotka–Volterra model with abiotic control presented in the same paper and is commonly used in ecological studies. Ricker model with abiotic control in discrete time for abundance of species Y at time step t has the following form:

$$ Y_t = \begin{cases} Y_{t-1} \exp \left( \Delta t \left( \sum_{X_{t-1} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-1} + \bar{Y} (-a_y) \exp \left( -\frac{(o_y - x)^2}{2\sigma_y^2} \right) \right) + \xi_t^y \right), & (preys) \\ Y_{t-1} \exp \left( \Delta t \left( \sum_{X_{t-1} \in Pa_{G^w}(Y_t, G^w)} a_{xy} X_{t-1} - \mu \right) + \xi_t^y \right), & (predators) \end{cases} $$

where $Y_t$ is the abundance of species Y at time t, the upper equation is related to preys and the lower to predator species, $\bar{Y}$ is the abundance of species Y in the stationary state, $a_{xy}$ is the strength of the effect of species X on species Y and $a_y$ is strength of the effect of species Y on itself, $\xi_t^y$ is an i.i.d Gaussian random variable with variance $\sigma_r$, $o_y$ is a niche optimum for species Y, $x$ is the environmental variable, $\mu$ is the extinction rate of the predator. We run this simulation for the number of species $S = \{5, 10\}$ with the following fixed parameters: fixed environment $x = 0.5$, number of time steps $T = 1000$, $\mu = 0.05$, $\sigma_r = 0.2$, for each species $o_y$ randomly sampled from $U([0.05, 0, 95])$, the interaction matrix related to coefficients $a_{xy}$ and $a_y$ is obtained through a randomly generated WCG $G^w$ which is compatible with an SCG $G^s$. The SCG is constrained to contain only bi-directed edges and to encompass precisely 3 trophic levels, representing the hierarchical positions of species in the food chain. Specifically, these levels include basal species or prey ($L_1$), their predators ($L_2$), and the predators of predators ($L_3$). The second constraint ensures that $\forall X \in L_i$ and $\forall Y \in L_j$, if $X \rightleftharpoons Y$ in $G^s$ then $|i - j| = 1$ [for more details see Poggiato et al., 2022]. The interaction strength is randomly sampled for all interactions. We generate 100 graphs and thus we obtain 100 datasets.

For the Lotka–Volterra(5) datasets, NBCB-e performs better than the other methods followed by CBNB-e and PCGCE as shown in Table 2. Close to them perform VarLiNGAM, NBCB-w and CBNB-w. Dynotears and GCMVL have the
<table><caption>Table 3: Results for real datasets of Section 5.4 using linear methods. We report the mean and the variance (when meaningful, see data description) of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Temp.</th><th>Veil1</th><th>Veil2</th><th>Dairy</th><th>Ingest.</th><th>Web1</th><th>Web2</th><th>Antivirus1</th><th>Antivirus2</th></tr></thead><tbody><tr><td>NBCB-w</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.47 ± 0.03</td><td>0.2</td><td>0.23</td><td>0.13</td><td>0.3</td></tr><tr><td>CBNB-w</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.46 ± 0.11</td><td><b>0.24</b></td><td>0.29</td><td>0.18</td><td>0.18</td></tr><tr><td>NBCB-e</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td>0.5 ± 0.05</td><td><b>0.24</b></td><td><b>0.42</b></td><td><b>0.29</b></td><td><b>0.38</b></td></tr><tr><td>CBNB-e</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td><b>0.52</b> ± 0.06</td><td>0.15</td><td><b>0.38</b></td><td><b>0.33</b></td><td>0.27</td></tr><tr><td>GCMVL</td><td>0.67</td><td><b>1</b></td><td><b>1</b></td><td>0.33</td><td>0.5 ± 0.03</td><td>0.19</td><td>0.0</td><td>0.08</td><td>0.0</td></tr><tr><td>PCMCI<sup>+</sup></td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.3 ± 0.08</td><td>0.17</td><td>0.32</td><td>0.04</td><td>0.11</td></tr><tr><td>PCGCE</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td><b>0.55</b> ± 0.03</td><td>0.21</td><td>0.34</td><td><b>0.29</b></td><td><b>0.36</b></td></tr><tr><td>Dynotears</td><td>0.67</td><td><b>1</b></td><td><b>1</b></td><td>0.33</td><td>0.25 ± 0.06</td><td>0.22</td><td>0.3</td><td>0.18</td><td>0.17</td></tr><tr><td>VarLiNGAM</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.5</b></td><td>0.49 ± 0.05</td><td><b>0.23</b></td><td>0.2</td><td>0.18</td><td>0.18</td></tr></tbody></table>

lowest results. For the Lotka-Volterra(10) datasets, all results saw a significant decrease, with NBCB-w performing the best, followed by CBNB-w and NBCB-e. It is worth noting that for all datasets, NBCB-e and CBNB-e perform either better or equally as well as PCGCE and VarLiNGAM, while NBCB-w and CBNB-w outperform PCMCI.

## 5.4 Real data

Nine different real datasets are considered in this study. Taking into account the limitations of certain methods in handling nonlinearity, we start by evaluating our algorithms using linear tests and linear regressions, alongside linear tests for constraint-based methods. Then, we proceed to compare the nonlinear counterparts of our methods with those of the constraint-based methods while providing a computational time analysis.

### 5.4.1 The linear case

We detail the performance of each method in the following paragraphs, while the results are summarized in Table 3. See Appendix D for URL links to the considered datasets.

**Temperature.** This is a bivariate time series of length 168 about indoor *I* and outdoor *O* measurements. As noted by Assaad et al. [2021], it is expected that *O* causes *I*.

NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI<sup>+</sup>, PCGCE, and VarLiNGAM correctly infer $O \to I$. GCMVL and Dynotears infer a bidirected causal relation.

**Veilleux.** We considered two datasets for Figure 11(a) and 12(a) [Jost and Ellner, 2000] from Veilleux [1979] which study interactions between predatory ciliate Dinidum nasutum and its prey Paramecium aurelia with different values of Cerophyl concentrations (CC): 0.375 and 0.5. The lengths of the time series are 71 and 65. These data were previously analyzed with causal discovery algorithms [Barraquand et al., 2021, Sugihara et al., 2021], which showed bidirectional relationships in both cases.

Here, NBCB-e, CBNB-e, GCMVL, PCGCE, Dynotears and VarLiNGAM discover bidirected relationships Paramecium $\Leftrightarrow$ Didinium in both datasets, which is consistent with Sugihara et al. [2012], Barraquand et al. [2021]. CBNB-w, NBCB-w, PCMCI<sup>+</sup> have detected bidirected relationships between Paramecium and Didinium only in the first dataset.

**Diary.** This dataset provides ten years (from 09/2008 to 12/2018) of monthly prices for milk *M*, butter *B*, and cheddar cheese *C*, so the three time series are of length 124. We expect that the price of milk is a common cause of the price of butter and the price of cheddar cheese: $B \leftarrow M \to C$.

NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI<sup>+</sup>, PCGCE and VarLiNGAM correctly inferred that $M \to B$. But NBCB-w and CBNB-w wrongly inferred that $M \leftarrow C \to B$, NBCB-e, CBNB-e, wrongly inferred that $M \leftarrow B \leftarrow C$, PCMCI<sup>+</sup> and PCGCE wrongly inferred that $M \leftarrow C \leftarrow B$ and VarLiNGAM wrongly inferred that $C \to B$. GCMVL wrongly infers $B \Leftrightarrow M \leftarrow C \to B$ and Dynotears wrongly infers $M \Leftrightarrow B \Leftrightarrow C$.

**Ingestion mini.** This benchmark provided by EasyVista consists of 24 datasets each containing three time series with 1000 timestamps collected from an IT monitoring system with a one-minute sampling rate. Half the datasets are compatible with the graph C.M.I $\leftarrow$ M.E $\to$ G.H.I [Assaad et al., 2022b] where M.E is the metric extraction which
<table><caption>Table 3: Results for real datasets of Section 5.4 using linear methods. We report the mean and the variance (when meaningful, see data description) of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Temp.</th><th>Veil1</th><th>Veil2</th><th>Dairy</th><th>Ingest.</th><th>Web1</th><th>Web2</th><th>Antivirus1</th><th>Antivirus2</th></tr></thead><tbody><tr><td>NBCB-w</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.47 ± 0.03</td><td>0.2</td><td>0.23</td><td>0.13</td><td>0.3</td></tr><tr><td>CBNB-w</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.46 ± 0.11</td><td><b>0.24</b></td><td>0.29</td><td>0.18</td><td>0.18</td></tr><tr><td>NBCB-e</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td>0.5 ± 0.05</td><td><b>0.24</b></td><td><b>0.42</b></td><td><b>0.29</b></td><td><b>0.38</b></td></tr><tr><td>CBNB-e</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td><b>0.52</b> ± 0.06</td><td>0.15</td><td><b>0.38</b></td><td><b>0.33</b></td><td>0.27</td></tr><tr><td>GCMVL</td><td>0.67</td><td><b>1</b></td><td><b>1</b></td><td>0.33</td><td>0.5 ± 0.03</td><td>0.19</td><td>0.0</td><td>0.08</td><td>0.0</td></tr><tr><td>PCMCI<sup>+</sup></td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td>0.3 ± 0.08</td><td>0.17</td><td>0.32</td><td>0.04</td><td>0.11</td></tr><tr><td>PCGCE</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td><b>0.55</b> ± 0.03</td><td>0.21</td><td>0.34</td><td><b>0.29</b></td><td><b>0.36</b></td></tr><tr><td>Dynotears</td><td>0.67</td><td><b>1</b></td><td><b>1</b></td><td>0.33</td><td>0.25 ± 0.06</td><td>0.22</td><td>0.3</td><td>0.18</td><td>0.17</td></tr><tr><td>VarLiNGAM</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.5</b></td><td>0.49 ± 0.05</td><td><b>0.23</b></td><td>0.2</td><td>0.18</td><td>0.18</td></tr></tbody></table>

lowest results. For the Lotka-Volterra(10) datasets, all results saw a significant decrease, with NBCB-w performing the best, followed by CBNB-w and NBCB-e. It is worth noting that for all datasets, NBCB-e and CBNB-e perform either better or equally as well as PCGCE and VarLiNGAM, while NBCB-w and CBNB-w outperform PCMCI.

## 5.4 Real data

Nine different real datasets are considered in this study. Taking into account the limitations of certain methods in handling nonlinearity, we start by evaluating our algorithms using linear tests and linear regressions, alongside linear tests for constraint-based methods. Then, we proceed to compare the nonlinear counterparts of our methods with those of the constraint-based methods while providing a computational time analysis.

### 5.4.1 The linear case

We detail the performance of each method in the following paragraphs, while the results are summarized in Table 3. See Appendix D for URL links to the considered datasets.

**Temperature.** This is a bivariate time series of length 168 about indoor *I* and outdoor *O* measurements. As noted by Assaad et al. [2021], it is expected that *O* causes *I*.

NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI<sup>+</sup>, PCGCE, and VarLiNGAM correctly infer $O \to I$. GCMVL and Dynotears infer a bidirected causal relation.

**Veilleux.** We considered two datasets for Figure 11(a) and 12(a) [Jost and Ellner, 2000] from Veilleux [1979] which study interactions between predatory ciliate Dinidum nasutum and its prey Paramecium aurelia with different values of Cerophyl concentrations (CC): 0.375 and 0.5. The lengths of the time series are 71 and 65. These data were previously analyzed with causal discovery algorithms [Barraquand et al., 2021, Sugihara et al., 2021], which showed bidirectional relationships in both cases.

Here, NBCB-e, CBNB-e, GCMVL, PCGCE, Dynotears and VarLiNGAM discover bidirected relationships Paramecium $\Leftrightarrow$ Didinium in both datasets, which is consistent with Sugihara et al. [2012], Barraquand et al. [2021]. CBNB-w, NBCB-w, PCMCI<sup>+</sup> have detected bidirected relationships between Paramecium and Didinium only in the first dataset.

**Diary.** This dataset provides ten years (from 09/2008 to 12/2018) of monthly prices for milk *M*, butter *B*, and cheddar cheese *C*, so the three time series are of length 124. We expect that the price of milk is a common cause of the price of butter and the price of cheddar cheese: $B \leftarrow M \to C$.

NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI<sup>+</sup>, PCGCE and VarLiNGAM correctly inferred that $M \to B$. But NBCB-w and CBNB-w wrongly inferred that $M \leftarrow C \to B$, NBCB-e, CBNB-e, wrongly inferred that $M \leftarrow B \leftarrow C$, PCMCI<sup>+</sup> and PCGCE wrongly inferred that $M \leftarrow C \leftarrow B$ and VarLiNGAM wrongly inferred that $C \to B$. GCMVL wrongly infers $B \Leftrightarrow M \leftarrow C \to B$ and Dynotears wrongly infers $M \Leftrightarrow B \Leftrightarrow C$.

**Ingestion mini.** This benchmark provided by EasyVista consists of 24 datasets each containing three time series with 1000 timestamps collected from an IT monitoring system with a one-minute sampling rate. Half the datasets are compatible with the graph C.M.I $\leftarrow$ M.E $\to$ G.H.I [Assaad et al., 2022b] where M.E is the metric extraction which
<figure><img src="image_15.png" /><figcaption>Figure 8: Summary causal graphs for (a) the Web activity datasets and (b) the Antivirus activity datasets. Those summary causal graphs are constructed by an IT monitoring system's experts.</figcaption></figure>

represents the activity of the extraction of the metrics from the messages; *G.H.I* is the group history insertion, which represents the activity of the insertion of the historical status in the database; and *C.M.I* is the collector monitoring information, which represents the activity of the updates in a given database. The second half of the datasets are compatible with the graph *M.D* → *M.E* → *M.I* where *M.D* is the metric dispatcher which represents the activity of a process that orients messages to other processes with respect to different types of messages; and *M.I* is the metric insertion which represents the activity of insertion of data in a database. Lags between time series are unknown, as well as the existence of self-causes.

From Table 3, we can see that PCGCE has the best results followed by CBNB-e then by NBCB-e and GCMVL. After that comes VarLiNGAM and NBCB-w and CBNB-w and PCMCI$^+$.

Finally, Dynotears has the worst result. This might suggest that the lags between causes and effects are not consistent over time, in the sense that if the lags between two-time series vary (while respecting the maximal temporal lag), the extended summary causal graph might remain the same however this is not true for the window causal graph. In this case, we might expect that methods inferring window causal graphs (such as PCMCI$^+$) would perform worse than methods inferring extended summary causal graphs (such as PCGCE).

**Web.** We consider a dataset that reflects the activity in a web server which is provided by EasyVista. This dataset contains ten time series collected with a one-minute sampling rate. The raw data of this case study were initially misaligned. In order to align them, we use the two pre-processing strategies described in Appendix C. We denote the dataset pre-processed using Strategy 1 as Web 1 and the dataset pre-processed using Strategy 2 as Web 2. The two processed datasets contain 3000 timestamps. The corresponding summary causal graph for Web dataset is presented in Figure 8a where *NetIn* represents the data received by the network interface card in Kbytes/second; *NetOut* represents the data transmitted out by the network interface card in Kbytes/second; *NPH* represents the number of HTTP processes; *NPP* represents the number of PHP processes; *NCM* represents the number of open MySQL connections which are started by PHP processes; *CpuH* represents the percentage of CPU used by all HTTP processes; *RamH* represents the percentage of RAM used by all HTTP processes; *CpuP* represents the percentage of CPU used by all PHP processes; *DiskW* represents the Disk write in Kbytes/second; *CpuG* represents the percentage of global CPU usage.

From Table 3, for the Web1 dataset, we can see that NBCB-e and CBNB-w demonstrate the highest performance, followed by VarLiNGAM. Next are Dynotears and PCGCE. CBNB-e exhibits the lowest score. For the Web2 dataset, NBCB-e attains the highest F1 score, with CBNB-e closely following. The subsequent competitive results are observed with PCGCE and PCMCI$^+$. This could imply that the data is noisy and the time lags are inconsistent throughout time, making the inference of the extended summary graph more robust.

**Antivirus.** Lastly, we consider a dataset which depicts the impacts of antivirus activity in servers which is again provided by EasyVista. This dataset contains 13 time series such that 3 of them are collected with a one-minute sampling rate and the rest with a five-minute sampling rate. The raw data of this case study were initially misaligned. To align them, we use the two pre-processing strategies described in Appendix C, leading to the dataset Antivirus 1 for Strategy 1 and Antivirus 2 for Strategy 2. The two processed datasets consist of 1321 timestamps. The corresponding
<figure><img src="image_15.png" /><figcaption>Figure 8: Summary causal graphs for (a) the Web activity datasets and (b) the Antivirus activity datasets. Those summary causal graphs are constructed by an IT monitoring system's experts.</figcaption></figure>

represents the activity of the extraction of the metrics from the messages; *G.H.I* is the group history insertion, which represents the activity of the insertion of the historical status in the database; and *C.M.I* is the collector monitoring information, which represents the activity of the updates in a given database. The second half of the datasets are compatible with the graph *M.D* → *M.E* → *M.I* where *M.D* is the metric dispatcher which represents the activity of a process that orients messages to other processes with respect to different types of messages; and *M.I* is the metric insertion which represents the activity of insertion of data in a database. Lags between time series are unknown, as well as the existence of self-causes.

From Table 3, we can see that PCGCE has the best results followed by CBNB-e then by NBCB-e and GCMVL. After that comes VarLiNGAM and NBCB-w and CBNB-w and PCMCI$^+$.

Finally, Dynotears has the worst result. This might suggest that the lags between causes and effects are not consistent over time, in the sense that if the lags between two-time series vary (while respecting the maximal temporal lag), the extended summary causal graph might remain the same however this is not true for the window causal graph. In this case, we might expect that methods inferring window causal graphs (such as PCMCI$^+$) would perform worse than methods inferring extended summary causal graphs (such as PCGCE).

**Web.** We consider a dataset that reflects the activity in a web server which is provided by EasyVista. This dataset contains ten time series collected with a one-minute sampling rate. The raw data of this case study were initially misaligned. In order to align them, we use the two pre-processing strategies described in Appendix C. We denote the dataset pre-processed using Strategy 1 as Web 1 and the dataset pre-processed using Strategy 2 as Web 2. The two processed datasets contain 3000 timestamps. The corresponding summary causal graph for Web dataset is presented in Figure 8a where *NetIn* represents the data received by the network interface card in Kbytes/second; *NetOut* represents the data transmitted out by the network interface card in Kbytes/second; *NPH* represents the number of HTTP processes; *NPP* represents the number of PHP processes; *NCM* represents the number of open MySQL connections which are started by PHP processes; *CpuH* represents the percentage of CPU used by all HTTP processes; *RamH* represents the percentage of RAM used by all HTTP processes; *CpuP* represents the percentage of CPU used by all PHP processes; *DiskW* represents the Disk write in Kbytes/second; *CpuG* represents the percentage of global CPU usage.

From Table 3, for the Web1 dataset, we can see that NBCB-e and CBNB-w demonstrate the highest performance, followed by VarLiNGAM. Next are Dynotears and PCGCE. CBNB-e exhibits the lowest score. For the Web2 dataset, NBCB-e attains the highest F1 score, with CBNB-e closely following. The subsequent competitive results are observed with PCGCE and PCMCI$^+$. This could imply that the data is noisy and the time lags are inconsistent throughout time, making the inference of the extended summary graph more robust.

**Antivirus.** Lastly, we consider a dataset which depicts the impacts of antivirus activity in servers which is again provided by EasyVista. This dataset contains 13 time series such that 3 of them are collected with a one-minute sampling rate and the rest with a five-minute sampling rate. The raw data of this case study were initially misaligned. To align them, we use the two pre-processing strategies described in Appendix C, leading to the dataset Antivirus 1 for Strategy 1 and Antivirus 2 for Strategy 2. The two processed datasets consist of 1321 timestamps. The corresponding
<table><caption>Table 4: Results for real datasets of Section 5.4 using nonlinear methods. We report the mean and the variance (when meaningful, see data description) of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Temp.</th><th>Veil1</th><th>Veil2</th><th>Dairy</th><th>Ingest.</th></tr></thead><tbody><tr><td>NBCB-w-nl</td><td>1</td><td>1</td><td>0</td><td>0.4</td><td>0.52 ± 0.02</td></tr><tr><td>CBNB-w-nl</td><td>1</td><td>1</td><td>0</td><td>0.4</td><td>0.54 ± 0.06</td></tr><tr><td>NBCB-e-nl</td><td>1</td><td>1</td><td>1</td><td>0.4</td><td>0.43 ± 0.04</td></tr><tr><td>CBNB-e-nl</td><td>1</td><td>1</td><td>1</td><td>0.4</td><td>0.49 ± 0.07</td></tr><tr><td>PCMCI<sup>+</sup>-nl</td><td>1</td><td>1</td><td>0</td><td>0.0</td><td>0.38 ± 0.09</td></tr><tr><td>PCGCE-nl</td><td>0.67</td><td>1</td><td>1</td><td>0.4</td><td>0.43 ± 0.03</td></tr></tbody></table>

summary causal graph for Antivirus dataset is presented in Figure 8b where CUV represents the percentage of CPU usage of antivirus processes in server V; CUGV represents the percentage of CPU usage of the global server V; MUV represents the percentage of memory usage of antivirus process; MUGV represents the percentage of global memory usage of the server; RV represents the Disk IO read in Kbytes/second; ChIE refers to the required duration in seconds to open an *IE browser* on server V; CUP represents the percentage of CPU usage of antivirus processes in server P; CUGP represents the percentage of CPU usage of the global server P; MUP represents the percentage of memory usage of antivirus process; MUGP represents the percentage of global memory usage of the server; RP represents the Disk IO read in Kbytes/second; ChP represents the required duration in seconds to open a CITRIX Portal on server P; T represents the global time in seconds required to open a CITRIX portal and open the IE browser.

From Table 3 for Antivirus 1, we can see that CBNB-e achieves the best result followed by PCGCE and NBCB-e and then by CBNB-w, VarLiNGAM and Dynotears. GCMVL and PCMCI⁺ have low performance, with PCMCI⁺ being the worst. For Antivirus 2, we can see that NBCB-e performs best followed by PCGCE and then by CBNB-e. The remaining methods have a substantial drop in performance. As in the case of Web datasets, improved performance of the PCGCE over PCMCI⁺ suggests that inference of the extended summary graph is more robust in this case. We can also note that NBCB-e and NBCB-w perform well, while VarLiNGAM does not, which could suggest that VarLiNGAM infers a more dense graph.

### 5.4.2 The nonlinear case

The results of the nonlinear counterparts of our hybrid methods and of PCMCI⁺ and PCGCE on real data are presented in Table 4. In this scenario, we excluded the Web and Antivirus datasets due to computational constraints (running the nonlinear counterparts of the algorithms on these datasets is prohibitively expensive due to their size).

In the temperature dataset, the Veilleux datasets and the Dairy dataset, all considered methods produced consistent results compared to their linear counterparts, except for PCGCE, which exhibited a decrease in performance in the temperature dataset. In the Ingestion mini datasets, we can clearly see that all methods that infer a WCG has a slight increase in performance and all methods that infer an ECG has a slight decrease in performance. In terms of ranking, in the nonlinear case, NBCB-w-nl demonstrates the best performance, followed by CBNB-w-nl, while PCMCI⁺ performs the worst despite its performance increase. The performance improvement observed in methods inferring a WCG can suggest that the dataset Ingestion mini contains nonlinear causal relations. Conversely, the decline seen in methods inferring an ECG could be attributed to the combination of the complexity of nonlinear tests and the necessity for these methods to conduct conditional independence tests with larger conditional sets compared to those considered by other methods.

We also conduct a computation time analysis of the nonlinear counterparts of the methods, as shown in Figure 9. Notably, the computation times for all methods are comparable across the temperature dataset, the Veilleux datasets, and the Dairy dataset, which are relatively small in size. However, in the Ingestion mini datasets, variations in computation time are evident: PCMCI⁺ exhibits the longest computation time, followed by CBNB-w then by PCGCE and NBCB-w. Our hybrid methods for inferring an ECG (NBCB-e and CBNB-e) exhibit the shortest computation time, with NBCB-e showing a lower computation time compared to CBNB-e. The NBCB methods achieve better time computation because the nonlinear conditional tests they employ are significantly more computationally expensive than learning Gaussian process regression. This implies that the constraint-based step is more costly than the noise-based step. Hence, providing additional knowledge to the constraint-based step, in the form of a causal order, results in a much greater reduction in computation time compared to when the additional knowledge is given to the noise-based step, in the form of the skeleton.
<table><caption>Table 4: Results for real datasets of Section 5.4 using nonlinear methods. We report the mean and the variance (when meaningful, see data description) of the F1 score of the orientations in the SCG. The best results are in blue bold and the second best results are in green bold.</caption><thead><tr><th></th><th>Temp.</th><th>Veil1</th><th>Veil2</th><th>Dairy</th><th>Ingest.</th></tr></thead><tbody><tr><td>NBCB-w-nl</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td><b>0.52</b> ± 0.02</td></tr><tr><td>CBNB-w-nl</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td><b>0.4</b></td><td><b>0.54</b> ± 0.06</td></tr><tr><td>NBCB-e-nl</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td>0.43 ± 0.04</td></tr><tr><td>CBNB-e-nl</td><td><b>1</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td>0.49 ± 0.07</td></tr><tr><td>PCMCI<sup>+</sup>-nl</td><td><b>1</b></td><td><b>1</b></td><td>0</td><td>0.0</td><td>0.38 ± 0.09</td></tr><tr><td>PCGCE-nl</td><td><b>0.67</b></td><td><b>1</b></td><td><b>1</b></td><td><b>0.4</b></td><td>0.43 ± 0.03</td></tr></tbody></table>

summary causal graph for Antivirus dataset is presented in Figure 8b where CUV represents the percentage of CPU usage of antivirus processes in server V; CUGV represents the percentage of CPU usage of the global server V; MUV represents the percentage of memory usage of antivirus process; MUGV represents the percentage of global memory usage of the server; RV represents the Disk IO read in Kbytes/second; ChIE refers to the required duration in seconds to open an *IE browser* on server V; CUP represents the percentage of CPU usage of antivirus processes in server P; CUGP represents the percentage of CPU usage of the global server P; MUP represents the percentage of memory usage of antivirus process; MUGP represents the percentage of global memory usage of the server; RP represents the Disk IO read in Kbytes/second; ChP represents the required duration in seconds to open a CITRIX Portal on server P; T represents the global time in seconds required to open a CITRIX portal and open the IE browser.

From Table 3 for Antivirus 1, we can see that CBNB-e achieves the best result followed by PCGCE and NBCB-e and then by CBNB-w, VarLiNGAM and Dynotears. GCMVL and PCMCI⁺ have low performance, with PCMCI⁺ being the worst. For Antivirus 2, we can see that NBCB-e performs best followed by PCGCE and then by CBNB-e. The remaining methods have a substantial drop in performance. As in the case of Web datasets, improved performance of the PCGCE over PCMCI⁺ suggests that inference of the extended summary graph is more robust in this case. We can also note that NBCB-e and NBCB-w perform well, while VarLiNGAM does not, which could suggest that VarLiNGAM infers a more dense graph.

### 5.4.2 The nonlinear case

The results of the nonlinear counterparts of our hybrid methods and of PCMCI⁺ and PCGCE on real data are presented in Table 4. In this scenario, we excluded the Web and Antivirus datasets due to computational constraints (running the nonlinear counterparts of the algorithms on these datasets is prohibitively expensive due to their size).

In the temperature dataset, the Veilleux datasets and the Dairy dataset, all considered methods produced consistent results compared to their linear counterparts, except for PCGCE, which exhibited a decrease in performance in the temperature dataset. In the Ingestion mini datasets, we can clearly see that all methods that infer a WCG has a slight increase in performance and all methods that infer an ECG has a slight decrease in performance. In terms of ranking, in the nonlinear case, NBCB-w-nl demonstrates the best performance, followed by CBNB-w-nl, while PCMCI⁺ performs the worst despite its performance increase. The performance improvement observed in methods inferring a WCG can suggest that the dataset Ingestion mini contains nonlinear causal relations. Conversely, the decline seen in methods inferring an ECG could be attributed to the combination of the complexity of nonlinear tests and the necessity for these methods to conduct conditional independence tests with larger conditional sets compared to those considered by other methods.

We also conduct a computation time analysis of the nonlinear counterparts of the methods, as shown in Figure 9. Notably, the computation times for all methods are comparable across the temperature dataset, the Veilleux datasets, and the Dairy dataset, which are relatively small in size. However, in the Ingestion mini datasets, variations in computation time are evident: PCMCI⁺ exhibits the longest computation time, followed by CBNB-w then by PCGCE and NBCB-w. Our hybrid methods for inferring an ECG (NBCB-e and CBNB-e) exhibit the shortest computation time, with NBCB-e showing a lower computation time compared to CBNB-e. The NBCB methods achieve better time computation because the nonlinear conditional tests they employ are significantly more computationally expensive than learning Gaussian process regression. This implies that the constraint-based step is more costly than the noise-based step. Hence, providing additional knowledge to the constraint-based step, in the form of a causal order, results in a much greater reduction in computation time compared to when the additional knowledge is given to the noise-based step, in the form of the skeleton.
<figure><img src="image_16.png" /><figcaption>Figure 9: Time computation (in second) for NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI$^{+}$, and PCGCE for real datasets of Section 5.4. We report the mean and the standard deviation.</figcaption></figure>

# 6 Discussion

Experiments on simulated data, realistic ecological data, and real data from various applications, show that our hybrid approaches are robust and yield overall good results over all datasets. Notice that for all results on real data, NBCB-w and CBNB-w (which are based on PCMCI$^{+}$ and VarLiNGAM) never perform simultaneously worse than PCMCI$^{+}$ and VarLiNGAM. Similarly, NBCB-e and CBNB-e (which are based on PCGCE and VarLiNGAM) never perform simultaneously worse than PCGCE and VarLiNGAM, except for the Web 1 dataset where CBNB-e has the lowest F1-score (but NBCB-e has the best F1-score). In general, NBCB-e and CCBB-e seem to be more reliable than NBCB-w and CCBB-w for the real data we considered, especially when assuming linearity. As mentioned before, the possible explanation is noisy data and inconsistent time lags. In summary, results on simulated data, realistic ecological data, and real data are coherent with the theoretical findings, showing that algorithms from CBNB and NBCB classes are trade-offs between the original methods, potentially exhibiting enhanced performance compared to the original methods when certain assumptions are violated.

<table><caption>Table 5: Summary of the two proposed classes of methods, NBCB and CBNB. Rows respectively indicate: their versions NBCB-e, NBCB-w, CBNB-e, CBNB-w; outputs; the algorithms used in each step; their theoretical guarantees, or lack thereof (X), under violation of some assumption; and their advantageous (✓) or limited (X) performances on the different data scenarios considered in Section 5. Ĝ* represents the inferred WCG or ECG and G* represents the true WCG or ECG.</caption><thead><tr><th rowspan="2">Class</th><th colspan="2">NBCB</th><th colspan="2">CBNB</th></tr><tr><th>NBCB-e<br>ECG</th><th>NBCB-w<br>WCG</th><th>CBNB-e<br>ECG</th><th>CBNB-w<br>WCG</th></tr></thead><tbody><tr><td>Version</td><td>NBCB-e</td><td>NBCB-w</td><td>CBNB-e</td><td>CBNB-w</td></tr><tr><td>Output</td><td>ECG</td><td>WCG</td><td>ECG</td><td>WCG</td></tr><tr><td>Step 1</td><td colspan="2">NB1 (VLINGAM)</td><td colspan="2">CB1 (PCGCE - PCMCI<sup>+</sup>)</td></tr><tr><td>Step 2</td><td colspan="2">CB1' (PCGCE - PCMCI<sup>+</sup>)</td><td colspan="2">NB1' (VLINGAM)</td></tr><tr><td>Violation of Assump. 3</td><td>(X<sub>t</sub> → Y<sub>t</sub> ∈ Ĝ<sup>*</sup> ⇒ (Y<sub>t</sub> ∝ X<sub>t</sub> ∈ G<sup>*</sup>)</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>Violation of Assump. 4</td><td>✗</td><td>✗</td><td>(X<sub>t*</sub> - Y<sub>t</sub> ∈ Ĝ<sup>*</sup>) ⇒ (X<sub>t*</sub> - Y<sub>t</sub> ∈ G<sup>*</sup>)</td><td>✗</td></tr><tr><td>Simulated data (Sec. 5.2)</td><td>✗</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Realistic sim. data (Sec. 5.3)</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Real data (Sec. 5.4)</td><td>✓</td><td>✗</td><td>✓</td><td>✗</td></tr></tbody></table>

In Table 5, we provide various theoretical and experimental criteria to distinguish between the algorithms within the NBCB and CBNB classes. In the second and third rows of Table 5, we present which algorithms infer a WCG and which ones infer an ECG. In the fourth and fifth rows, we detail the steps of the NBCB and CBNB classes of methods. The sixth row indicates that NBCB-w and NBCB-e can still be applied even when Assumption 3 is violated. In such cases, the true graph may not be fully retrieved, but it is guaranteed that if the algorithms infer X → Y, then in the true graph, we can be certain that Y does not cause X (see Proposition 1). Similarly, the seventh row shows that CBNB-w and CBNB-e can be utilized even when Assumption 4 is violated. In this scenario, algorithms are capable of inferring accurate skeleton of the graph, but not the orientations. The last three rows highlight the scenarios where each algorithm outperformed others. Interestingly, in the experimental section, we observed that methods inferring a WCG perform better with simulated data, while methods inferring an ECG excel with real data, assuming linearity.
<figure><img src="image_16.png" /><figcaption>Figure 9: Time computation (in second) for NBCB-w, CBNB-w, NBCB-e, CBNB-e, PCMCI$^{+}$, and PCGCE for real datasets of Section 5.4. We report the mean and the standard deviation.</figcaption></figure>

# 6 Discussion

Experiments on simulated data, realistic ecological data, and real data from various applications, show that our hybrid approaches are robust and yield overall good results over all datasets. Notice that for all results on real data, NBCB-w and CBNB-w (which are based on PCMCI$^{+}$ and VarLiNGAM) never perform simultaneously worse than PCMCI$^{+}$ and VarLiNGAM. Similarly, NBCB-e and CBNB-e (which are based on PCGCE and VarLiNGAM) never perform simultaneously worse than PCGCE and VarLiNGAM, except for the Web 1 dataset where CBNB-e has the lowest F1-score (but NBCB-e has the best F1-score). In general, NBCB-e and CCBB-e seem to be more reliable than NBCB-w and CCBB-w for the real data we considered, especially when assuming linearity. As mentioned before, the possible explanation is noisy data and inconsistent time lags. In summary, results on simulated data, realistic ecological data, and real data are coherent with the theoretical findings, showing that algorithms from CBNB and NBCB classes are trade-offs between the original methods, potentially exhibiting enhanced performance compared to the original methods when certain assumptions are violated.

<table><caption>Table 5: Summary of the two proposed classes of methods, NBCB and CBNB. Rows respectively indicate: their versions NBCB-e, NBCB-w, CBNB-e, CBNB-w; outputs; the algorithms used in each step; their theoretical guarantees, or lack thereof (X), under violation of some assumption; and their advantageous (✓) or limited (X) performances on the different data scenarios considered in Section 5. Ĝ* represents the inferred WCG or ECG and G* represents the true WCG or ECG.</caption><thead><tr><th rowspan="2">Class</th><th colspan="2">NBCB</th><th colspan="2">CBNB</th></tr><tr><th>NBCB-e<br>ECG</th><th>NBCB-w<br>WCG</th><th>CBNB-e<br>ECG</th><th>CBNB-w<br>WCG</th></tr></thead><tbody><tr><td>Version</td><td>NBCB-e</td><td>NBCB-w</td><td>CBNB-e</td><td>CBNB-w</td></tr><tr><td>Output</td><td>ECG</td><td>WCG</td><td>ECG</td><td>WCG</td></tr><tr><td>Step 1</td><td colspan="2">NB1 (VLINGAM)</td><td colspan="2">CB1 (PCGCE - PCMCI<sup>+</sup>)</td></tr><tr><td>Step 2</td><td colspan="2">CB1' (PCGCE - PCMCI<sup>+</sup>)</td><td colspan="2">NB1' (VLINGAM)</td></tr><tr><td>Violation of Assump. 3</td><td>(X<sub>t</sub> → Y<sub>t</sub> ∈ Ĝ<sup>*</sup> ⇒ (Y<sub>t</sub> ⊢ X<sub>t</sub> ∈ G<sup>*</sup>)</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>Violation of Assump. 4</td><td>✗</td><td>✗</td><td>(X<sub>t*</sub> - Y<sub>t</sub> ∈ Ĝ<sup>*</sup>) ⇒ (X<sub>t*</sub> - Y<sub>t</sub> ∈ G<sup>*</sup>)</td><td>✗</td></tr><tr><td>Simulated data (Sec. 5.2)</td><td>✗</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Realistic sim. data (Sec. 5.3)</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Real data (Sec. 5.4)</td><td>✓</td><td>✗</td><td>✓</td><td>✗</td></tr></tbody></table>

In Table 5, we provide various theoretical and experimental criteria to distinguish between the algorithms within the NBCB and CBNB classes. In the second and third rows of Table 5, we present which algorithms infer a WCG and which ones infer an ECG. In the fourth and fifth rows, we detail the steps of the NBCB and CBNB classes of methods. The sixth row indicates that NBCB-w and NBCB-e can still be applied even when Assumption 3 is violated. In such cases, the true graph may not be fully retrieved, but it is guaranteed that if the algorithms infer X → Y, then in the true graph, we can be certain that Y does not cause X (see Proposition 1). Similarly, the seventh row shows that CBNB-w and CBNB-e can be utilized even when Assumption 4 is violated. In this scenario, algorithms are capable of inferring accurate skeleton of the graph, but not the orientations. The last three rows highlight the scenarios where each algorithm outperformed others. Interestingly, in the experimental section, we observed that methods inferring a WCG perform better with simulated data, while methods inferring an ECG excel with real data, assuming linearity.
One of the key limitations of the CBNB and NBCB classes for real-world applications is their reliance on the restrictive assumption that there are no hidden confounders (Assumption 1), which could be often violated. For CBNB, we could consider an extension of the FRITL algorithm [Chen et al., 2021] for time-series, which is based on building the skeleton using FCI extension and refining after using noise-based methods. However, it is not clear how to adapt the necessary conditions for undirected cycle groups. For NBCB it could be more complicated to relax the no hidden confounders assumption.

Additionally, it would be interesting to adapt the CBNB and NBCB classes to the cases when there is a violation of consistency over time or stationarity. One potential direction could be to integrate our methods in the strategy proposed by Saggioro et al. [2020] which combines a causal discovery with a regime learning optimisation approach. However, this direction might require additional assumptions and while it appears promising for CBNB, it is less obvious for NBCB. Another direction could be to assume the presence of an observed contextual variable that explains the non-stationarity [Mooij et al., 2020, Günther et al., 2023].

Finally, adapting these methods for mixed data can be important in many applications. This can be straightforward for the constraint-based part of our algorithms, given the availability of conditional independence tests for mixed data [Zan et al., 2022]. However, the adaptation is more challenging for the noise-based part.

# 7 Conclusion

In this paper, we introduced a framework for hybrids of noise-based and constraint-based methods that can discover causal graphs from temporal data. Algorithms in the first class, denoted NBCB, start with ordering instantaneous relations, and then prune edges of the fully oriented graph. On the other hand, the algorithms from the second class, denoted CBNB, start by finding the skeleton and the orientation of lagged relations using a constraint-based method and temporal priority, then orient instantaneous relations by ordering nodes in each cycle group of instantaneous relations. Overall, the performance of our algorithms is a trade-off between the performance of constraint-based and noise-based algorithms when all assumptions are satisfied, and they outperform other methods in the cases where some of the assumptions are violated.

For future works, it would be interesting to extend these approaches to cases involving hidden confounders, non-stationarity, and mixed data.

# Acknowledgments

We thank Ali Aït-Bachir, Christophe de Bignicourt and Rachid Mokhtari from EasyVista for providing the IT monitoring data along with the underlying causal graphs. We also thank Giovanni Poggiato for several discussions about the realistic ecological simulated data. Finally, we thank the anonymous reviewers for their helpful comments. This work was partially supported by the CIPHOD project (ANR-23-CPJ1-0212-01), by MIAI@Grenoble Alpes (ANR-19-P3IA-0003), by the Horizon Europe Obsession project (No: 101134954) and the FRB-CESAB through the IMPACT working group.

# References

Hugo Arlegui, Emmanuel Mahé, M.-A Richard, Yann Rycke, M. Viguier, Marie Beylot-Barry, A. Dupuy, N. Beneton, P. Joly, Denis Jullien, H. Bachelez, É Sbidian, O. Chosidow, Carle Paul, and F. Tubach. Impact of the first wave of the covid-19 pandemic on the treatment of psoriasis with systemic therapies in france: Results from the psobioteq cohort. *Annales de Dermatologie et de Vénéréologie*, 150, 02 2023. doi: 10.1016/j.annder.2023.01.005.

Andrew Arnold, Yan Liu, and Naoki Abe. Temporal causal modeling with graphical granger methods. In *Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, KDD '07, page 66–75, New York, NY, USA, 2007. Association for Computing Machinery. ISBN 9781595936097. doi: 10.1145/1281192.1281203.

Charles K. Assaad, Emilie Devijver, Eric Gaussier, and Ali Aït-Bachir. A mixed noise and constraint-based approach to causal inference in time series. In *Machine Learning and Knowledge Discovery in Databases. Research Track*, pages 453–468, Cham, 2021. Springer International Publishing. ISBN 978-3-030-86486-6.

Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. *J. Artif. Int. Res.*, 73, may 2022a. doi: 10.1613/jair.1.13428.

Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Entropy-based discovery of summary causal graphs in time series. *Entropy*, 24(8), 2022b. ISSN 1099-4300. doi: 10.3390/e24081156.
One of the key limitations of the CBNB and NBCB classes for real-world applications is their reliance on the restrictive assumption that there are no hidden confounders (Assumption 1), which could be often violated. For CBNB, we could consider an extension of the FRITL algorithm [Chen et al., 2021] for time-series, which is based on building the skeleton using FCI extension and refining after using noise-based methods. However, it is not clear how to adapt the necessary conditions for undirected cycle groups. For NBCB it could be more complicated to relax the no hidden confounders assumption.

Additionally, it would be interesting to adapt the CBNB and NBCB classes to the cases when there is a violation of consistency over time or stationarity. One potential direction could be to integrate our methods in the strategy proposed by Saggioro et al. [2020] which combines a causal discovery with a regime learning optimisation approach. However, this direction might require additional assumptions and while it appears promising for CBNB, it is less obvious for NBCB. Another direction could be to assume the presence of an observed contextual variable that explains the non-stationarity [Mooij et al., 2020, Günther et al., 2023].

Finally, adapting these methods for mixed data can be important in many applications. This can be straightforward for the constraint-based part of our algorithms, given the availability of conditional independence tests for mixed data [Zan et al., 2022]. However, the adaptation is more challenging for the noise-based part.

# 7 Conclusion

In this paper, we introduced a framework for hybrids of noise-based and constraint-based methods that can discover causal graphs from temporal data. Algorithms in the first class, denoted NBCB, start with ordering instantaneous relations, and then prune edges of the fully oriented graph. On the other hand, the algorithms from the second class, denoted CBNB, start by finding the skeleton and the orientation of lagged relations using a constraint-based method and temporal priority, then orient instantaneous relations by ordering nodes in each cycle group of instantaneous relations. Overall, the performance of our algorithms is a trade-off between the performance of constraint-based and noise-based algorithms when all assumptions are satisfied, and they outperform other methods in the cases where some of the assumptions are violated.

For future works, it would be interesting to extend these approaches to cases involving hidden confounders, non-stationarity, and mixed data.

# Acknowledgments

We thank Ali Aït-Bachir, Christophe de Bignicourt and Rachid Mokhtari from EasyVista for providing the IT monitoring data along with the underlying causal graphs. We also thank Giovanni Poggiato for several discussions about the realistic ecological simulated data. Finally, we thank the anonymous reviewers for their helpful comments. This work was partially supported by the CIPHOD project (ANR-23-CPJ1-0212-01), by MIAI@Grenoble Alpes (ANR-19-P3IA-0003), by the Horizon Europe Obsession project (No: 101134954) and the FRB-CESAB through the IMPACT working group.

# References

* Hugo Arlegui, Emmanuel Mahé, M.-A Richard, Yann Rycke, M. Viguier, Marie Beylot-Barry, A. Dupuy, N. Beneton, P. Joly, Denis Jullien, H. Bachelez, É Sbidian, O. Chosidow, Carle Paul, and F. Tubach. Impact of the first wave of the covid-19 pandemic on the treatment of psoriasis with systemic therapies in france: Results from the psobioteq cohort. *Annales de Dermatologie et de Vénéréologie*, 150, 02 2023. doi: 10.1016/j.annder.2023.01.005.

* Andrew Arnold, Yan Liu, and Naoki Abe. Temporal causal modeling with graphical granger methods. In *Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, KDD '07, page 66–75, New York, NY, USA, 2007. Association for Computing Machinery. ISBN 9781595936097. doi: 10.1145/1281192.1281203.

* Charles K. Assaad, Emilie Devijver, Eric Gaussier, and Ali Aït-Bachir. A mixed noise and constraint-based approach to causal inference in time series. In *Machine Learning and Knowledge Discovery in Databases. Research Track*, pages 453–468, Cham, 2021. Springer International Publishing. ISBN 978-3-030-86486-6.

* Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. *J. Artif. Int. Res.*, 73, may 2022a. doi: 10.1613/jair.1.13428.

* Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Entropy-based discovery of summary causal graphs in time series. *Entropy*, 24(8), 2022b. ISSN 1099-4300. doi: 10.3390/e24081156.
Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Discovery of extended summary graphs in time series. In James Cussens and Kun Zhang, editors, *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, volume 180 of *Proceedings of Machine Learning Research*, pages 96–106. PMLR, 01–05 Aug 2022c.

Charles K. Assaad, Imad Ez-Zejjari, and Lei Zan. Root cause identification for collective anomalies in time series given an acyclic summary causal graph with loops. In Francisco Ruiz, Jennifer Dy, and Jan-Willem van de Meent, editors, *Proceedings of The 26th International Conference on Artificial Intelligence and Statistics*, volume 206 of *Proceedings of Machine Learning Research*, pages 8395–8404. PMLR, 25–27 Apr 2023.

Mélanie Bales, Elodie Pambrun, Charlotte Maguet, Judith van der Waerden, Nine Glangeaud-Freudenthal, Marie-Aline Charles, Corinne Bois, Maria Melchior, Jeannette Milgrom, Bruno Falissard, Helene Verdoux, and Anne Laure Sutter. Pathways between risk/protective factors and maternal postnatal depressive symptoms: The elfe cohort. *Journal of Clinical Medicine*, 12:3204, 04 2023. doi: 10.3390/jcm12093204.

Frédéric Barraquand, Coralie Picoche, Matteo Detto, and Florian Hartig. Inferring species interactions using granger causality and convergent cross mapping. *Theoretical Ecology*, 14(1):87–105, 2021.

Wei Chen, Kun Zhang, Ruichu Cai, Biwei Huang, Joseph Ramsey, Zhifeng Hao, and Clark Glymour. Fritl: A hybrid method for causal discovery in the presence of latent confounders. *arXiv preprint arXiv:2103.14238*, 2021.

David Maxwell Chickering. Learning equivalence classes of bayesian-network structures. *J. Mach. Learn. Res.*, 2: 445–498, mar 2002a. ISSN 1532-4435. doi: 10.1162/153244302760200696.

David Maxwell Chickering. Optimal structure identification with greedy search. *Journal of machine learning research*, 3(Nov):507–554, 2002b.

David Maxwell Chickering and Christopher Meek. Selective greedy equivalence search: Finding optimal bayesian networks using a polynomial number of score evaluations. *arXiv preprint arXiv:1506.02113*, 2015.

Maria Dornelas, Laura H Antao, Faye Moyes, Amanda E Bates, Anne E Magurran, Dušan Adam, Asem A Akhmet-zhanova, Ward Appeltans, Jose Manuel Arcos, Haley Arnold, et al. Biotime: A database of biodiversity time series for the anthropocene. *Global Ecology and Biogeography*, 27(7):760–786, 2018.

Doris Entner and Patrik O Hoyer. On causal discovery from time series data using fci. *Probabilistic graphical models*, pages 121–128, 2010.

Chang Gong, Di Yao, Chuzhe Zhang, Wenbin Li, Jingping Bi, Lun Du, and Jin Wang. Causal discovery from temporal data. In *Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, KDD '23, page 5803–5804, New York, NY, USA, 2023. Association for Computing Machinery. ISBN 9798400701030. doi: 10.1145/3580305.3599552.

Clive Granger. Investigating causal relations by econometric models and cross-spectral methods. *Econometrica*, 37 (3):424–38, 1969.

Clive W. J. Granger. Time series analysis, cointegration, and applications. *The American Economic Review*, 94(3): 421–425, 2004. ISSN 00028282.

Wiebke Günther, Urmi Ninad, and Jakob Runge. Causal discovery for time series from multiple datasets with latent contexts. In Robin J. Evans and Ilya Shpitser, editors, *Proceedings of the Thirty-Ninth Conference on Uncertainty in Artificial Intelligence*, volume 216 of *Proceedings of Machine Learning Research*, pages 766–776. PMLR, 31 Jul–04 Aug 2023.

Uzma Hasan, Emam Hossain, and Md Osman Gani. A survey on causal discovery methods for iid and time series data. *Transactions on Machine Learning Research*, 2023.

Patrik O. Hoyer, Aapo Hyvärinen, Richard Scheines, Peter L. Spirtes, Joseph Ramsey, Gustavo Lacerda, and Shohei Shimizu. Causal discovery of linear acyclic models with arbitrary distributions. In *Conference on Uncertainty in Artificial Intelligence*, 2008.

Aapo Hyvärinen, Shohei Shimizu, and Patrik O. Hoyer. Causal modelling combining instantaneous and lagged effects: An identifiable model based on non-gaussianity. In *Proceedings of the 25th International Conference on Machine Learning*, ICML '08, pages 424–431, New York, NY, USA, 2008. ACM. ISBN 978-1-60558-205-4.

Christian Jost and Stephen P Ellner. Testing for predator dependence in predator-prey dynamics: a non-parametric approach. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 267(1453):1611–1620, 2000.

Marcus Kaiser and Maksim Sipos. Unsuitability of notears for causal graph discovery when dealing with dimensional quantities. *Neural Processing Letters*, 54:1587–1595, 2021.

Yan Li, Rui Xia, Chunchen Liu, and Liang Sun. A hybrid causal structure learning algorithm for mixed-type data. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 36, pages 7435–7443, 2022.
Charles K. Assaad, Emilie Devijver, and Eric Gaussier. Discovery of extended summary graphs in time series. In James Cussens and Kun Zhang, editors, *Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence*, volume 180 of *Proceedings of Machine Learning Research*, pages 96–106. PMLR, 01–05 Aug 2022c.

Charles K. Assaad, Imad Ez-Zejjari, and Lei Zan. Root cause identification for collective anomalies in time series given an acyclic summary causal graph with loops. In Francisco Ruiz, Jennifer Dy, and Jan-Willem van de Meent, editors, *Proceedings of The 26th International Conference on Artificial Intelligence and Statistics*, volume 206 of *Proceedings of Machine Learning Research*, pages 8395–8404. PMLR, 25–27 Apr 2023.

Mélanie Bales, Elodie Pambrun, Charlotte Maguet, Judith van der Waerden, Nine Glangeaud-Freudenthal, Marie-Aline Charles, Corinne Bois, Maria Melchior, Jeannette Milgrom, Bruno Falissard, Helene Verdoux, and Anne Laure Sutter. Pathways between risk/protective factors and maternal postnatal depressive symptoms: The elfe cohort. *Journal of Clinical Medicine*, 12:3204, 04 2023. doi: 10.3390/jcm12093204.

Frédéric Barraquand, Coralie Picoche, Matteo Detto, and Florian Hartig. Inferring species interactions using granger causality and convergent cross mapping. *Theoretical Ecology*, 14(1):87–105, 2021.

Wei Chen, Kun Zhang, Ruichu Cai, Biwei Huang, Joseph Ramsey, Zhifeng Hao, and Clark Glymour. Fritl: A hybrid method for causal discovery in the presence of latent confounders. *arXiv preprint arXiv:2103.14238*, 2021.

David Maxwell Chickering. Learning equivalence classes of bayesian-network structures. *J. Mach. Learn. Res.*, 2: 445–498, mar 2002a. ISSN 1532-4435. doi: 10.1162/153244302760200696.

David Maxwell Chickering. Optimal structure identification with greedy search. *Journal of machine learning research*, 3(Nov):507–554, 2002b.

David Maxwell Chickering and Christopher Meek. Selective greedy equivalence search: Finding optimal bayesian networks using a polynomial number of score evaluations. *arXiv preprint arXiv:1506.02113*, 2015.

Maria Dornelas, Laura H Antao, Faye Moyes, Amanda E Bates, Anne E Magurran, Dušan Adam, Asem A Akhmet-zhanova, Ward Appeltans, Jose Manuel Arcos, Haley Arnold, et al. Biotime: A database of biodiversity time series for the anthropocene. *Global Ecology and Biogeography*, 27(7):760–786, 2018.

Doris Entner and Patrik O Hoyer. On causal discovery from time series data using fci. *Probabilistic graphical models*, pages 121–128, 2010.

Chang Gong, Di Yao, Chuzhe Zhang, Wenbin Li, Jingping Bi, Lun Du, and Jin Wang. Causal discovery from temporal data. In *Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, KDD '23, page 5803–5804, New York, NY, USA, 2023. Association for Computing Machinery. ISBN 9798400701030. doi: 10.1145/3580305.3599552.

Clive Granger. Investigating causal relations by econometric models and cross-spectral methods. *Econometrica*, 37 (3):424–38, 1969.

Clive W. J. Granger. Time series analysis, cointegration, and applications. *The American Economic Review*, 94(3): 421–425, 2004. ISSN 00028282.

Wiebke Günther, Urmi Ninad, and Jakob Runge. Causal discovery for time series from multiple datasets with latent contexts. In Robin J. Evans and Ilya Shpitser, editors, *Proceedings of the Thirty-Ninth Conference on Uncertainty in Artificial Intelligence*, volume 216 of *Proceedings of Machine Learning Research*, pages 766–776. PMLR, 31 Jul–04 Aug 2023.

Uzma Hasan, Emam Hossain, and Md Osman Gani. A survey on causal discovery methods for iid and time series data. *Transactions on Machine Learning Research*, 2023.

Patrik O. Hoyer, Aapo Hyvärinen, Richard Scheines, Peter L. Spirtes, Joseph Ramsey, Gustavo Lacerda, and Shohei Shimizu. Causal discovery of linear acyclic models with arbitrary distributions. In *Conference on Uncertainty in Artificial Intelligence*, 2008.

Aapo Hyvärinen, Shohei Shimizu, and Patrik O. Hoyer. Causal modelling combining instantaneous and lagged effects: An identifiable model based on non-gaussianity. In *Proceedings of the 25th International Conference on Machine Learning*, ICML '08, pages 424–431, New York, NY, USA, 2008. ACM. ISBN 978-1-60558-205-4.

Christian Jost and Stephen P Ellner. Testing for predator dependence in predator-prey dynamics: a non-parametric approach. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 267(1453):1611–1620, 2000.

Marcus Kaiser and Maksim Sipos. Unsuitability of notears for causal graph discovery when dealing with dimensional quantities. *Neural Processing Letters*, 54:1587–1595, 2021.

Yan Li, Rui Xia, Chunchen Liu, and Liang Sun. A hybrid causal structure learning algorithm for mixed-type data. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 36, pages 7435–7443, 2022.