equation(#1)

# Structural Hawkes Processes for Learning Causal Structure from Discrete-Time
Event Sequences

Jie Qiao1 Ruichu Cai1,2 Corresponding author. Siyu Wu1 Yu Xiang1 Keli
Zhang3&Zhifeng Hao4 1School of Computer Science, Guangdong University of
Technology, Guangzhou 510006, China
2Peng Cheng Laboratory, Shenzhen 518066, China
3Huawei Noah’s Ark Lab, Huawei, Shenzhen 518116, China
4College of Science, Shantou University, Shantou 515063, China
qiaojie.chn@gmail.com,cairuichu<EMAIL_ADDRESS>thexiang2000@gmail.com,zhangkeli1<EMAIL_ADDRESS>
###### Abstract

Learning causal structure among event types from discrete-time event sequences
is a particularly important but challenging task. Existing methods, such as
the multivariate Hawkes processes based methods, mostly boil down to learning
the so-called Granger causality which assumes that the cause event happens
strictly prior to its effect event. Such an assumption is often untenable
beyond applications, especially when dealing with discrete-time event
sequences in low-resolution; and typical discrete Hawkes processes mainly
suffer from identifiability issues raised by the instantaneous effect, i.e.,
the causal relationship that occurred simultaneously due to the low-resolution
data will not be captured by Granger causality. In this work, we propose
Structure Hawkes Processes (SHPs) that leverage the instantaneous effect for
learning the causal structure among events type in discrete-time event
sequence. The proposed method is featured with the minorization-maximization
of the likelihood function and a sparse optimization scheme. Theoretical
results show that the instantaneous effect is a blessing rather than a curse,
and the causal structure is identifiable under the existence of the
instantaneous effect. Experiments on synthetic and real-world data verify the
effectiveness of the proposed method.

##  1 Introduction

Learning causal structure among event types on multi-type event sequences is
an important and challenging task, and has recently found applications in
social science Zhou et al. (2013a), economic Bacry et al. (2015), network
operation maintenance Cai et al. (2022), etc. Existing methods, such as the
multivariate Hawkes processes based methods Xu et al. (2016); Bhattacharjya et
al. (2018); Salehi et al. (2019), mostly boil down to learning the so-called
Granger causality Granger (1969) which implicitly assumes that all events are
recorded instantaneously and accurately such that the cause event happens
strictly prior to its effect event (known as temporal precedence assumption).
However, due to the limited recording capabilities and storage capacities,
retaining event’s occurred times with high-resolution is expensive or
practically impossible in many real-world applications, and we usually only
can access the corresponding discrete-time event sequences. For example, in
large wireless networks, the event sequences are usually logged at a certain
frequency by different devices whose time might not be accurately
synchronized. As a result, low-resolution discrete-time event sequences are
obtained and the temporal precedence assumption will be frequently violated in
discrete-time event sequences, which raises a serious identifiability issue of
causal discovery. For example, as shown in Fig. 1, there are three event
sequences produced by three event types $v_{1}$, $v_{2}$, and $v_{3}$,
respectively. Let $v_{1}$ be the cause of $v_{2}$ and $v_{3}$, Fig. 1(a) shows
the accurate continuous-time event sequences such that each event occurred
time is recorded by $t_{1},...,t_{6}$. However, such high-resolution sequences
are usually not accessible and we can only observe the discrete-time event
sequences as shown in Fig. 1(b). As a result, $v_{1}$ and $v_{2}$ will be
considered simultaneous events due to the low-resolution, which violates the
temporal precedence assumption. Consequently, many existing point process
based methods will fail to capture the causal relationship $v_{1}\to v_{2}$ as
only the events that occur earlier are considered as causes. In contrast, the
causal relationship $v_{1}\to v_{3}$ can still be captured using the point
process based method because the cause event $v_{1}$ occurs before $v_{3}$.
However, one can imagine that as the resolution becomes lower the causal
relationship between $v_{1}$ and $v_{3}$ might no longer be identified as they
might become occurs at the same time. Thus, in this paper, we aim to answer
the following two questions: 1) How to design and learn a Hawkes process that
leverages the instantaneous effect in discrete time? 2) Can we identify the
causal relationship in event sequences under the existence of the
instantaneous effect?

These two questions have to do with the point processes and the causal
discovery, respectively. The former question is about the design of the
discrete-time Hawkes processes. The effect of discretization of Hawkes
processes has been widely discussed in many ways Foufoula-Georgiou and
Lettenmaier (1986); Kirchner (2016); Shlomovich et al. (2022); Jacod and
Todorov (2009). However, how to design a Hawkes process that leverages the
instantaneous effect in discrete Hawkes process is still unknown. The latter
question is related to the causal discovery Zhang et al. (2018); Qiao et al.
(2021); Cai et al. (2018a, b); Yu et al. (2019, 2021); Yang et al. (2021).
Some methods have been developed for continuous-value time series with
instantaneous effect. For example, the structural vector autoregression model
Swanson and Granger (1997) is an extension of the vector autoregression model
with instantaneous effect, and it has been shown that under the linear non-
Gaussian model, the causal relationship is identifiable under the existence of
the instantaneous effect Hyvärinen et al. (2010). However, we are not aware of
any method that leverages the instantaneous effect to identify the causal
relationship in event sequences with point processes. Thus, in this work, we
propose Structural Hawkes Processes (SHPs) that leverage the instantaneous
effect in discrete-time event sequences. We will proceed as follows. In
Section 2, we review the related work. In Section 3.1, we show how to design
and learn the structural Hawkes processes. In Section 4, we investigate the
identification of structural Hawkes processes and theoretically show that the
causal relationship among the instantaneous effect is indeed identifiable. In
Section 5, we perform extensive empirical evaluations over both synthetic and
real-world data.

(a) Continuous-time event sequences.

(b) Discrete-time event sequences.

Figure 1: Toy example of the continuous-time and discrete-time event
sequences. We consider three types of event $v_{1},v_{2},v_{3}$ and the causal
relationship satisfying $v_{2}\leftarrow v_{1}\rightarrow v_{3}$. While, in
continuous-time, each event’s occur times $(t_{1},...,t_{6})$ can be observed
accurately, the discrete-time event sequences have lower resolution, such that
types $v_{1},v_{2}$ will be considered as simultaneous and the Granger-based
methods will fail to capture the causal relationship $v_{1}\to v_{2}$.

##  2 Related Work

This work is closely related to two topics: point processes and learning
causal structure from event sequence.

#### Point processes.

There are mainly two types of point processes for modeling the event sequence.
Most existing methods focus on developing the continues-time Hawkes processes
Hawkes (1971) which assume that the event in the past can influence the event
in the future. Different variants have been developed with different types of
intensity in the Hawkes processes, e.g., the parametric functions Farajtabar
et al. (2014); Zhou et al. (2013a); Rasmussen (2013); Cai et al. (2022), the
non-parametric functions Lewis and Mohler (2011); Zhou et al. (2013b); Achab
et al. (2017), and the recent deep learning based functions Du et al. (2016);
Mei and Eisner (2017); Shang and Sun (2019). Another line of research focuses
on the discrete-time Hawkes processes which are more appealing for certain
applications. Seol (2015) study the limit theorem for the discrete Hawkes
processes and it has been extended to the discrete marked Hawkes processes by
Wang (2022). Shlomovich et al. (2022) further discusses the estimation method
for 1-dimensional discrete Hawkes processes. However, none of the above
methods leverage the instantaneous effect for learning causal structure in
discrete-time event sequences.

#### Learning causal structure from event sequences.

Most approaches in learning causal structure from event sequences are based on
Granger causality Granger (1969). The basic idea of Granger causality is to
constrain that the effect cannot precede the cause and the causal analysis can
be conducted using predictability. In particular, many methods have been
developed for learning Granger causality from the continuous-time event
sequences based on multivariate Hawkes processes. Xu et al. (2016) proposes a
nonparametric Hawkes processes model with group sparsity regularization, while
Zhou et al. (2013a) proposes to use a nuclear and $\ell_{1}$ norm and Idé et
al. (2021) consider an $\ell_{0}$ norm as the sparse regularization. In
addition, Achab et al. (2017) proposes to use cumulant for learning the causal
structure without optimization. Recently, some deep point process based
methods have been proposed, e.g., Zhang et al. (2020) introduce an attribution
method to uncover the Granger causality. However, the Granger causal relations
could be misled in low-resolution which is also pointed out by Spirtes and
Zhang (2016). One remedy is to extend Granger causality with instantaneous
effects. It has been found that in continuous-value time series, one is able
to incorporate the instantaneous effect under the linear relation with non-
Gaussian noise Hyvärinen et al. (2010). Runge (2020) further proposes to use a
constraint-based method for learning the causal structure under the
instantaneous effects. However, the extensions are only applicable in a
restricted case and it is still very challenging to learn causal structure
from event sequences with instantaneous effects.

##  3 Structural Hawkes Processes

We begin with a brief introduction to the general continuous-time multivariate
point processes. We then develop the structural Hawkes processes that leverage
the instantaneous effect in discrete-time.

###  3.1 Multivariate Point Processes

A multivariate point process is a random process that can be presented via a
$|\mathbf{V}|$-dimensional counting process
$\mathbf{N}=\\{N_{v}(t)|t\in[0,T],v\in\mathbf{V}\\}$ where $N_{v}(t)=N((0,t])$
measures the number of events that occur before time $t$ in the event type
$v$. Each counting process $N_{v}(t)$ can be characterized by the conditional
intensity function $\lambda_{v}(t)$ satisfying:

$\lambda_{v}(t)dt=\mathbb{E}[dN_{v}(t)|{\mathcal{F}_{t-}}],$

where $\lambda_{v}(t)$ characterizes the (conditional) expected number of
jumps per unit of time, $dN_{v}(t)=N_{v}(t+dt)-N_{v}(t)$ measures the
increment of the jump, $\mathcal{F}_{t-}=\bigcup_{0\leq
s<t,v\in\mathbf{V}}\mathcal{F}^{v}_{s}$ in which $\mathcal{F}^{v}_{s}$ is the
canonical filtration of sub-processes $N_{v}(t)$ up to time $s$.

In particular, Hawkes process is a counting process with an intensity function
that has the following form:

$\displaystyle\lambda_{v}(t)$
$\displaystyle=\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\int_{0}^{t}\phi_{v^{\prime},v}(t-t^{\prime})dN_{v^{\prime}}(t^{\prime}),$
(1)

where $\mu_{v}$ is the immigration intensity and
$\sum_{v^{\prime}\in\mathbf{V}}\int_{0}^{t}\phi_{v^{\prime},v}(t-t^{\prime})dN_{v^{\prime}}(t^{\prime})$
is the endogenous intensity aiming to model the influence from other event
types occurring near time $t$ Farajtabar et al. (2014).
$\phi_{v^{\prime},v}(t)$ is a reproduction function characterizing the time
variation of the causal influence from event $v^{\prime}$ to $v$. Intuitively,
one can imagine that events can be seen to arrive either via immigration
according to $\mu_{v}$ or via reproduction from past events according to the
endogenous intensity.

###  3.2 Design of Structural Hawkes Processes

To model the instantaneous effect, we begin with extending the continues-time
counting processes into the discrete-time in which the event sequences are
observed or collected in $\mathbf{T}=\\{\Delta,2\Delta,...,K\Delta\\}$, for
$K=\lfloor T/\Delta\rfloor$ where $\Delta>0$ is the length of time interval at
each observed time. Then the multivariate counting processes in discrete-time
can defined as
$\mathbf{N}^{(\Delta)}=\\{N_{v}^{(\Delta)}(k)|k\in\\{0,\dots,K\\},v\in\mathbf{V}\\}$,
where $N_{v}^{(\Delta)}(k)=N_{v}((0,k\Delta])$ measures the number of events
that occurs not later than $k\Delta$. We further let
$\mathbf{X}=\\{X_{v,t}|v\in\mathbf{V},t\in\\{0,\dots,K\\}]\\}$ denote the set
of observations at each time interval where $X_{v,t}\coloneq
N_{v}(t\Delta)-N_{v}((t-1)\Delta)$ is an analogy to $dN_{v}(t)$.

The discrete-time counting processes, however, as we discussed earlier, will
ignore the instantaneous effect which could lead to a misleading result. To
tackle this issue, we propose structural counting processes with a new type of
conditional intensity function that leverages the events that occur at the
same period of time.

######  Definition 1 (Structural counting processes).

A structural counting process is a multivariate counting process
$\mathbf{N}^{(\Delta)}$ in discrete-time with the conditional intensity of
$N^{(\Delta)}_{v}$ for each $v\in\mathbf{V}$ satisfying:

$\lambda_{v}(k\Delta)\Delta=\mathbb{E}[X_{v,k}|{\mathcal{F}_{(k-1)\Delta}\cup\mathcal{F}^{-v}_{k\Delta}}],$
(2)

where $\mathcal{F}_{(k-1)\Delta}=\bigcup_{0\leq s\leq
k-1,v\in\mathbf{V}}\mathcal{F}^{v}_{s\Delta}$ is the filtration with discrete-
time in the past and
$\mathcal{F}^{-v}_{k\Delta}\coloneq\left\\{\mathcal{F}^{v^{\prime}}_{k\Delta}|v^{\prime}\in\mathbf{V}\backslash
v\right\\}$ is the filtration that except for type-$v$ event.

Note that the exclusion in $\mathcal{F}^{-v}_{k\Delta}$ is necessary since it
makes no sense to use the current number of events to predict itself. Based on
the structural counting processes in Definition 1, the structural Hawkes
processes can be designed as follows:

######  Definition 2 (Structural Hawkes processes).

A structural Hawkes process is a structural counting process such that for
each $v\in\mathbf{V}$, the intensity of $N^{(\Delta)}_{v}$ can be written as:

$\lambda_{v}(k\Delta)=\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum_{i=1}^{k}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i},$
(3)

where $\phi_{v,v}(0)\equiv 0$ ensures the exclusive of type-$v$ event at time
$k\Delta$.

We can see that the intensity in Eq. 3 is not only influenced by the events
that occur in the past $(k-1)\Delta$ but also the events that occur at the
same period of time $k\Delta$.

###  3.3 Learning of Structural Hawkes Processes

The goal of this work is to identify a proper directed acyclic graph (DAG)
$\mathcal{G}(\mathbf{V},\mathbf{E})$ among event types $\mathbf{V}$ such that
for a type-$v$ event, the intensity $\lambda_{v}(t)$ only depends on the event
from its cause $\mathbf{Pa}_{v}$ where $\mathbf{Pa}_{v}$ is the parent set of
$v$ according to the edge set $\mathbf{E}$, i.e., $\\{v^{\prime}\to
v|v^{\prime}\in\mathbf{Pa}_{v}\\}\subseteq\mathbf{E}$. Note that a DAG
constraint for the instantaneous causal relation is necessary due to the lack
of temporal precedence but it can be easily relaxed by taking lagged and
instantaneous relations separately which we omit here for simplicity.

To learn causal structure, we parameterize the causal influence into the
reproduction function $\phi_{v^{\prime},v}(t)=\alpha_{v^{\prime},v}\kappa(t)$
where $\alpha_{v^{\prime},v}$ denotes the causal strength of causal relation
$v^{\prime}\to v$ and $\kappa(t)$ characterizes the time-decay of the causal
influence which is usually set as an exponential form $\kappa(t)=\exp(\beta
t)$ for $t\geq 0$ with the hyper-parameter $\beta$. That is, the causal
relationship can be encoded by the impact function such that for any pairs
$v^{\prime},v\in\mathbf{V}$ if $\phi_{v^{\prime},v}(t)=0$ we have
$v^{\prime}\to v\not\in\mathbf{E}$.

Thus, one of the challenges of learning causal structure is to constrain the
sparsity of the reproduction function. To this end, we devise the objective
function with sparsity constraint using the $\ell_{0}$ norm under the
likelihood framework.

Given a collection of discrete-time sequences $\mathbf{X}$, the log-likelihood
of parameters
$\Theta=\\{\mathbf{A}=[\alpha_{v^{\prime},v}]\in\mathbb{R}^{|\mathbf{V}|\times|\mathbf{V}|},\mathbf{\mu}=[\mu_{v}]\in\mathbb{R}^{|\mathbf{V}|}\\}$
of SHP can be expressed as follows:

$\displaystyle\mathcal{L}(\mathcal{G},\Theta;\mathbf{X})$ (4) $\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}[-\lambda_{v}(k\Delta)\Delta+X_{v,k}\log(\lambda_{v}(k\Delta))]$
$\displaystyle+\underbrace{\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}[-\log(X_{v,k}!)+X_{v,k}\log(\Delta)]}_{\text{Const.}}$
$\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v,i}\right)\Delta\right.$
$\displaystyle\left.+X_{v,k}\log\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v,i}\right)\right]+\text{Const.}$

Without further constraint, the log-likelihood function will tend to produce
excessive redundant causal edges. Thus, we further penalize the model with
$\ell_{0}$ norm to enforce the sparsity and we obtain the objective function
as follows:

$\mathcal{L}_{p}(\mathcal{G},\Theta;\mathbf{X})=\mathcal{L}(\mathcal{G},\Theta;\mathbf{X})+\alpha_{S}\|\mathbf{A}\|_{0},$
(5)

where $\alpha_{S}$ controls the strength of $\ell_{0}$ norm. Although there
are different forms of sparse regularization for the constraint, e.g., Xu et
al. (2016) proposed to use the $\ell_{1}$ and $\ell_{2}$ norm, the $\ell_{0}$
can provide better sparsity performance and it is easy to be extended, for
example, considering the directed acyclic constraint Tsamardinos et al.
(2006), as many causal structures in real-world are often directed acyclic.

###  3.4 A Minorization-Maximization-based Algorithm for Learning Causal
Structure

However, estimation of the parameters by maximizing the likelihood in Eq. 5
has two obstacles. First, the likelihood of the point processes model is known
as flat making the optimization unstable and computationally intensive Veen
and Schoenberg (2008). Second, learning DAGs from observational data is a
combinatorial problem and, without any assumption, it has been shown to be NP-
hard Chickering et al. (2004). To tackle the issues above, following Lewis and
Mohler (2011); Chickering (2002), we propose to use minorization-maximization
(MM) based optimization Hunter and Lange (2004) with a two-step causal
structure learning algorithm.

For the first step, the parameters are estimated given a fixed causal graph
using the MM algorithm which leverages the additive structure in Eq. 4 to
apply Jensen’s inequality similar to the EM algorithm and obtain the following
lower bound:

$\displaystyle Q(\Theta|\Theta^{(j)})$ (6) $\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}\right)\Delta\right.$
$\displaystyle+X_{v,k}\left({\displaystyle
q^{\mu}_{v,k}}\log\left(\frac{\mu_{v}}{{\displaystyle
q^{\mu}_{v,k}}}\right)\right.$
$\displaystyle\left.\left.+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}\log\left(\frac{\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}}\right)\right)\right]$

where $q^{\mu}_{v,k}=\frac{\mu^{(j)}_{v}}{\lambda^{(j)}_{v}(k\Delta)}$ and
$q^{\alpha}_{v,k}(v^{\prime},i)=\frac{\phi^{(j)}_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}}{\lambda^{(j)}_{v}(k\Delta)}$.
$\lambda^{(j)}_{v}(k\Delta)$ is the conditional intensity function with
parameters $\Theta^{(j)}$ in the $j$-th iteration. Then, by setting
$\frac{\partial Q(\Theta|\Theta^{(j)})}{\partial\mu_{v}}=0$, $\frac{\partial
Q(\Theta|\Theta^{(j)})}{\partial\alpha_{v^{\prime},v}}=0$, we obtain the
close-form iteration formulas:

$\displaystyle\mu^{(j+1)}_{v}$
$\displaystyle=\frac{\sum\nolimits^{K}_{k=1}X_{v,k}q^{\mu}_{v,k}}{K\Delta}$
(7) $\displaystyle\alpha^{(j+1)}_{v^{\prime},v}$
$\displaystyle=\begin{cases}\frac{\sum\nolimits^{K}_{k=1}\sum^{k}_{i=1}q^{\alpha}_{v,k}(v^{\prime},i)X_{v,k}}{\sum\nolimits^{K}_{k=1}\sum^{k}_{i=1}\kappa((k-i)\Delta)X_{v^{\prime},i}\Delta}&v^{\prime}\neq
v\\\
\frac{\sum\nolimits^{K}_{k=1}\sum\nolimits^{k-1}_{i=1}q^{\alpha}_{v,k}(v^{\prime},i)X_{v,k}}{\sum\nolimits^{K}_{k=1}\sum\nolimits^{k-1}_{i=1}\kappa((k-i)\Delta)X_{v^{\prime},i}\Delta}&v^{\prime}=v\end{cases}$

Algorithm 1 Learning causal structure using SHP

Input: Data set $\mathbf{X}$
Output: $G^{*},\Theta^{*}$

1: $G^{\prime}\leftarrow empty~{}graph$,
$\mathcal{L}^{*}_{p}\leftarrow-\infty$

2: while
$\mathcal{L}^{*}_{p}(G^{*},\Theta^{*};\mathbf{X})<\mathcal{L}^{\prime}_{p}(G^{\prime},\Theta^{\prime};\mathbf{X})$
do

3: $G^{*},\Theta^{*}\leftarrow G^{\prime},\Theta^{\prime}$ with largest
$\mathcal{L}^{\prime}_{p}(G^{\prime},\Theta^{\prime};\mathbf{X})$

4: for every $G^{\prime}\in\mathcal{V}(G^{*})$ do

5: Update $\Theta^{\prime}$ via iteration in Eq. 7

6: Record score
$\mathcal{L}^{\prime}_{p}(G^{\prime},\Theta^{\prime};\mathbf{X})$

7: end for

8: end while

9: return $G^{*},\Theta^{*}$

With the MM algorithm, we then search the causal structure by using a Hill-
Climbing-based algorithm as shown in Algorithm 1. It mainly consists of two
phases. First, we perform a structure searching scheme by taking one step
adding, deleting, and reversing the graph $G^{*}$ in the last iteration, i.e.,
in Line 4, $\mathcal{V}(G^{*})$ represents a collection of the one-step
modified graph of $\mathcal{V}(G^{*})$. Furthermore, the acyclic constraint is
implemented by eliminating all cyclic causal graphs in $\mathcal{V}(G^{*})$.
Second, by fixing the graph $G^{\prime}$, we optimize the log-likelihood using
the MM algorithm in Line 5. Iterating the two steps above until the likelihood
no longer increases.

##  4 Identifiability

In this section, we aim to answer the question that whether we can identify
the causal relationship under the existence of the instantaneous effect. In
answering this question, one will need to explore the property of Hawkes
processes in discrete-time. Based on the discrete-time likelihood in Eq. 4,
each interval is modeled by conditional Poisson distribution with the linear
structure in Eq. 3. As such, one may alternatively represent the relation by
integer-valued autoregressive (INAR) processes according to Kirchner (2016)
for analyzing the identification of structural Hawkes Processes. Furthermore,
we assume that the causal sufficiency holds, i.e., all relevant variables have
been observed Spirtes et al. (2000). Following Kirchner (2016), the
INAR($\infty$) processes can be defined as follows:

######  Definition 3 (INAR($\infty$)).

For $\theta_{k}\geq 0,k\in\mathbb{N}_{0}$, let
$\epsilon_{t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\operatorname{Pois}(\theta_{0}),t\in\mathbb{N}$, and
$\xi_{i}^{(t,k)}\sim\operatorname{Pois}(\theta_{k})$. An Integer-valued
autoregressive time series of infinite order (INAR($\infty$)) process
${X_{t}},t\in\mathbb{N}$ is defined by

$X_{t}=\sum_{k=1}^{\infty}\theta_{k}\circ X_{t-k}+\epsilon_{t}$ (8)

where $\circ$ is a reproduction operator given by $\theta_{k}\circ
X_{t-k}\equiv\sum_{i=1}^{X_{t-k}}\xi_{i}^{(t,k)}$ with $\xi_{i}^{(t,k)}$ be a
sequence of i.i.d. non-negative integer-valued random variables that depends
on the reproduction coefficients $\theta_{k}$,
$\\{\epsilon_{t}\\}_{t\in\mathbb{N}}$ is an i.i.d. integer-valued immigration
sequence that are independent of $\\{\xi_{i}^{(t,k)}\\}$, and $X_{t-k}$ is
independent of $\epsilon_{t}$ for all $k$.

In general, INAR is a discrete analogy of the continuous autoregressive model.
Note that the distribution of the independent variables in INAR could be
different and different choices of the distribution would lead to different
INAR models Guerrero et al. (2022). Here, we use the Poisson choice to adopt
the conditional Poisson distribution in Eq. 4, and it has been shown that such
INAR($\infty$) model will converge to the continuous-time Hawkes processes as
the time interval $\Delta\to 0$:

######  Theorem 1 (Kirchner (2016)).

Let $N$ be a Hawkes process with immigration intensity $\mu$ and let
$\phi:\mathbb{R}\to\mathbb{R}_{0}^{+}$ be a reproduction intensity that is
piecewise continuous with $\phi(t)=0,t\leq 0$ and $\int\phi(t)dt<1$. For
$\Delta\in(0,\delta)$, let $\left(X_{t}^{(\Delta)}\right)$ be an
INAR($\infty$) sequence with immigration parameter $\Delta\mu$ and
reproduction coefficients $\Delta\phi(k\Delta),k\in\mathbb{N}$. From the
sequences ${(X_{t}^{(\Delta)})}_{\Delta\in(0,\delta)}$, we define a family of
point processes by

$N^{(\Delta)}(A)\coloneq\sum_{k:k\Delta\in A}X_{k}^{(\Delta)},\quad
A\in\mathcal{B},\Delta\in(0,\delta),$ (9)

where $\mathcal{B}\coloneq\mathcal{B}(\mathbb{R})$ is the Borel set in
$\mathbb{R}$. Then, we have that

$N^{(\Delta)}\stackrel{{\scriptstyle\mathrm{w}}}{{\longrightarrow}}N\quad\text{
for }\Delta\rightarrow 0.$ (10)

Theorem 1 implies that the properties of INAR can be utilized for analyzing
the Hawkes processes in discrete-time. Intuitively, the reproduction function
$\phi$ in the discrete-time Hawkes processes can be represented by a series of
reproduction coefficients $\theta_{k}$ for each time period, and the
immigration intensity $\mu$ in the discrete-time Hawkes processes can be
represented by the immigration parameters of $\epsilon_{t}$.

Specifically, given the property of INAR processes, the analysis of the
identifiability of structural Hawkes processes can be typically performed in
two folds—the identifiability of the temporal structural Hawkes processes and
the instantaneous structural Hawkes processes Hyvärinen et al. (2010). For the
former, the temporal resolution of the measurement that is high enough to
capture the former and latter relationship of the events, the identifiability
can be derived by local independence, which has been well explored by Mogensen
and Hansen (2020). For example, in Fig. 1, the causal relationship $v_{1}\to
v_{3}$ can be simply identified by the independence. Thus, in this work, we
are more interested in the instantaneous structural Hawkes processes—the
measurements have lower time resolution such that the causal influences are
instantaneous. In such a case, one can use a model in which the influences are
instantaneous, leading to Bayesian networks (BNs) or structural equation
models (SEMs), e.g., the causal relationship $v_{1}\to v_{2}$ in Fig. 1 belong
to this class.

Thus, to analyze the identification of SHP, based on the INAR model, we
consider the instantaneous causal structure in the structural Hawkes process:

######  Definition 4 (Instantaneous causal structure in structural Hawkes
process).

Let $\epsilon_{v,t}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\operatorname{Pois}(\mu_{v})$, and
$\xi_{i}^{(v^{\prime},v)}\sim\operatorname{Pois}(\alpha_{v^{\prime},v})$. The
instantaneous causal structure in the structural Hawkes process consists of a
set of equations of the form

$X_{v,t}=\sum_{v^{\prime}\in\mathbf{V}}\alpha_{v^{\prime},v}\circ
X_{v^{\prime},t}+\epsilon_{v},\quad v\in\mathbf{V},$ (11)

where $\alpha_{v^{\prime},v}>0$ for $v^{\prime}\in\mathbf{Pa}_{v}$ and
$\alpha_{v^{\prime},v}=0$ for $v^{\prime}\not\in\mathbf{Pa}_{v}$, with
$\alpha_{v^{\prime},v}\circ
X_{v^{\prime},t}\equiv\sum_{i=1}^{X_{v^{\prime},t}}\xi_{i}^{(v,v^{\prime})}$,
and $\mu_{v}>0$ if $v$ is a root variable i.e., the variable whose parent set
is empty $\mathbf{Pa}_{v}=\phi$. The random variables $X_{v^{\prime},t}$,
$\xi_{i}^{(v,v^{\prime})}$, and $\epsilon_{v}$ are independent of each other.

The identifiability of the instantaneous causal structure, however, is much
more difficult compared with the temporal causal relationship and most of them
suffer from lack of identifiability. For example, without any constraint, the
instantaneous causal structure may not be identified due to the Markov
equivalence class Pearl (2009) in which all graphs in the equivalence class
encode the same conditional independence. In such a case, different model with
different causal graph and parameter, $M_{1}=(\mathcal{G}_{1},\Theta_{1})$,
$M_{2}=(\mathcal{G}_{2},\Theta_{2})$ will produce the same distribution. Such
non-identifiability also exists even considering an additional constraint,
e.g., the linear Gaussian SEM is also non-identifiable Spirtes et al. (2000).

Thus, an essential goal of this work is to investigate the identifiability of
the instantaneous causal structure, i.e., whether there exists another causal
structure that entails the same distribution to the underlying causal model.
Such identifiability basically boils down to bivariate cases that whether
there exists a backward model that is distribution equivalent to the causal
direction. It can be easily extended to multivariate cases by incorporating
the conditional independent constraint in the causal structure with the causal
faithfulness assumption, i.e., the causal graph faithfully displays every
dependency Pearl (1988). Specifically, let $X$ be the cause variable and $Y$
be the effect variable, and surprisingly, the following theorem shows that the
bivariate instantaneous causal pair is indeed identifiable:

######  Theorem 2.

Let $X\to Y$ be the correct causal direction that follows

$\displaystyle Y$ $\displaystyle=\sum\nolimits_{i=1}^{X}\xi_{i}+\epsilon,$
$\displaystyle\quad X,\xi_{i},\textrm{ and }\epsilon\textrm{~{}are
independent},$ (12)

where $\xi_{i}\sim\operatorname{Pois}(\alpha_{X,Y})$,
$\epsilon\sim\operatorname{Pois}(\mu_{Y})$,
$X\sim\operatorname{Pois}(\mu_{X})$. Then, there does not exist a backward
model that admits the following equation:

$\displaystyle X$
$\displaystyle=\sum\nolimits_{i=1}^{Y}\hat{\xi}_{i}+\hat{\epsilon},$
$\displaystyle\quad Y,\hat{\xi}_{i},\textrm{ and }\hat{\epsilon}\textrm{~{}are
independent},$ (13)

where $\hat{\xi}_{i}\sim\operatorname{Pois}(\hat{\alpha}_{Y,X})$,
$\hat{\epsilon}\sim\operatorname{Pois}(\hat{\mu}_{X})$,
$Y\sim\operatorname{Pois}(\hat{\mu}_{Y})$.

To better understand the identifiability, we provide two alternative proofs
for this theorem. The details of the proofs are given in the supplementary
material. The first proof show that it is impossible to have a backward model
that the probability distribution of the two directions is equivalence. While
the second proof shows that the distribution of $Y$ must not admit the Poisson
distribution and therefore the distribution is not equivalent. We also provide
an empirical study of the identifiability for the bivariate causal pair in
Section 5. Finally, we generalize the result of Theorem 2 to the multivariate
causal structure and show that the multivariate instantaneous causal structure
is also identifiable:

######  Theorem 3.

With the causal faithfulness assumption and causal sufficiency assumption, the
multivariate instantaneous causal structure is identifiable.

The main idea behind the proof is that for any causal structure
$(\mathcal{G},\Theta)$. there does not exist another causal structure
$(\hat{\mathcal{G}},\hat{\Theta})$ that is distribution equivalent.
Specifically, with the causal faithfulness and causal sufficiency assumption,
the Markov equivalent class is identified. Since all Markov equivalent classes
share the same skeleton, we only need to show that each causal pair will not
admit another causal structure that has the revered causal direction while
having the same likelihood (i.e., distribution equivalent).

##  5 Experiments

In this section, we test the proposed SHP and the baselines on both synthetic
and real-world data. The baseline methods include ADM4 Zhou et al. (2013a),
NPHC Achab et al. (2017), MLE_SGL Xu et al. (2016), and PCMCI Plus Runge et
al. (2019). We further develop SHP_NH, an ablation study of SHP that removes
the hill-climb based searching scheme and uses a threshold to determine the
causal structure following the work of ADM4 and MLE_SGL. In all following
experiments, Recall, Precision, and F1 are used as the evaluation metrics. The
results of recall and precision are provided in the supplementary material.

###  5.1 Synthetic Experiments

Figure 2: The likelihood difference between two causal pairs in Theorem 2 with
different causal influence and immigration intensity.

(a) Sensitivity to Time Interval

(b) Sensitivity to Range of $\alpha$

(c) Sensitivity to Range of $\mu$

(d) Sensitivity to Sample Size

(e) Sensitivity to Num. of Event Types

(f) Sensitivity to Avg. Indegree

Figure 3: F1 in the Sensitivity Experiments

In this part, we design extensive control experiments using synthetic data to
test the correctness of the theory and the sensitivity of sample size, length
of time interval, number of event types, and different ranges of $\alpha$ and
$\mu$. In the sensitivity experiment, we synthesize data with fixed parameters
while traversing the target parameter as shown in Fig. 3. The default settings
are listed below, sample size=20000, time interval=5, number of event
types=20, range of $\alpha\in[0.3,0.5]$, range of $\mu\in[0.0005,0.0001]$. All
experimental results are averaged over 100 randomly generated causal
structures.

The generating process proceeds as follows: 1) randomly generate a directed
causal graph $\mathcal{G}$ with a certain average in degree; 2) generate the
events according to randomly generated parameters
$\alpha_{v^{\prime},v},\mu_{v}$ from $\mathcal{G}$, and 3) aggregate the
number of counts at each interval for each event according to the length of
time interval to synthesize the discrete-time process.

#### Two-variable case.

We first conduct a simple two-variable experiment to verify the
identifiability in Theorem 2, by computing the likelihood difference between
the causal direction and the reversed direction on a simulated causal pair
with different causal influence $\alpha_{X,Y}$. That is, we simulated data
using the model $Y_{t}=\alpha_{X,Y}\circ X_{t}+\epsilon_{t}$ with
$X_{t}\sim\operatorname{Pois}(\mu_{x})$ and
$\epsilon_{t}\sim\operatorname{Pois}(\mu_{y})$. Each experiment is conducted
100 times with random causal pairs. As shown in Fig. 2, we can see that the
likelihood difference is always greater than zero which means that it is
always identifiable unless the degenerate cases and it verifies the
correctness of Theorem 2. In addition, as the causal influence $\alpha_{X,Y}$
decreases, the likelihood difference also decreases and tends to zero. This is
reasonable as the causal influence becomes zero, and the two variables will
also become independent making the model non-identifiable. Moreover, the level
of immigration intensity also would affect the likelihood difference, which is
reasonable as the lower the $\mu_{x}$, the weaker the causal relation.
Similarly, the higher the $\mu_{y}$, the stronger the noise making the causal
relation weaker. The more general multivariate case experiments will be
implied in the following sensitivity analysis.

#### Sensitivity analysis.

As shown in Fig. 3, we conduct six different control experiments for SHP. In
general, our proposed SHP method outperforms all the baseline methods in all
six control experiments.

In the control experiments of time interval given in Fig. 3(a), as the time
interval controls the temporal resolution of the measurement sequence, the
larger the time interval, the lower the temporal resolution, and the more
instantaneous causal influences will occur leading to the decrease or
insensitive of performance of the baseline methods. In the contract, SHP keeps
giving the best results at all intervals. We also notice that the performance
will decrease as the time interval increase because it reduces the sample size
and the sequence become less informative. In addition, the performance of
SHP_NH stresses the effectiveness of the searching algorithm which has lower
precision and stability than SHP but also outperforms the baseline methods if
the event sequences contain sufficient information on the instantaneous
effect, which also verifies the ability to capture the instantaneous effect in
SHP.

For the controlled experiments of causal strength and immigration intensity in
Fig. 3(b) and Fig. 3(c), we can see that a reasonable causal strength and
immigration intensity are required for the baseline methods while SHP is
insensitive to both of them which shows the robustness of SHP.

In the sample size controlled experiments given in Fig. 3(d), all methods are
robust to the sample size, and in particular, a 1000 sample size is enough to
produce a reasonable performance for SHP, which demonstrates the practicality
of SHP.

For the causal structure controlled experiments given in Fig. 3(e) and Fig.
3(f), SHP performs well in both experiments, and the average indegree is also
important as it increases the performance of most methods decreases but SHP
has the smallest decline.

###  5.2 Real World Experiments

Figure 4: Real-world experiment on different temporal resolutions.

We also test the proposed SHP on a very challenging real-world
dataset111https://competition.huaweicloud.com/informations/mobile/1000041487/dataset
from real telecommunication networks. The dataset records eight months of
alarms that occurred in a real metropolitan cellular network. The alarms are
generated by fifty-five devices in the metropolitan cellular network, which
consists of eighteen types of alarms. Our goal is to learn the causal
structure among the alarm types. Note that the causal impact among alarms is
not deterministic but probabilistic in this system, which is similar to our
model setting. All experiments have been conducted with different random seeds
and the results are significant according to Wilcoxon signed-rank test. In
addition, the ground truths of the causal relationships among alarm types are
provided by domain experts.

As shown in Fig. 4, we conduct the real-world experiments by manually setting
different temporal resolutions of the observed sequence (i.e., from one second
to nine seconds) and the original temporal resolution is one second.
Interestingly, as the time interval increases, the performance of our method
also increases instead of decreasing, while the baseline methods all decrease.
The reason is that in the real-world scenario, there exists a communication
latency in the logging system such that the recorded timestamp might be not
fully accurate, and this error can be mitigated by decreasing the temporal
resolution. In the contrast, after decreasing the temporal resolution, other
methods fail to capture the causal relationship under the instantaneous
effect. This verifies the effectiveness of SHP and stresses the importance of
the instantaneous effect.

##  6 Conclusion

In this work, we study how to model and leverage the instantaneous effect for
learning causal structure in discrete-time event sequences. We propose
structural Hawkes processes that leverage the instantaneous effects and a
practical algorithm for learning the causal structure among event types.
Theoretical results show that the instantaneous causal structure in structural
Hawkes processes is indeed identifiable. To the best of our knowledge, this is
the first causal structure learning method for event sequences with
instantaneous effects. The success of SHP not only provides an effective
solution for learning causal structure from real-world event sequences but
also shows a promising direction for causal discovery from the discrete-time
event sequences. In the future, we plan to extend our work to a general point
process with a more general intensity function.

## Acknowledgements

This research was supported in part by National Key R&D Program of China
(2021ZD0111501), National Science Fund for Excellent Young Scholars
(62122022), Natural Science Foundation of China (61876043, 61976052), the
major key project of PCL (PCL2021A12). We appreciate the comments from
anonymous reviewers, which greatly helped to improve the paper.

## References

  * Achab et al. [2017] Massil Achab, Emmanuel Bacry, Stéphane Gaïffas, Iacopo Mastromatteo, and Jean-François Muzy.  Uncovering causality from multivariate Hawkes integrated cumulants.  Journal of Machine Learning Research, 18(1):6998–7025, 2017.
  * Bacry et al. [2015] Emmanuel Bacry, Iacopo Mastromatteo, and Jean-François Muzy.  Hawkes processes in finance.  Market Microstructure and Liquidity, 1(01):1550005, 2015.
  * Bacry et al. [2018] Emmanuel Bacry, Martin Bompaire, Philip Deegan, Stéphane Gaïffas, and Søren V. Poulsen.  tick: a python library for statistical learning, with an emphasis on hawkes processes and time-dependent models.  Journal of Machine Learning Research, 18(214):1–5, 2018.
  * Bhattacharjya et al. [2018] Debarun Bhattacharjya, Dharmashankar Subramanian, and Tian Gao.  Proximal graphical event models.  In Advances in Neural Information Processing Systems, pages 8136–8145, 2018.
  * Cai et al. [2018a] Ruichu Cai, Jie Qiao, Kun Zhang, Zhenjie Zhang, and Zhifeng Hao.  Causal discovery from discrete data using hidden compact representation.  In Advances in Neural Information Processing Systems, pages 2666–2674, 2018.
  * Cai et al. [2018b] Ruichu Cai, Jie Qiao, Zhenjie Zhang, and Zhifeng Hao.  SELF: Structural equational embedded likelihood framework for causal discovery.  In AAAI, 2018.
  * Cai et al. [2022] Ruichu Cai, Siyu Wu, Jie Qiao, Zhifeng Hao, Keli Zhang, and Xi Zhang.  THPs: Topological Hawkes Processes for Learning Causal Structure on Event Sequences.  IEEE Transactions on Neural Networks and Learning Systems, 2022\.
  * Chickering et al. [2004] Max Chickering, David Heckerman, and Chris Meek.  Large-sample learning of Bayesian networks is NP-hard.  Journal of Machine Learning Research, 5:1287–1330, 2004.
  * Chickering [2002] David Maxwell Chickering.  Optimal structure identification with greedy search.  Journal of machine learning research, 3(Nov):507–554, 2002.
  * Du et al. [2016] Nan Du, Hanjun Dai, Rakshit Trivedi, Utkarsh Upadhyay, Manuel Gomez-Rodriguez, and Le Song.  Recurrent marked temporal point processes: Embedding event history to vector.  In International Conference on Knowledge Discovery and Data Mining, pages 1555–1564, 2016.
  * Farajtabar et al. [2014] Mehrdad Farajtabar, Nan Du, Manuel Gomez Rodriguez, Isabel Valera, Hongyuan Zha, and Le Song.  Shaping social activity by incentivizing users.  In Advances in Neural Information Processing Systems, pages 2474–2482, 2014.
  * Foufoula-Georgiou and Lettenmaier [1986] Efi Foufoula-Georgiou and Dennis P Lettenmaier.  Continuous-time versus discrete-time point process models for rainfall occurrence series.  Water Resources Research, 22(4):531–542, 1986.
  * Granger [1969] Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica: Journal of the Econometric Society, pages 424–438, 1969.
  * Guerrero et al. [2022] Matheus B Guerrero, Wagner Barreto-Souza, and Hernando Ombao.  Integer-valued autoregressive processes with prespecified marginal and innovation distributions: A novel perspective.  Stochastic Models, 38(1):70–90, 2022.
  * Hawkes [1971] Alan G Hawkes.  Spectra of some self-exciting and mutually exciting point processes.  Biometrika, 58(1):83–90, 1971.
  * Hunter and Lange [2004] David R Hunter and Kenneth Lange.  A tutorial on MM algorithms.  The American Statistician, 58(1):30–37, 2004.
  * Hyvärinen et al. [2010] Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer.  Estimation of a structural vector autoregression model using non-gaussianity.  Journal of Machine Learning Research, 11(5), 2010.
  * Idé et al. [2021] Tsuyoshi Idé, Georgios Kollias, Dzung Phan, and Naoki Abe.  Cardinality-Regularized Hawkes-Granger Model.  Advances in Neural Information Processing Systems, 34:2682–2694, 2021.
  * Jacod and Todorov [2009] Jean Jacod and Viktor Todorov.  Testing for common arrivals of jumps for discretely observed multidimensional processes.  The Annals of Statistics, 37(4):1792–1838, 2009.
  * Kirchner [2016] Matthias Kirchner.  Hawkes and INAR ($\infty$) processes.  Stochastic Processes and their Applications, 126(8):2494–2525, 2016\.
  * Lewis and Mohler [2011] Erik Lewis and George Mohler.  A nonparametric EM algorithm for multiscale Hawkes processes.  Journal of Nonparametric Statistics, 1(1):1–20, 2011.
  * Mei and Eisner [2017] Hongyuan Mei and Jason M Eisner.  The neural hawkes process: A neurally self-modulating multivariate point process.  In Advances in Neural Information Processing Systems, pages 6754–6764, 2017.
  * Mogensen and Hansen [2020] Søren Wengel Mogensen and Niels Richard Hansen.  Markov equivalence of marginalized local independence graphs.  The Annals of Statistics, 48(1):539–559, 2020.
  * Pearl [1988] Judea Pearl.  Probabilistic reasoning in intelligent systems: networks of plausible inference.  Morgan kaufmann, 1988.
  * Pearl [2009] Judea Pearl.  Causality: Models, Reasoning and Inference.  Cambridge university press, 2009.
  * Qiao et al. [2021] Jie Qiao, Ruichu Cai, Kun Zhang, Zhenjie Zhang, and Zhifeng Hao.  Causal discovery with confounding cascade nonlinear additive noise models.  ACM Transactions on Intelligent Systems and Technology (TIST), 12(6):1–28, 2021.
  * Rasmussen [2013] Jakob Gulddahl Rasmussen.  Bayesian inference for Hawkes processes.  Methodology and Computing in Applied Probability, 15(3):623–642, 2013.
  * Runge et al. [2019] Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic.  Detecting and quantifying causal associations in large nonlinear time series datasets.  Science Advances, 5(11):eaau4996, 2019.
  * Runge [2020] Jakob Runge.  Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets.  In Uncertainty in Artificial Intelligence, pages 1388–1397. PMLR, 2020.
  * Salehi et al. [2019] Farnood Salehi, William Trouleau, Matthias Grossglauser, and Patrick Thiran.  Learning hawkes processes from a handful of events.  In Advances in Neural Information Processing Systems, pages 12694–12704, 2019.
  * Seol [2015] Youngsoo Seol.  Limit theorems for discrete Hawkes processes.  Statistics & Probability Letters, 99:223–229, 2015.
  * Shang and Sun [2019] Jin Shang and Mingxuan Sun.  Geometric hawkes processes with graph convolutional recurrent neural networks.  In AAAI Conference on Artificial Intelligence, volume 33, pages 4878–4885, 2019.
  * Shlomovich et al. [2022] Leigh Shlomovich, Edward AK Cohen, Niall Adams, and Lekha Patel.  Parameter estimation of binned hawkes processes.  Journal of Computational and Graphical Statistics, pages 1–11, 2022\.
  * Spirtes and Zhang [2016] Peter Spirtes and Kun Zhang.  Causal discovery and inference: concepts and recent methodological advances.  In Applied informatics, volume 3, pages 1–28. SpringerOpen, 2016\.
  * Spirtes et al. [2000] Peter Spirtes, Clark N Glymour, and Richard Scheines.  Causation, Prediction, and Search.  MIT press, 2000.
  * Swanson and Granger [1997] Norman R Swanson and Clive WJ Granger.  Impulse response functions based on a causal approach to residual orthogonalization in vector autoregressions.  Journal of the American Statistical Association, 92(437):357–367, 1997.
  * Tsamardinos et al. [2006] Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis.  The max-min hill-climbing Bayesian network structure learning algorithm.  Machine learning, 65(1):31–78, 2006.
  * Veen and Schoenberg [2008] Alejandro Veen and Frederic P Schoenberg.  Estimation of space–time branching process models in seismology using an em–type algorithm.  Journal of the American Statistical Association, 103(482):614–624, 2008.
  * Wang [2022] Haixu Wang.  Limit theorems for a discrete-time marked Hawkes process.  Statistics & Probability Letters, 184:109368, 2022.
  * Xu et al. [2016] Hongteng Xu, Mehrdad Farajtabar, and Hongyuan Zha.  Learning granger causality for hawkes processes.  In International Conference on Machine Learning, pages 1717–1726, 2016.
  * Yang et al. [2021] Shuai Yang, Kui Yu, Fuyuan Cao, Lin Liu, Hao Wang, and Jiuyong Li.  Learning causal representations for robust domain adaptation.  IEEE Transactions on Knowledge and Data Engineering, 2021.
  * Yu et al. [2019] Kui Yu, Lin Liu, Jiuyong Li, Wei Ding, and Thuc Duy Le.  Multi-source causal feature selection.  IEEE transactions on pattern analysis and machine intelligence, 42(9):2240–2256, 2019.
  * Yu et al. [2021] Kui Yu, Lin Liu, and Jiuyong Li.  A unified view of causal and non-causal feature selection.  ACM Transactions on Knowledge Discovery from Data (TKDD), 15(4):1–46, 2021.
  * Zhang et al. [2018] Kun Zhang, Bernhard Schölkopf, Peter Spirtes, and Clark Glymour.  Learning causality and causality-related learning: some recent progress.  National science review, 5(1):26–29, 2018.
  * Zhang et al. [2020] Wei Zhang, Thomas K. Panum, Somesh Jha, Prasad Chalasani, and David Page.  CAUSE: Learning granger causality from event sequences using attribution methods.  In International Conference on Machine Learning, 2020.
  * Zhou et al. [2013a] Ke Zhou, Hongyuan Zha, and Le Song.  Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes.  In Artificial Intelligence and Statistics, pages 641–649, 2013\.
  * Zhou et al. [2013b] Ke Zhou, Hongyuan Zha, and Le Song.  Learning triggering kernels for multi-dimensional hawkes processes.  In International Conference on Machine Learning, pages 1301–1309, 2013.

## Supplementary Material

In this supplementary material, we provide the derivation of minorization-
maximization algorithm in Section A, the proofs of Theorem 2 in Section B, the
proof of Theorem 3 in Section C, and the additional experiments results in
Section D.

##  Appendix A Derivation of minorization-maximization algorithm

In this section, we provide the detailed derivation of the proposed MM
algorithm. Given a set of observations $\mathbf{X}$, the log-likelihood of the
causal graph $\mathcal{G}$ and parameters $\Theta$, the log-likelihood can be
derived as follows:

$\displaystyle\mathcal{L}(\mathcal{G},\Theta;\mathbf{X})$ (A.1)
$\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\log\left[\frac{e^{-\lambda_{v}(k\Delta)\Delta}}{X_{v,k}!}(\lambda_{v}(k\Delta))\Delta)^{X_{v,k}}\right]$
$\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}[-\lambda_{v}(k\Delta)\Delta+X_{v,k}\log(\lambda_{v}(k\Delta))]$
$\displaystyle+\underbrace{\sum_{v\in\mathbf{V}}\sum\nolimits^{K}_{k=1}[-\log(X_{v,k}!)+X_{v,k}\log(\Delta)]}_{\coloneq\text{Const}}$
$\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}\right)\Delta\right.$
$\displaystyle\left.+X_{v,k}\log\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}\right)\right]+\text{Const}$
$\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}+\sum_{v^{\prime}\in\mathbf{V}\backslash
v}\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}\right)\Delta\right.$
$\displaystyle\left.+X_{v,k}\log\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}+\sum_{v^{\prime}\in\mathbf{V}\backslash
v}\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}\right)\right]+\text{Const},$

where $\displaystyle\phi_{v^{\prime},v}(t)=\begin{cases}0&v^{\prime}=v\
\text{and}\ t=0\\\
\alpha_{v^{\prime},v}\kappa(t)&\text{otherwise}\end{cases}$. By applying the
Jensen inequality to $\log(\lambda_{v}(t))$, we obtain the lower bound of the
intensity function:

$\displaystyle\log(\lambda_{v}(k\Delta))$ (A.2) $\displaystyle=$
$\displaystyle\log\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}+\sum_{v^{\prime}\in\mathbf{V}\backslash
v}\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}\right)$ $\displaystyle=$
$\displaystyle\log\left({\displaystyle
q^{\mu}_{v,k}\frac{\mu_{v}}{{\displaystyle
q^{\mu}_{v,k}}}}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)\frac{\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}}}+{\displaystyle\sum_{v^{\prime}\in\mathbf{V}\backslash
v}q^{\alpha}_{v,k}(v^{\prime},k)\frac{\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},k)}}}\right)$ $\displaystyle\geqslant$
$\displaystyle{\displaystyle
q^{\mu}_{v,k}}\log\left(\frac{\mu_{v}}{{\displaystyle
q^{\mu}_{v,k}}}\right)+{\displaystyle\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}q^{\alpha}_{v,k}(v^{\prime},i)}\log\left({\frac{\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}}}\right)+{\displaystyle\sum_{v^{\prime}\in\mathbf{V}\backslash
v}q^{\alpha}_{v,k}(v^{\prime},k)}\log\left({\frac{\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},k)}}}\right),$

where ${\displaystyle
q^{\mu}_{v,k}=\frac{\mu^{(j)}_{v}}{\lambda^{(j)}_{v}(k\Delta)}}$ and
${\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)=\frac{\phi^{(j)}_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}}{\lambda^{(j)}_{v}(k\Delta)}}$,
and $\lambda^{(j)}_{v}(k\Delta)$ is the conditional intensity function in the
$j$-th iteration. By subsisting Eq. A.2 into A.1, we obtain the objective
function in our work:

$\displaystyle Q(\Theta|\Theta^{(j)})$ (A.3) $\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}\right)\Delta\right.$
$\displaystyle+$ $\displaystyle\left.X_{v,k}\left({\displaystyle
q^{\mu}_{v,k}}\log\left(\frac{\mu_{v}}{{\displaystyle
q^{\mu}_{v,k}}}\right)+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k}_{i=1}{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}\log\left(\frac{\phi_{v^{\prime},v}((k-i)\Delta)X_{v^{\prime},i}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}}\right)\right)\right]$ $\displaystyle=$
$\displaystyle\sum\limits_{v\in\mathbf{V}}\sum\limits^{K}_{k=1}\left[-\left(\mu_{v}+\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}+\sum_{v^{\prime}\in\mathbf{V}\backslash
v}\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}\right)\Delta\right.$
$\displaystyle+$ $\displaystyle\left.X_{v,k}{\displaystyle
q^{\mu}_{v,k}}\log\left(\frac{\mu_{v}}{{\displaystyle
q^{\mu}_{v,k}}}\right)+{\displaystyle
X_{v,k}}\sum_{v^{\prime}\in\mathbf{V}}\sum^{k-1}_{i=1}{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}\log\left(\frac{\alpha_{v^{\prime},v}\kappa((k-i)\Delta)X_{v^{\prime},i}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},i)}}\right)+{\displaystyle
X_{v,k}}\sum_{v^{\prime}\in\mathbf{V}\backslash v}{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},k)}\log\left(\frac{\alpha_{v^{\prime},v}\kappa(0)X_{v^{\prime},k}}{{\displaystyle
q^{\alpha}_{v,k}(v^{\prime},k)}}\right)\right].$

Then, by setting $\frac{\partial Q(\Theta|\Theta^{(j)})}{\partial\mu_{v}}=0$
and $\frac{\partial Q(\Theta|\Theta^{(j)})}{\partial\alpha_{v^{\prime},v}}=0$,
we obtain the close-form iteration formulas:

$\displaystyle\mu^{(j+1)}_{v}$
$\displaystyle=\frac{\sum\nolimits^{K}_{k=1}X_{v,k}q^{\mu}_{v,k}}{K\Delta}$
(A.4) $\displaystyle\alpha^{(j+1)}_{v^{\prime},v}$
$\displaystyle=\begin{cases}\frac{\sum\nolimits^{K}_{k=1}\sum^{k}_{i=1}q^{\alpha}_{v,k}(v^{\prime},i)X_{v,k}}{\sum\nolimits^{K}_{k=1}\sum^{k}_{i=1}\kappa((k-i)\Delta)X_{v^{\prime},i}\Delta}&v^{\prime}\neq
v\\\
\frac{\sum\nolimits^{K}_{k=1}\sum\nolimits^{k-1}_{i=1}q^{\alpha}_{v,k}(v^{\prime},i)X_{v,k}}{\sum\nolimits^{K}_{k=1}\sum\nolimits^{k-1}_{i=1}\kappa((k-i)\Delta)X_{v^{\prime},i}\Delta}&v^{\prime}=v\end{cases}$

##  Appendix B Proof of Theorem 2

Here, to present a better understanding of the identifiability of the
instantaneous causal structure, we provide two different proofs for Theorem 2
from the perspective of the likelihood and the distribution of $Y$,
respectively.

######  Theorem 2.

Let $X\to Y$ be the correct causal direction that follows

$\displaystyle Y$ $\displaystyle=\sum\nolimits_{i=1}^{X}\xi_{i}+\epsilon,$
$\displaystyle\quad X,\xi_{i},\textrm{ and }\epsilon\textrm{~{}are
independent},$ (B.1)

where $\xi_{i}\sim\operatorname{Pois}(\alpha_{X,Y})$,
$\epsilon\sim\operatorname{Pois}(\mu_{Y})$,
$X\sim\operatorname{Pois}(\mu_{X})$. Then, there does not exist a backward
model that admits the following equation:

$\displaystyle X$
$\displaystyle=\sum\nolimits_{i=1}^{Y}\hat{\xi}_{i}+\hat{\epsilon},$
$\displaystyle\quad Y,\hat{\xi}_{i},\textrm{ and }\hat{\epsilon}\textrm{~{}are
independent},$ (B.2)

where $\hat{\xi}_{i}\sim\operatorname{Pois}(\hat{\alpha}_{Y,X})$,
$\hat{\epsilon}\sim\operatorname{Pois}(\hat{\mu}_{X})$,
$Y\sim\operatorname{Pois}(\hat{\mu}_{Y})$.

###  B.1 Proof by Likelihood Function

###### Proof.

Let $\pi(x,y):=\log p(x,y)$ denote the log-likelihood of the causal model. We
will prove by contradiction that there does not exist a backward model that
has the same log-likelihood as the causal direction, i.e., distribution
equivalent.

Suppose that there exists a backward model that has the same log-likelihood as
the causal direction. Then based on the model given in Eq. B.1, the log-
likelihood of the reversed direction can be written as follows. For the causal
direction,

$\displaystyle\pi(x,y)$ $\displaystyle=\log p(x)+\log p(y|x)$ (B.3)
$\displaystyle=x\log\mu_{X}-\mu_{X}-\log
x!+y\log(x\alpha_{X,Y}+\mu_{Y})-x\alpha_{X,Y}-\mu_{Y}-\log y!,$

and for the reverse direction,

$\displaystyle\pi(x,y)$ $\displaystyle=\log p(y)+\log p(x|y)$ (B.4)
$\displaystyle=y\log\hat{\mu}_{Y}-\hat{\mu}_{Y}-\log
y!+x\log(y\hat{\alpha}_{Y,X}+\hat{\mu}_{X})-y\hat{\alpha}_{Y,X}-\hat{\mu}_{X}-\log
x!.$

Let $\Delta_{x}\pi(x,y):=\pi(x+1,y)-\pi(x,y)$,
$\Delta_{x}^{2}\pi(x,y):=\Delta_{x}\pi(x+1,y)-\Delta_{x}\pi(x,y)$ denotes the
first and the second order of the difference of $X$, respectively.

If there exists a backward model, then Eq. B.4 holds, implying:

$\Delta_{x}\pi(x,y)=\log(y\hat{\alpha}_{Y,X}+\hat{\mu}_{X})-\log(x+1)$ (B.5)

and the second order of difference:

$\Delta_{x}^{2}\pi(x,y)=-\log(x+2)+\log(x+1)$ (B.6)

Using Eq. B.3, we have

$\Delta_{x}\pi(x,y)=\log\mu_{X}-\log(x+1)+y\log(\alpha_{X,Y}(x+1)+\mu_{Y})-y\log(\alpha_{X,Y}(x+1)+\mu_{Y})-\alpha_{X,Y},$
(B.7)

and

$\Delta^{2}_{x}\pi(x,y)=-\log(x+2)+\log(x+1)+y\log(\alpha_{X,Y}(x+2)+\mu_{Y})-2y\log(\alpha_{X,Y}(x+1)+\mu_{Y})+y\log(x\alpha_{X,Y}+\mu_{Y}).$
(B.8)

Combining Eq. B.6 and Eq. B.8, yields

$y\log(\alpha_{X,Y}(x+2)+\mu_{Y})-2y\log(\alpha_{X,Y}(x+1)+\mu_{Y})+y\log(x\alpha_{X,Y}+\mu_{Y})=0$
(B.9)

for all $x,y$ holds. The necessary condition for Eq. B.9 holds for all
$x,y\geq 0$ is that $\alpha_{X,Y}=0$ which contradicts the model assumption
that $\alpha_{X,Y}\neq 0$. This completes the proof. ∎

###  B.2 Proof by Probability Generating Function

###### Proof.

If there exists a backward causal model following Eq. B.2, then $Y$ must be
the Poisson distribution, otherwise, the distribution will not be equivalent,
and a simple distribution test would identify the causal direction. Thus, to
show the identifiability of such a model, we only need to prove that the
distribution of $Y$ can not be Poisson.

Let $\Psi_{X}(s)=E\left(s^{X}\right)$ be the probability generating function
(PGF) of the discrete random variable $X$, where $s$ belongs to some interval
containing 1. Specifically, for $X\sim\operatorname{Pois}(\mu_{X})$, the PGF
of $X$ has the form $\Psi_{X}(s)=\exp[\mu_{X}(s-1)]$, and similarly
$\Psi_{\xi_{i}}(s)=\exp[\alpha_{X,Y}(s-1)]$,
$\Psi_{\epsilon}(s)=\exp[\mu_{Y}(s-1)]$. Then, based on the causal model in
Eq. B.1, the probability generating function of $Y$ can be written as follows:

$\Psi_{Y}(s)=\Psi_{X}(\Psi_{\xi_{i}}(s))\Psi_{\epsilon}(s),$ (B.10)

which yields

$\Psi_{Y}(s)=\exp[\mu_{X}(\exp[\alpha_{X,Y}(s-1)]-1)]\exp[\mu_{Y}(s-1)],$
(B.11)

Suppose the backward model holds, and the desired Poisson PGF of $Y$ is that

$\Psi_{Y}(s)=\exp[\hat{\mu}_{Y}(s-1)],$ (B.12)

and the necessary condition for Eq. B.11 and Eq. B.12 has the same form for
all $s$ is that $\mu_{X}=0$ or $\alpha_{X,Y}=0$. If $\mu_{X}=0$, the causal
variable $X$ is a constant that contradicts the model assumption. Similarly,
$\alpha_{X,Y}=0$ also contradicts the assumption that $\alpha_{X,Y}\neq 0$.
Thus, Eq. B.11 and Eq. B.12 do not have the same form and $Y$ can not be
Poisson, which completes the proof. ∎

##  Appendix C Proof of Theorem 3

######  Theorem 3.

With the causal faithfulness assumption and causal sufficiency assumption, the
multivariate instantaneous causal structure is identifiable.

###### Proof.

With the causal faithfulness and the causal sufficiency assumption, we can
identify the causal structure through conditional independence or the sparsity
up to the Markov equivalent class since all Markov equivalent classes share
the same skeleton. We therefore only need to show that any graph that has the
same skeleton will not admit another causal structure that has a different
causal direction while having the same likelihood (i.e., distribution
equivalent). To show this, for a multivariate instantaneous causal structure,
the log-likelihood of the causal graph $\displaystyle\mathcal{G}$ with
parameters $\displaystyle\Theta$ is given as follows:

$L(\mathcal{G},\Theta;\mathbf{X})=\sum_{v\in\mathbf{V}}\log
p(X_{v}|X_{\mathbf{Pa}^{\mathcal{G}}_{v}}).$ (C.1)

(a) Correct causal pair in graph $\mathcal{G}$

(b) Reversed causal pair in graph $\hat{\mathcal{G}}$

Figure 5: An example for two causal graphs that has only one reversed causal
pair. We let $S=\\{v_{1},v_{2}\\}$ denote the set of reversed nodes in
$\hat{\mathcal{G}}$ compared with $\mathcal{G}$

Because the likelihood can be decomposed according to the graph, we can
instead analyze the identifiability of each causal pair as shown in Fig. 5,
that for the correct causal pair in Fig. 5(a) will not admit the reversed
causal direction in Fig. 5(b). Thus, the full likelihood will also not admits
another causal graph among the Markov equivalent class.

The proof is similar to Theorem 2, and we prove it by contradiction. Without
loss of generality, as shown in Fig. 5, we only focus on the graph that only
one causal pair has reversed edges denoted by
$\displaystyle\hat{\mathcal{G}}$, while we denote the correct causal graph as
$\displaystyle\mathcal{G}$. The main difference between these two graphs is
the reversed nodes $\displaystyle S=\left\\{v_{i}|Y\rightarrow
v_{i}\in\hat{\mathcal{G}}\land Y\rightarrow
v_{i}\not{\in}\mathcal{G}\right\\}$ such that
$\displaystyle\mathbf{Pa}^{\mathcal{G}}_{Y}=\mathbf{Pa}^{\hat{\mathcal{G}}}_{Y}\cup
S$ and
$\displaystyle\mathbf{Ch}^{\mathcal{G}}_{S}=\mathbf{Ch}^{\hat{\mathcal{G}}}_{S}\cup
Y$. For example, in Fig. 5, $\displaystyle S=\\{X_{1},X_{2}\\}$.

Suppose that the instantaneous causal structure holds with graph
$\displaystyle\hat{\mathcal{G}}$ such that $\displaystyle
L(\hat{\mathcal{G}},\hat{\Theta};\mathbf{X})=L(\mathcal{G},\Theta;\mathbf{X})$
then the log-likelihood of $\displaystyle\hat{\mathcal{G}}$ can be written as
follows:

$\pi(X_{v_{1}},X_{v_{2}},...,X_{v_{n}}):=\log
p(X_{v_{1}},X_{v_{2}},...,X_{v_{n}})=\sum_{v\in\mathbf{V}}\log
p(X_{v}|X_{\mathbf{Pa}^{\hat{\mathcal{G}}}_{v}}).$ (C.2)

For some set $\displaystyle S=\\{v_{i},v_{j},...\\}$, we denote the first-
order difference of $\displaystyle\pi$ as follows:

$\Delta_{S}\pi(X_{v_{1}},X_{v_{2}},...,X_{v_{n}})=\pi(X_{v_{1}},...,X_{v_{i}}+1,...,X_{v_{j}}+1,...,X_{v_{n}})-\pi(X_{v_{1}},...,X_{v_{i}},...,X_{v_{j}},...,X_{v_{n}}),$
(C.3)

and the second-order difference of $\displaystyle\pi$:

$\Delta^{2}_{S}\pi(X_{v_{1}},X_{v_{2}},...,X_{v_{n}})=\Delta_{S}\pi(X_{v_{1}},...,X_{v_{i}}+1,...,X_{v_{j}}+1,...,X_{v_{n}})-\Delta_{S}\pi(X_{v_{1}},...,X_{v_{i}},...,X_{v_{j}},...,X_{v_{n}}).$
(C.4)

Then, by taking the first-order difference on $\displaystyle S$, we have

$\displaystyle\Delta_{S}\pi(X_{v_{1}},X_{v_{2}},\dotsc,X_{v_{n}})=$
$\displaystyle\sum_{S_{i}\in S}[\log P(X_{S_{i}}+1\mid
X_{\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S},X_{\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap S}+1)-\log P(X_{S_{i}}\mid
X_{\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}})]$ (C.5)
$\displaystyle+\underbrace{\sum_{v\in\mathbf{Ch}_{S}^{\hat{\mathcal{G}}}\backslash
S}[\log P(X_{v}\mid X_{\mathbf{Pa}_{v}^{{\mathcal{G}}}\backslash
S},X_{\mathbf{Pa}_{v}^{{\mathcal{G}}}\cap S}+1)-\log P(X_{v}\mid
X_{\mathbf{Pa}_{v}^{{\mathcal{G}}}})]}_{:=R_{1}}$ $\displaystyle=$
$\displaystyle\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$ $\displaystyle-
X_{S_{i}}\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}\right]$
$\displaystyle-\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}-\log(X_{S_{i}}+1)\Biggr{\\}}+R_{1}$

Note that in the first equality, because only one causal pair has reversed
edges, the parents of variables in
$\displaystyle\mathbf{Ch}^{\hat{\mathcal{G}}}_{S_{i}}\backslash S$ are the
same as in graph $\displaystyle\mathcal{G}$, and therefore we write
$\displaystyle\mathbf{Pa}^{\mathcal{G}}_{v}$ instead of
$\displaystyle\mathbf{Pa}^{\hat{\mathcal{G}}}_{v}$. Similarly, in the second
equality, since there is only one causal pair difference, we have
$\displaystyle\mathbf{Pa}^{\hat{\mathcal{G}}}_{S_{i}}\cap
S=\mathbf{Pa}^{\mathcal{G}}_{S_{i}}\cap S$.

Furthermore, we have

$\displaystyle\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{i}}+1,\dotsc,X_{v_{j}}+1,X_{v_{n}})$
(C.6) $\displaystyle=\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+2)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+2)\right]$
$\displaystyle-(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$
$\displaystyle-\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}-\log(X_{S_{i}}+2)\Biggr{\\}}$
$\displaystyle+\underbrace{\sum_{v\in\mathbf{Ch}_{S}^{\hat{\mathcal{G}}}\backslash
S}[\log P(X_{v}\mid X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+2)-\log P(X_{v}\mid
X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+1)]}_{:=R_{2}}$

then for the second-order difference, we obtain

$\displaystyle\Delta_{S}^{2}\pi(X_{v_{1}},X_{v_{2}},\dotsc,X_{v_{n}})$
$\displaystyle=\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+2)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{t,v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+2)\right]$ (C.7)
$\displaystyle-2(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$
$\displaystyle+X_{S_{i}}\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}\right]$
$\displaystyle-\log(X_{S_{i}}+2)+\log(X_{S_{i}}+1)\Biggr{\\}}+R_{2}-R_{1}$

For the causal direction,

$\displaystyle\Delta_{S}\pi(X_{v_{1}},X_{v_{2}},\dotsc,X_{v_{n}})=$
$\displaystyle\sum_{S_{i}\in S}[\log P(X_{S_{i}}+1\mid
X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap S}+1)-\log P(X_{S_{i}}\mid
X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}})]$ (C.8)
$\displaystyle+\sum_{v\in\mathbf{Ch}_{S}^{\mathcal{G}}\backslash S}[\log
P(X_{v}\mid X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+1)-\log P(X_{v}\mid
X_{Pa_{v}^{\mathcal{G}}})]$ $\displaystyle=$ $\displaystyle\sum_{S_{i}\in
S}[\log P(X_{S_{i}}+1\mid X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap S}+1)-\log P(X_{S_{i}}\mid
X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}})]$
$\displaystyle+\underbrace{\sum_{v\in\mathbf{Ch}_{S}^{\hat{\mathcal{G}}}\backslash
S}[\log P(X_{v}\mid X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+1)-\log P(X_{v}\mid
X_{Pa_{v}^{\mathcal{G}}})]}_{:=R_{1}}$ $\displaystyle+[\log P(X_{Y}\mid
X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\cap S}+1)-\log P(X_{Y}\mid
X_{\mathbf{Pa}_{Y}^{\mathcal{G}}})]$

where in the second equality, because
$\displaystyle\mathbf{Ch}^{\mathcal{G}}_{S_{i}}=\mathbf{Ch}^{\hat{\mathcal{G}}}_{S_{i}}\cup
Y$, we can decompose the sum into
$\displaystyle\mathbf{Ch}^{\hat{\mathcal{G}}}_{S_{i}}$ and $\displaystyle Y$.

Similarly, we have

$\displaystyle\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{i}}+1,\dotsc,X_{v_{j}}+1,\dotsc,X_{v_{n}})$
(C.9) $\displaystyle=$ $\displaystyle\sum_{S_{i}\in S}[\log p(X_{S_{i}}+2\mid
X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap S}+2)-\log p(X_{S_{i}}+1\mid
X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap S}+1)]$
$\displaystyle+\underbrace{\sum_{v\in\mathbf{Ch}_{S}^{\hat{\mathcal{G}}}\backslash
S}[\log P(X_{v}\mid X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+2)-\log P(X_{v}\mid
X_{\mathbf{Pa}_{v}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{v}^{\mathcal{G}}\cap S}+1)]}_{:=R_{2}}$
$\displaystyle+[\log P(X_{Y}\mid X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\cap S}+2)-\log P(X_{Y}\mid
X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S},X_{\mathbf{Pa}_{Y}^{\mathcal{G}}\cap S}+1)]$

By expanding Eq. C.8 and Eq. C.9, we have

$\displaystyle\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{n}})$
$\displaystyle=\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$ (C.10) $\displaystyle-
X_{t,S_{i}}\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}\right]$
$\displaystyle-\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}-\log(X_{S_{i}}+1)\Biggr{\\}}+R_{1}$
$\displaystyle+\Biggl{\\{}X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}(X_{v^{\prime}}+1)\right]$ $\displaystyle-
X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}\right]-\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}-\log X_{Y}\Biggr{\\}},$

and

$\displaystyle\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{i}}+1,\dotsc,X_{v_{j}}+1,\dotsc,X_{v_{n}})$
(C.11) $\displaystyle=$ $\displaystyle\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+2)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+2)\right]$
$\displaystyle-(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\hat{\mathcal{G}}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$
$\displaystyle-\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}-\log(X_{S_{i}}+2)\Biggr{\\}}+R_{2}$
$\displaystyle+\Biggl{\\{}X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+2)\right]$ $\displaystyle-
X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]-\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}-\log X_{Y}\Biggr{\\}},$

respectively.

Then for the second-order difference, we have

$\displaystyle\Delta_{S}^{2}\pi(X_{v_{1}},X_{v_{2}},\dotsc,X_{v_{n}})$ (C.12)
$\displaystyle=\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{i}}+1,\dotsc,X_{v_{j}}+1,\dotsc,X_{v_{n}})-\Delta_{S}\pi(X_{v_{1}},\dotsc,X_{v_{i}},\dotsc,X_{v_{j}},\dotsc,X_{v_{n}})$
$\displaystyle=\sum_{S_{i}\in
S}\Biggl{\\{}(X_{S_{i}}+2)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+2)\right]$
$\displaystyle-2(X_{S_{i}}+1)\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}(X_{v^{\prime}}+1)\right]$
$\displaystyle+X_{S_{i}}\log\left[\mu_{S_{i}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{S_{i}}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},s_{i}}X_{v^{\prime}}\right]$
$\displaystyle-\log(X_{S_{i}}+2)+\log(X_{S_{i}}+1)\Biggr{\\}}+R_{2}-R_{1}$
$\displaystyle+\Biggl{\\{}X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}(X_{v^{\prime}}+2)\right]$
$\displaystyle-2X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}(X_{v^{\prime}}+1)\right]$
$\displaystyle+X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}\right]\Biggr{\\}}$

If the model is not identifiable, the second-order difference should have the
same value in both causal direction and the reverse direction. Thus, combining
Eq. C.7 and Eq. C.12, we obtain the following equation

$\displaystyle
X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}(X_{v^{\prime}}+2)\right]$ (C.13)
$\displaystyle-2X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}(X_{v^{\prime}}+1)\right]+X_{Y}\log\left[\mu_{Y}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\backslash
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}+\sum_{v^{\prime}\in\mathbf{Pa}_{Y}^{\mathcal{G}}\cap
S}\alpha_{v^{\prime},Y}X_{v^{\prime}}\right]=0$

which must hold for all possible values $\displaystyle X_{v}\in\mathbb{N}$.
The necessary condition for Eq. C.13 holds is that for all $\displaystyle
v^{\prime}\in S$, $\displaystyle\alpha_{v^{\prime},Y}=0$ which contradicts to
the assumption that $\alpha_{v^{\prime},Y}\neq 0$. This finishes the proof. ∎

##  Appendix D Additional Experiments

The main paper has shown the F1 scores and other baselines in both synthetic
data and real-world experiments. Here, we further provide the Precision,
Recall, and Structural Hamming Distance (SHD) in these experiments, as shown
in Fig. 6 and Fig. 7, Fig. 8, and Fig. 9. Note that most of the parameters are
based on the default setting in the tick packages Bacry et al. [2018].

(a) Sensitivity to Time Interval

(b) Sensitivity to Range of $\alpha$

(c) Sensitivity to Range of $\mu$

(d) Sensitivity to Sample Size

(e) Sensitivity to Num. of Event Types

(f) Sensitivity to Avg. Indegree

Figure 6: Precision in the Sensitivity Experiments

(a) Sensitivity to Time Interval

(b) Sensitivity to Range of $\alpha$

(c) Sensitivity to Range of $\mu$

(d) Sensitivity to Sample Size

(e) Sensitivity to Num. of Event Types

(f) Sensitivity to Avg. Indegree

Figure 7: Recall in the Sensitivity Experiments

(a) Sensitivity to Time Interval

(b) Sensitivity to Range of $\alpha$

(c) Sensitivity to Range of $\mu$

(d) Sensitivity to Sample Size

(e) Sensitivity to Num. of Event Types

(f) Sensitivity to Avg. Indegree

Figure 8: SHD in the Sensitivity Experiments

(a) Precision

(b) Recall

(c) SHD

Figure 9: Real-World Experiment on Different Temporal Resolutions