# Bayesian Estimation of Multivariate Hawkes Processes with Inhibition and
Sparsity

###### Abstract

Hawkes processes are point processes that model data where events occur in
clusters through the self-exciting property of the intensity function. We
consider a multivariate setting where multiple dimensions can influence each
other with intensity function to allow for excitation and inhibition, both
within and across dimensions. We discuss how such a model can be implemented
and highlight challenges in the estimation procedure induced by a potentially
negative intensity function. Furthermore, we introduce a new, stronger
condition for stability that encompasses current approaches established in the
literature. Finally, we examine the total number of offsprings to
reparametrise the model and subsequently use Normal and sparsity-inducing
priors in a Bayesian estimation procedure on simulated data.

Keywords: Intensity Function, Interpretation, Prior Choice, Stability

##  1 Introduction

Classic Hawkes processes (Hawkes,, 1971) are point processes capturing the
self-exciting behaviours of events that occur in clusters or bursts. A
prominent application is found in the earthquake literature, where large
earthquakes trigger subsequent aftershocks (Ogata,, 1988). Hawkes processes
are successfully used in a variety of other settings that share this self-
excitation nature: crime and terror modelling (Mohler,, 2013; Shelton et al.,,
2018; Tucker et al.,, 2019), accidents (Kalair et al.,, 2021), stock market
tradings (Rambaldi et al.,, 2017), and social media (Lai et al.,, 2016).

For each event from a point process process we record the time when it
happened $Y=\left(t_{1}\dots t_{N}\right)\in\left[0,T_{max}\right]$, such that
$t_{i}$ is the event time at which the $i^{th}$ event took place. The self-
exciting linear Hawkes process is defined by its conditional intensity
function $\lambda(t|Y^{t},\Theta)$, which at time $t$ is conditional on the
previous events $Y^{t}=\\{t_{i}:t>t_{i}\\}$ (Hawkes,, 1971) and parameters
$\Theta$ used to specify the parametric form of intensity. For convenience of
notation the dependence on $Y^{t}$ and $\Theta$ is suppressed further on.

The following specification encompasses the self-exciting nature:

$\lambda(t|Y^{t},\Theta)=\mu(t)+\sum_{i:t>t_{i}}K\,g(t-t_{i})$ (1)

Here, $\mu(\cdot)>0$ is the background rate that can capture seasonality and
underlying trends. We call $K\,g(\cdot)$ the excitation kernel, where $K\geq
0$, $g(x)\geq 0$ for $x\geq 0$ and $\int_{0}^{\infty}g(x)\,dx=1$. Each
observation prior to $t$ contributes to the intensity at time $t$ as governed
by the kernel. This drives the self-exciting behaviour of the Hawkes process.
Both the background rate and the kernel depend on parameters
$\Theta=(\Theta_{\mu},K,\Theta_{g})$ where $\Theta_{\mu}$ and $\Theta_{g}$
contains all parameters from $\mu(\cdot)>0$ and $g(\cdot)$, respectively. To
estimate them, a plethora of methods is successfully employed in the
literature, for example maximum likelihood approaches and Bayesian methods
(Veen and Schoenberg,, 2008; Rasmussen,, 2013; Guo et al.,, 2018; Chen and
Stindl,, 2018).

###  1.1 Paper Overview

The remainder of this paper is organised as follows. Section 2 reviews the
$M$-dimensional Hawkes process with corresponding sparsity and the branching
structure interpretation. In Section 3 we describe the multivariate Hawkes
process with inhibition alongside considerations to ensure a non-negative
intensity. Section 4 discusses two aspects regarding the model implementation.
First, we introduce a new, stronger conditions for stability, which is
compared to two current approaches from the literature. Second, we highlight
challenges when integrating the intensity function. In Section 5 we examine
Normal and sparsity priors based on a reparametrisation of the multivariate
Hawkes processes. Section 6 summarises the posterior inference and Section 7
presents an example on simulated data using the aforementioned aspects. We
finish with a discussion and comments on future work in Section 8.

##  2 Background

This section considers the $M$-dimensional Hawkes process and summarises
approaches for sparse estimation procedures. Moreover, we examine how the
Hawkes process can be written as a superposition of Poisson processes that
allow for a branching interpretation.

###  2.1 Multivariate Hawkes Process

The model presented in Equation 1 extends to an $M$-dimensional linear Hawkes
process incorporating both self-excitation and cross-excitation. This is
called a multivariate linear Hawkes process. Assume that there are $M$
dimensions with event times $Y_{1}=\left(t_{1\,1}\dots t_{1\,N_{1}}\right)$ in
dimension $1$ to event times $Y_{M}=\left(t_{M\,1}\dots t_{M\,N_{M}}\right)$
in dimension $M$. At time $t$ the intensity in dimension $i$ is :

$\lambda_{i}(t)=\mu_{i}(t)+\sum_{j=1}^{M}\sum_{l:t>t_{j\,l}}K_{ji}\,g_{ji}(t-t_{j\,l})$
(2)

For the remainder of this paper we consider the background rate to be constant
in each dimension: $\mu_{i}(x)=\mu_{i}$ and write
$\mathbf{\mu}=(\mu_{1}\dots\mu_{M})$. Moreover, we assume the following form
for the excitation kernel for all $i,j$: $g_{ij}(x)\geq 0$ for $x\geq 0$ and
$\int_{0}^{\infty}g_{ij}(x)\,dx=1$. Here, $K_{ij}\geq 0$ describes the
excitation effect an event in dimension $i$ has on dimension $j$. We write
$\mathbf{K}=\\{K_{ij}\\}$ where $i,j=1\dots M$.

For any multivariate point process the likelihood for data $Y$ is (Daley and
Vere-Jones,, 2003, p. 23):

$p(Y|\Theta)=\prod_{m=1}^{M}\prod_{i=1}^{N_{m}}\lambda_{m}(t_{i}|\theta)\,exp\left(-\Lambda\right)$
(3)

where $\Lambda=\sum_{m=1}^{M}\int_{0}^{T_{max}}\lambda_{m}(x|\Theta)\,dx$ and
$\Theta=(\mu,\mathbf{K},\Theta_{g})$. Here, $\Theta_{g}$ contains all
parameters from $g_{ij}(\cdot)$ for all $i,j=1\dots M$. Note that the
evaluation of both the likelihood and the log-likelihood require $\Lambda$,
which is the integral of the intensity function.

###  2.2 Sparsity

Many real life problems of multivariate Hawkes processes also contain sparsity
in $\mathbf{K}$ (e.g. Eichler et al.,, 2017). For our example (Section 7) we
make use of this concept through sparsity-inducting priors (Section 5.3).
There are multiple successful cases of sparsity application available in the
literature. Xu et al., (2016) use group Lasso and the EM-Algorithm to uncover
the entries of $\mathbf{K}$ that are truly zero. Two-step procedure are a
common workaround to uncover those structural zeros. Yu et al., (2020)
separate the learning of the structure and the estimation of the non-zero
parameters. Similarly, Miscouridou et al., (2018) use elaborate hierarchical
structures and prior specifications to uncover sparsity, heterogeneity,
reciprocity, and community structure in their multivariate Bayesian estimation
procedure. Both methods only provide an approximation to the parameters (or
their posterior distributions).

###  2.3 Branching Interpretation

We now examine how data from a one-dimensional linear Hawkes process can be
related to an underlying branching structure as we utilise this concept in
Section 5.1 to motivate a new parametrisation. We discuss here the common
branching structure interpretation when $K>0$.

The model in Equation 1 can be written as the superposition of Poisson
processes (Hawkes and Oakes,, 1974) such that the intensity function is a sum
of independent Poisson processes. Suppose that $j$ events $(t_{1}\dots t_{j})$
happened before time $t$. Then the intensity at time $t$ is the sum of the
background process $\mu(t)$ and $j$ offspring processes with intensities
$K\,g(t-t_{j})$, where each offspring process was triggered by a previous
event. This gives rise to the following branching structure interpretation of
a Hawkes process (Daley and Vere-Jones,, 2003, p. 202). Data generated from an
excitation-only Hawkes process consists of two types of events that come from
distinct Poisson processes, where each event is generated by exactly one
process (Rasmussen,, 2013):

  1. 1.

Immigrant events come from the background process with intensity $\mu(t)$.

  2. 2.

Offspring events come from an offspring process which had been triggered by a
previous event. Here, each event has an average of $K$ direct offsprings if
$T_{max}\to\infty$.

Note that offspring events trigger offspring processes as well. This can lead
to cascades started by an immigrant event that has offspring events, which, in
turn, has offspring events etc. Figure 1 visualises this branching structure.
This concept will be utilised in Section 5.1 to introduce a new
parametrisation based on the total number of offsprings.

Figure 1: Sketch of the branching structure interpretation in one dimension.
In this example $(t_{1},t_{6},t_{9})$ are the immigrant events, indicated by a
star on the top and the summarised timeline below. The event $t_{1}$ has two
direct offsprings, namely $t_{2}$ and $t_{3}$. The former has further
offsprings. For example $t_{1}$ has a total of five offsprings (number of
offspring events in the cascade that started in $t_{1}$), whereas $t_{9}$ has
none.

While it is not known whether an events is an immigrant of offspring for a
given data set, this perspective nevertheless can be exploited for sampling,
inference, and interpretation. These interpretations hold asymptotically for
$T_{max}\to\infty$. For finite $T_{max}$ some offspring events may be larger
than $T_{max}$ and therefore would not be included in a simulated data set.
Such edge effects are common in the Hawkes literature and diminish when
$T_{max}$ is large (Daley and Vere-Jones,, 2003, p. 275).

This interpretation extends to the excitation multivariate case (all entries
of $\mathbf{K}$ positive) where an event in dimension $i$ triggers offspring
processes in all dimensions $j$. In turn, resulting offspring events can
trigger offspring processes in all dimensions. Hence, cascades can potentially
lead through different dimensions, which all contribute to the total number of
offsprings of the immigrant event that started the cascade, which is discussed
in Section 5.1.

##  3 Multivariate Hawkes Process with Inhibition

This section discusses the multivariate Hawkes process under both excitation
and inhibition. After the model set up we review approaches to ensure a non-
negative intensity function and focus on the specifications of a link
function.

While most examples of a Hawkes process just modelling excitation (each
$K_{ij}\geq 0$), Hawkes processes also capture both excitation and inhibition
when each $K_{ij}\in\mathbb{R}$. Hence, the occurrence of an event in
dimension $i$ can increase or decrease the intensity function in each
dimension $j=1\dots M$. This is a substantial extension of the specification
outlined in Equation 2, which has received less attention in the literature.

A prominent application of Hawkes processes with inhibition can be found in
neural spike trains, where the inhibition captures the period of decreased
intensity after a neural spike event (Eichler et al.,, 2017). Moreover, such
models have been successful employed in a variety of fields, such as financial
applications (Lu and Abergel,, 2018; Tetereva,, 2018) or in adequately
modelling the popular MemeTracker data set (Lemonnier and Vayatis,, 2014).

Other approaches for inhibition are also suggested in the literature, for
example Apostolopoulou et al., (2019) base their approach on a multiplicative
thinning factor, which leads to a separate strand of work beyond of the
standard Hawkes Process set up. For this paper we only consider inhibition
through negative $K_{ij}$.

###  3.1 Model

The intensity of a $M$-dimensional Hawkes process with inhibition at time $t$
in dimension $i$ is:

$\lambda_{i}(t)=\mu_{i}+\sum_{j=1}^{M}\sum_{l:t>t_{j\,l}}K_{ji}\,g_{ji}(t-t_{j\,l})$
(4)

Here, each $K_{ij}\in\mathbb{R}$, and $\mathbf{K}=\\{K_{ij}\\}$ where
$i,j=1\dots M$. When $K_{ij}>0$ we talk about excitation as above, when
$K_{ij}<0$ we call this inhibition. However, any point process requires, by
design, a non-negative intensity function at every $t\in[0,T_{max}]$. We
therefore now discuss two approaches to ensure this in Section 3.2.

###  3.2 Non-Negative Intensity

As the $K_{ij}$ in Equation 4 can be negative it is not guaranteed that the
intensity always stays above zero. There are two main approaches in the
literature to tackle this issue to ensure non-negativity: restricting the
parameter space or using a link function to ensure positivity.

####  3.2.1 Restricting the Parameter Space

We can restrict the parameter space to those parameters that lead to a
positive intensity function for the observed data set. In a Bayesian frame
work this can be incorporated into the prior. However, this method has three
fundamental shortcomings: potential non-consistency, sampling issues, and
data-dependency.

Firstly, in many applications we see an intensity function that exhibits many
spikes and drops with short, rapid changes. Excluding all parameter
combinations that lead to a decrease below zero somewhere would substantially
restrict the parameter space, where many parameter combinations are
prohibited. This could potentially exclude the true parameters if inhibition
is truly present, which would lead to a non-consistent estimation procedure.
Secondly, restricting the parameter space can cause sampling issues as MCMC
samplers are prone to boundary artefacts under such conditions. Thirdly, the
‘permissibility’ of a parameter combination depends on the data. A set of
parameters causing a non-negative intensity on a data set does not guarantee
this property when more data is collected.

####  3.2.2 Link Function

Instead, a link function $\phi(\cdot)$ is commonly used to ensure a non-
negative intensity, which leads to the following intensity:

$\lambda_{i}(t)=\phi\left(\mu_{i}+\sum_{j=1}^{M}\sum_{l:t>t_{j\,l}}K_{ji}\,g_{ji}(t-t_{j\,l})\right)$
(5)

The link function $\phi(\cdot):\mathbb{R}\to\mathbb{R}^{+}$ is needed to
ensures a non-negative intensity at every $t$. For example, Mei and Eisner,
(2017) use the softplus function $\phi(x)=s\,\,log(1+exp(x/s))$ with parameter
$s$. One straightforward choice of $\phi$ is the ReLU function where
$\phi(x)=\max(a,x)$ for a small, non-negative $a$. A popular approach in the
literature is to set $a=0$ (Lemonnier and Vayatis,, 2014; Lu and Abergel,,
2018; Costa et al.,, 2020). Crucially, this choice of a link function
preservers the interpretation of $K_{ij}$ as the average number of direct
offsprings (as described in Section 2.3). Moreover, we also shown in Appendix
C that it retains the branching structure, even under inhibition, which is a
novel contribution to the field. Hence, for the remainder of this paper we use
the ReLU function with $\phi(x)=\max(0,x)$.

##  4 Implementation

The model introduced in Equation 5 with $\phi(\cdot)=\max(0,x)$ is
straightforward to define, but difficult to implement as two fundamental
challenges need to be overcome. First, we review approaches from the
literature to assess stability (as defined in Section 4.1) and subsequently
introduce a new, stronger condition. Second, we discuss the challenge of
integrating the intensity function under the given model.

###  4.1 Stability

Due to the self-exciting behaviour of the Hawkes process it is possible that
an infinite number of events take place in finite time. For example, this can
happen in a one-dimensional Hawkes-process if each event has more on average
one or more offsprings. This behaviour is called supercritical (Helmstetter
and Sornette,, 2002) or explosive (Browning et al.,, 2021). However, for real-
life applications it can be desirable to limit the parameter space to non-
explosive instances (Kolev and Ross,, 2019). This is referred to as stability
(Bremaud and Massoulie,, 1996; Bacry et al.,, 2020), i.e. there exists a
unique stationary distribution with finite average intensity (Sulem et al.,,
2021).

Two conditions (LABEL:c1eigen, LABEL:c2colsum defined below) have been used in
the literature to determine stability, as outlines in Section 4.1.1. In
Section 4.1.2 we introduce a new condition that is stronger than both of the
previously used ones. We proof that when at least one of (LABEL:c1eigen,
LABEL:c2colsum) hold, so does our condition. Moreover, there exist parameters
$\mathbf{K}$ for which our condition holds when neither of (LABEL:c1eigen,
LABEL:c2colsum) do. This permits us to have a unified approach for checking
stability and to classify a larger set of parameters as stable. We also
provide a toy-example to illustrate this usefulness of our suggested
condition.

We introduce the following notations: $\operatorname{\text{abs}}(\mathbf{A})$
is an $M\times M$ matrix where each entry is $|A_{ij}|$, the absolute value of
$A_{ij}$. Moreover, $\mathbf{A}^{+}$ is the matrix with entries
$\max(A_{ij},0)$. We write $\rho(\mathbf{A})$ for the spectral radius of
matrix $\mathbf{A}$, i.e. the largest absolute eigenvalue of matrix
$\mathbf{A}$.

####  4.1.1 Stability Conditions in the Literature

To start, we state two condition to assess stability that have been introduced
in the literature:

  1. C1

A Hawkes process is stable if $\rho(\operatorname{\text{abs}}(\mathbf{K}))<1$
(Bremaud and Massoulie,, 1996).

  2. C2

A Hawkes process is stable if $max_{j}\sum_{i=1}^{M}K^{+}_{ij}<1$ (Sulem et
al.,, 2021).

Note that LABEL:c1eigen uses the spectral radius of the absolute value matrix,
whereas LABEL:c2colsum utilises only the positive part of the matrix
$\mathbf{K}$. Both conditions are sufficient, but not necessary. Sulem et al.,
(2021) discuss both LABEL:c1eigen and LABEL:c2colsum, but do not compare them
as neither is stronger than the other. As LABEL:c1eigen uses the absolute
value matrix, negative entries (inhibition) are converted into excitation.
Hence, while the intensity is decreased by certain events, stability is check
as if those events increased the intensity. Of course, this procedure is
sufficient, but the intensity
$\lambda_{m}(\cdot|\operatorname{\text{abs}}(\mathbf{K}))$ may be rather
different than the original $\lambda_{m}(\cdot|\mathbf{K})$, in particular
when strong inhibition is present. Hence, we found there to be scope to
develop a new criterion that is more closely tailored to processes that
contain inhibition.

####  4.1.2 Introducing a New Condition

We propose a new condition to assess stability for a given parameter
$\mathbf{K}$:

######  Theorem 1 (C3).

If $\rho(\mathbf{K}^{+})<1$, then the process with intensities
$\lambda_{m}(\cdot\,|\,\mathbf{K}$), $m=1\dots M$, is stable.

To prove it, we state an auxiliary lemma.

######  Lemma 2.

Suppose a process with intensity $\bar{\lambda}(\cdot)$ is stable. Then a
process with intensity $\lambda(\cdot)$ such that
$\bar{\lambda}(x)\geq\lambda(x)$ for all $x\geq 0$ is also stable.

###### Proof.

The stable process with intensity $\bar{\lambda}(\cdot)$ has a finite average
intensity and $\bar{\lambda}(x)\geq\lambda(x)$ holds for all $x\geq 0$. The
average intensity of $\lambda(\cdot)$ is therefore at most as large as the
average intensity of $\bar{\lambda}(\cdot)$. Hence, the average intensity of
$\lambda(x)$ is also finite and the process with intensity $\lambda(x)$ is
also stable. ∎

Using this lemma we now prove Theorem 1.

###### Proof.

Since $K^{+}_{ij}\geq K_{ij}$ for each $i,j=1\dots M$, it follows that
$\lambda_{m}(t\,|\,\mathbf{K}^{+})\geq\lambda_{m}(t\,|\,\mathbf{K)}$ for all
$m=1\dots M$ and all $t\in[0,T_{max}]$. By Lemma 2, if a process with
intensity $\lambda_{m}(\cdot\,|\,\mathbf{K}^{+})$ is stable then a process
with intensity $\lambda_{m}(\cdot\,|\,\mathbf{K)}$ is stable as well. A
multivariate process is stable if the intensity in each dimension $m=1\dots M$
has finite average. If $\rho(\mathbf{K}^{+})<1$ then the multivariate process
with intensities $\lambda_{m}(\cdot\,|\,\mathbf{K}^{+})$, $m=1\dots M$, is
stable, and therefore the multivariate process with intensities
$\lambda_{m}(\cdot\,|\,\mathbf{K)}$, $m=1\dots M$, is stable as well. ∎

####  4.1.3 Comparison

We now compare LABEL:c3new to either of the existing conditions
(LABEL:c1eigen, LABEL:c2colsum) and show that if at least one of them holds,
so does LABEL:c3new. Moreover, there are examples where only LABEL:c3new
holds. This implies that LABEL:c3new can confirm stability for more
parameters, which is useful when fitting a multivariate Hawkes process.

######  Theorem 3.

When LABEL:c1eigen holds, then LABEL:c3new holds as well.

We first state the following lemma.

######  Lemma 4.

Let $X$ be a $N\times N$ matrix with entry $X_{ij}$ in row $i$ and colum $j$.
Then $\rho(X)\leq\max_{j}\sum_{i=1}^{N}X_{ij}$.

This is a direct consequence of the Gelfand formula (Gelfand,, 1941).

With that, we can now provide the proof for the stated Theorem 3.

###### Proof.

First we compare LABEL:c3new to LABEL:c1eigen when all entries $K_{ij}$ are
non-negative (i.e. excitation only). It is trivial to see that when
LABEL:c1eigen holds, then LABEL:c3new holds as well since
$\operatorname{\text{abs}}(\mathbf{K})=\mathbf{K}=\mathbf{K}^{+}$ and
therefore $\rho(\mathbf{K}^{+})<1$.

When we do not restrict the entries $K_{ij}$ to non-negative, we note that
each entry of $\operatorname{\text{abs}}(\mathbf{K})$ is as least as large as
the corresponding entry in $\mathbf{K}^{+}$. Within the entry-wise positive
matrices, the spectral radius is monotonous (Lemma 12, p. 153 Serre,, 2002).
Hence, if $\rho(\operatorname{\text{abs}}(\mathbf{K}))<$ then also
$\rho(\mathbf{K}^{+})<1$.

Therefore, when LABEL:c1eigen holds, then LABEL:c3new holds as well. ∎

######  Theorem 5.

When LABEL:c2colsum holds, then LABEL:c3new holds as well.

###### Proof.

By Lemma 4, if $\max_{j}\sum_{i=1}^{N}K_{ij}<1$ then also
$\rho(\mathbf{K}^{+})<1$ and therefore LABEL:c3new holds if LABEL:c2colsum
holds. ∎

Hence, we have shown that if at least one of (LABEL:c1eigen, LABEL:c2colsum)
holds, so does LABEL:c3new. In addition, there are examples where neither of
the existing conditions could confirm stability, but by using LABEL:c3new we
can verify that the process is stable. Let us examine a two-dimensional Hawkes
process with

$\displaystyle\mathbf{K}=\left(\begin{array}[]{cc}0.5&1\\\
-2&0.5\end{array}\right)$

as an illustrative example. Note that
$\rho(\operatorname{\text{abs}}(\mathbf{K}))>1$ and the maximum column sum of
$\mathbf{K}^{+}$ is also larger than $1$, hence neither LABEL:c1eigen nor
LABEL:c2colsum hold. Hence, by just using the two existing conditions it is
not possible to assess whether a process using $\mathbf{K}$ would be stable.
However we can make use of LABEL:c3new, as $\rho(\mathbf{K}^{+})<1$ and
confirm that a process using $\mathbf{K}$ is stable.

In summary, LABEL:c3new not only provides a unified approach to assess
stability, it also permits us to determine stability for more parameters than
by just using (LABEL:c1eigen, LABEL:c2colsum).

###  4.2 Integrating the Intensity

Many parameter estimation methods, such as our subsequent Bayesian procedures,
need to evaluate the likelihood (Equation 3). However, this requires
$\Lambda$, the integral of the intensity function. While this is
straightforward task in an excitation-only Hawkes process, more deliberation
is required for the more general model introduced in Section 3.1. We discuss
different approaches from the literature, present exact solutions for
particular choices of the influence kernel, and suggest an approximate
solution that will be subsequently utilised.

When all $K_{ij}$ are non-negative and $\phi(x)=x$, the integral of the
intensity can be computed by integrating each segment between events (as well
as the ones between $0$ and the first event and between the last event and
$T_{max}$). This reduces to an easy computation of the integral that sums over
the background rate and the contributions of each event (Ogata,, 1981).
However, it becomes difficult to integrate the intensity function under
inhibition in this scenario as the intensity might drop below zero. Lu and
Abergel, (2018) choose to calculate the integral as in the excitation-only
case, even though the parts of the integral below zero contributing negatively
to the integral.

The use of a link function $\phi(\cdot)$ guarantees a non-negative intensity
function (see Section 3.2.2). However, any choice of link that is not the
identity prohibits a straightforward calculation of the integral as the
approach by Ogata, (1981) cannot be employed anymore. Mei and Eisner, (2017)
explain in their supplementary materials that they approximate $\Lambda$ for a
one-dimensional Hawkes process by sampling one $t^{*}$ uniformly from
$[0,T_{max}]$ and then use $\hat{\Lambda}=T_{max}\lambda(t^{*})$. While
$\hat{\Lambda}$ is indeed an unbiased estimator for $\Lambda$ it has a large
variance. The intensity functions often exhibit rapid ups and downs, such that
just evaluating it once cannot capture all aspects of the function. Ertekin et
al., (2015) use Approximate Bayesian Computation to circumvent the problem as
this method does not require an evaluation of the likelihood and therefore the
integral is not required.

Another possibility is to identify the intervals of the intensity function
that are non-negative and integrate those only, which corresponds to setting
$\phi(x)=\max(0,x)$. While this is an exact solution, this approach requires
finding the roots of the intensity function, which depend on the exact
specifications of the model. A popular choice for the influence kernel is the
exponential kernel $g_{ij}(x)=\beta_{ij}\,exp\left(-\beta_{ij}\,x\right)$
(e.g. Blundell et al.,, 2012; Shelton et al.,, 2018). In Appendix A.2 we show
that the roots of the intensity function with this kernel are the solutions to
a high-order polynomial. The computational complexity of this root finding
approach is high since the root-finding needs to be employed between each
event time in every dimension. When setting all $\beta_{ij}$ to be equal, we
find a simple expression for the roots when using the exponential kernel, as
shown in Appendix A.1. However, this assumption may be too restrictive in
application.

For our example we resort to a numerical approximation of $\Lambda$ without
placing any restrictions on $\beta_{ij}$. In each segment between events in
every dimension we use a Simpson’s rule approximation, for details see
Appendix B. In trial runs on simulated data we found this cubic approximation
(‘3/8 Simpson’s rule’) to provide a good approximation, albeit with four
intensity evaluations per segment. While this is a costly approximation in
terms of computational complexity, any lower level approximation introduced
too much bias, and yet the cubic approximation is faster than solving a
polynomial in each interval.

##  5 Prior Choice

To implement the proposed model in a Bayesian framework we investigate prior
choices. Throughout the literature Hawkes processes are parameterised in terms
of $\mathbf{K}$. For example, Browning et al., (2021) choose a uniform prior
for each of its entries. However, when the dimension $M$ is large the non-
negative entries of $\mathbf{K}$ have to be smaller in order to retain
stability in accordance the criteria outlined in Section 4.1. For example,
when all entries of $\mathbf{K}$ are $0.4$, a two-dimensional process is
stable, whereas a three-dimensional $\mathbf{K}$ has an eigenvalue larger than
one. Hence, in a Bayesian framework priors on $\mathbf{K}$ would have to be
adapted according to the dimension $M$.

Instead, we suggest to reparametrise the model to circumvent this problem
using the total number of offsprings, which does not suffer from this
dimension-dependency. In the following sections we define this quantity and
subsequently use it to parametrise the problem, such that normal priors and
sparsity-inducing horseshoe priors can be employed.

###  5.1 Reparametrise the Model Using the Total Number of Offsprings

We now formally define the total number of offsprings. First, we assume that
all entries of $\mathbf{K}$ are non-negative. Section 2.3 introduces the
notion of direct offspring, i.e. an event in dimension $j$ from the immigrant
process of an event in dimension $i$. However, when estimating the parameters
it can become difficult to distinguish between and event in $i$ triggering an
immigrant process in dimension $j$ (‘$i\to j$’) and an event in $i$ triggering
an immigrant process in dimension $k$, which produces an event that, in turn,
triggers an immigrant process in $j$ (‘$i\to k\to j$‘). This problem is
perpetuated further in higher dimensions (Eichler,, 2013). We therefore
propose to investigate the total number of offsprings that circumvents this
issue. For this, two kinds of events are define as an indirect offspring of
event $t_{i}$:

  1. 1.

an event from an immigrant process that was triggered by a direct offspring of
$t_{i}$

  2. 2.

an event from an immigrant process that was triggered by an indirect offspring
of $t_{i}$

We define $K^{*}_{ij}$ as the total number of offsprings an event in dimension
$i$ has in dimension $j$, which is calculated as the sum of direct and
indirect offsprings. We write $\mathbf{K}^{*}=\\{K^{*}_{ij}\\}$ were
$i,j=1\dots M$.

######  Theorem 6.

The total number of offsprings is

$\mathbf{K}^{*}=\left(I-\mathbf{K}\right)^{-1}-I$

The proof is available in Appendix D. It is important to note that the idea
and definition of $\mathbf{K}^{*}$ have already been presented in the
literature, but not in a unified concept towards use in application. On one
hand, Bacry and Muzy, (2016) introduce this formula without its interpretation
in the context of a matrix convolution. On the other hand, Bacry et al.,
(2016) define the concept of total offsprings but do not provide a closed form
expression. Crucially, neither of them make extensive use of the concept.

If a process is stable (see Section 4.1), the calculation remains the same
when entries of $\mathbf{K}$ are negative and $\mathbf{K}^{*}$ still provides
meaningful interpretation. While positive entries of $\mathbf{K}^{*}$ describe
the average number of total offsprings, a negative entry summarises the
negative contributions to the intensity function across dimensions. The number
of actually inhibited events depends on the number of events in the process.
Nevertheless, $\mathbf{K}^{*}$ retains its attractive interpretation and can
be used to place priors without having to consider the dimensions $M$.

Crucially, the entries of $\mathbf{K}^{*}$ for a stable process are not
dependent on the dimension $M$. For the example in Section 7 we therefore
reparameterise the multivariate Hawkes process in terms of $\mathbf{K}^{*}$
such that the intensity for dimension $i$ is written as:

$\lambda_{i}(t)=\phi\left(\mu_{i}+\sum_{j=1}^{M}\sum_{l:t>t_{j\,l}}\left\\{f(\mathbf{K}^{*})\right\\}_{ji}\,g_{ji}(t-t_{j\,l})\right)$
(6)

where $f(\mathbf{X})=I-(\mathbf{X}-I)^{-1}$.

###  5.2 Normal Priors

This reparametresation in Equation 6 permits us to use priors directly for
$\mathbf{K}^{*}$ and write the model already in the quantity of interest. We
restrict the parameter space of $\mathbf{K}^{*}$ such that only stable
parameters (according to LABEL:c3new) are allowed. One option is to place
independent normal priors on each entry of $\mathbf{K}^{*}$:

$\displaystyle K^{*}_{ij}$ $\displaystyle\sim\mathcal{N}(0,1)$
$\displaystyle\text{ for }i,j=1\dots M$ (stable only) (7)

###  5.3 Sparsity Priors

Many real life problems of multivariate Hawkes processes contain sparsity (see
Section 2.2, which we wish to reflect in our model. As mentioned in Section 5,
when $M$ is large, the positive entries of $\mathbf{K}$ have to be small in
order to retain stability (according to LABEL:c3new in Section 4.1). When
employing shrinkage strategies in larger dimensions it can therefore be
difficult to distinguish between small non-zero entries of $\mathbf{K}$ and
truly zero ones.

For our example we therefore choose to place horseshoe priors independently on
each $K^{*}_{ij}$, which shrink small values towards zero but leave larger
values untouched (Carvalho et al.,, 2009):

$\displaystyle\xi_{ij}$ $\displaystyle\sim\text{Cauchy}(0,1)$ (8)
$\displaystyle K^{*}_{ij}$ $\displaystyle\sim\mathcal{N}(0,\xi_{ij})$
$\displaystyle\text{ for }i,j=1\dots M$ (stable only) (9)

##  6 Posterior Inference

To summarise the above, we fit the following multivariate Hawkes process in a
Bayesian manner. The intensity in dimension $i$ is

$\lambda_{i}(t)=\phi\left(\mu_{i}+\sum_{j=1}^{M}\sum_{l:t>t_{j\,l}}\left\\{g(\mathbf{K}^{*})\right\\}_{ji}\,g_{ji}(t-t_{j\,l})\right)$
(10)

with link function

$\displaystyle\phi(x)=\max(0,x)$ (11)

The likelihood is

$\displaystyle
p(Y)=\prod_{m=1}^{M}\prod_{i=1}^{N_{m}}\lambda_{m}(t_{i}|\theta)\,exp\left(-\Lambda\right)$
(12)

where $\Lambda=\sum_{m=1}^{M}\int_{0}^{T_{max}}\lambda_{m}(x)\,dx$. The
following priors are chosen for the background rate:

$\displaystyle\mu_{i}$ $\displaystyle\sim\mathcal{U}(0,10)$
$\displaystyle\text{ for }i=1\dots M$ (13)

We found that in practice prior knowledge is often available to find sensible
upper bounds for these uniform priors. For example the magnitude of number of
events can help guide the priors for $\mu_{i}$.

For the influence kernels we utilise the popular exponential kernel
$g_{ij}(x)=\beta_{ij}\,exp\left(-\beta_{ij}\,x\right)$ (e.g. Blundell et al.,,
2012; Shelton et al.,, 2018). Here, we assume that all
$\beta_{ii}=\beta_{diag}$ and $\beta_{ij}=\beta_{off}$ when $i\neq j$. Hence,
the values for are the same for all self-influences, as well as cross-
influences, respectively. For $\beta_{diag}$ and $\beta_{off}$ we choose the
following priors.

$\displaystyle\beta_{diag}$ $\displaystyle\sim\mathcal{U}(0,3)$ (14)
$\displaystyle\beta_{off}$ $\displaystyle\sim\mathcal{U}(0,3)$ (15)

In line with the arguments from Section 5.1 we place priors on the entries of
$\mathbf{K}^{*}$. The estimation is carried out both using Normal and sparsity
inducing priors (using the horseshoe prior from Section 5.3). First, we use
Normal priors:

$\displaystyle K^{*}_{ij}$ $\displaystyle\sim\mathcal{N}(0,1)$
$\displaystyle\text{ for }i,j=1\dots M$ (stable only) (16)

In addition, we also fit the model with the horseshoe prior (Carvalho et al.,,
2009):

$\displaystyle\xi_{ij}$ $\displaystyle\sim\text{Cauchy}(0,1)$ (17)
$\displaystyle K^{*}_{ij}$ $\displaystyle\sim\mathcal{N}(0,\xi_{ij})$
$\displaystyle\text{ for }i,j=1\dots M$ (stable only) (18)

We use Stan (Stan Development Team,, 2019) to obtain posterior samples from
this probabilistic model. As each evaluation of the likelihood requires
$\Lambda$, the integration of the likelihood, we implement a cubic Simpson’s
rule to approximate the integral of the intensity function between each event
(see Section 4.2 and Appendix B). To analyse the posterior we generate density
plots and summary statistics.

##  7 Example

This section showcases one example of a multivariate Hawkes process with
excitation and inhibition using the set up from Section 6. We simulate a data
set for $M=3$ dimensions using the following parameter values for
$\Theta=(\mu,\mathbf{K}^{*},\beta_{diag},\beta_{off})$.

$\displaystyle\mu$ $\displaystyle=(0.15,0.15,0.15)$
$\displaystyle\mathbf{K}^{*}$
$\displaystyle=\left(\begin{array}[]{ccc}0.3&-0.3&0\\\ 0&0.3&0.3\\\
0&-0.3&0\end{array}\right)$ $\displaystyle\beta_{diag}$
$\displaystyle=\beta_{off}=0.5$

Data is generated for the interval $[0,T_{max}=1,500]$, resulting in $714$
observations, ($275$, $162$, and $277$ for each dimension). Figure 2 displays
the intensity functions for each dimension under the generating parameters
alongside the events themselves.

Figure 2: Intensity functions for each dimension under the true parameters in
the first three panels and observations on the bottom panel.

The estimation is carried out both using normal priors and sparsity-inducing
horseshoe prior from Section 5.3. Figure 3 displays the posterior estimates
based on $750$ samples obtained by Stan (Stan Development Team,, 2019). Table
1 summarises the posterior means (and standard deviations) for both models.
Most posterior distributions are centred around the parameter used for
sampling, only $\mathbf{K}^{*}_{22}$ is not covered well. Generally, both
models are successful in recovering the parameters, given that it is a rather
challenging task to estimate $14$ heavily intertwined parameters with just
$714$ observations.

The horseshoe prior gives similar results as the Normal priors. As expected,
the horseshoe posteriors have less variance when the true value is actually
zero. For example the posterior mean in the Normal case for $K^{*}_{31}=0$ is
$-0.03$ ($0.07$), whereas the horseshoe posterior mean is $-0.01$ ($0.06$).
This can help to identify truly zero entries in the $\mathbf{K}^{*}$ matrix.

Figure 3: Posterior estimates for all parameters. Black lines are for Normal priors, grey for horseshoe priors. Vertical line indicates the true value. Table 1:  Posterior mean (and standard deviation) rounded to two decimal places for each parameter for two models: using Normal priors and a shrinkage approach (using the horseshoe prior). | true | Normal | horseshoe
---|---|---|---
$\mu_{1}$ | 0.15 | 0.14 (0.02) | 0.14 (0.02)
$\mu_{2}$ | 0.15 | 0.12 (0.02) | 0.12 (0.02)
$\mu_{2}$ | 0.15 | 0.15 (0.02) | 0.16 (0.02)
$K^{*}_{11}$ | 0.30 | 0.36 (0.12) | 0.33 (0.12)
$K^{*}_{21}$ | 0.00 | 0.07 (0.12) | 0.05 (0.09)
$K^{*}_{31}$ | 0.00 | -0.03 (0.07) | -0.01 (0.06)
$K^{*}_{12}$ | -0.30 | -0.39 (0.11) | -0.34 (0.10)
$K^{*}_{22}$ | 0.30 | 0.77 (0.24) | 0.70 (0.22)
$K^{*}_{32}$ | -0.30 | -0.31 (0.09) | -0.28 (0.08)
$K^{*}_{13}$ | 0.00 | -0.03 (0.06) | -0.02 (0.05)
$K^{*}_{23}$ | 0.30 | 0.33 (0.11) | 0.29 (0.11)
$K^{*}_{33}$ | 0.00 | 0.00 (0.07) | 0.00 (0.06)
$\beta_{diag}$ | 0.50 | 0.33 (0.08) | 0.37 (0.09)
$\beta_{off}$ | 0.50 | 0.50 (0.11) | 0.51 (0.11)

##  8 Discussion

In this paper we focused on the multivariate Hawkes process under both
excitation and inhibition. We discussed methodological aspects of the model
and introduced a new, stronger stability criterion. We also showed that the
branching structure can be retained even under inhibition. We examined prior
selection in a Bayesian framework using a reparametrisation based on the total
number of events. The resulting Normal and sparsity priors where then
showcased on an example on simulated data.

Our work can be extended by considering different link functions and a variety
of influence kernels, or differently structured inhibition altogether (for
example along the lines of Apostolopoulou et al.,, 2019). Moreover, all
stability conditions, including our newly proposed LABEL:c3new, are only
sufficient to check stability. It would be of great interest to develop a
criterion that was both necessary and sufficient.

Future work will be carried out on real-world data for one of the many
attractive applications of Hawkes processes with inhibition. Moreover, we
would like to extend the applications to higher dimensions and to incorporate
additional covariate information.

## References

  * Apostolopoulou et al., (2019) Apostolopoulou, I., Linderman, S., Miller, K., and Dubrawski, A. (2019).  Mutually regressive point processes.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
  * Bacry et al., (2020) Bacry, E., Bompaire, M., Gaïffas, S., and Muzy, J.-F. (2020).  Sparse and low-rank multivariate hawkes processes.  Journal of Machine Learning Research, 21(50):1–32.
  * Bacry et al., (2016) Bacry, E., Jaisson, T., and Muzy, J. (2016).  Estimation of slowly decreasing Hawkes kernels: application to high-frequency order book dynamics.  Quantitative Finance, 16(8):1179–1201.
  * Bacry and Muzy, (2016) Bacry, E. and Muzy, J.-F. (2016).  First-and second-order statistics characterization of hawkes processes and non-parametric estimation.  IEEE Transactions on Information Theory, 62(4):2184–2202.
  * Blundell et al., (2012) Blundell, C., Beck, J., and Heller, K. A. (2012).  Modelling reciprocating relationships with Hawkes processes.  Advances in Neural Information Processing Systems, 25:2600–2608.
  * Bremaud and Massoulie, (1996) Bremaud, P. and Massoulie, L. (1996).  Stability of nonlinear Hawkes processes.  The Annals of Probability, 24(3):1563–1588.
  * Browning et al., (2021) Browning, R., Sulem, D., Mengersen, K., Rivoirard, V., and Rousseau, J. (2021).  Simple discrete-time self-exciting models can describe complex dynamic processes: A case study of COVID-19.  PLOS ONE, 16(4):e0250015.
  * Carvalho et al., (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In Artificial Intelligence and Statistics, pages 73–80. PMLR.
  * Chen and Stindl, (2018) Chen, F. and Stindl, T. (2018).  Direct likelihood evaluation for the renewal hawkes process.  Journal of Computational and Graphical Statistics, 27(1):119–131.
  * Costa et al., (2020) Costa, M., Graham, C., Marsalle, L., and Tran, V.-C. (2020).  Renewal in Hawkes processes with self-excitation and inhibition.  Advances in Applied Probability, 52(3):879–915.
  * Daley and Vere-Jones, (2003) Daley, D. J. and Vere-Jones, D. (2003).  An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods.  Probability and Its Applications, An Introduction to the Theory of Point Processes. Springer-Verlag, New York, NY, 2 edition.
  * Eichler, (2013) Eichler, M. (2013).  Causal inference with multiple time series: principles and problems.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1997):20110613.
  * Eichler et al., (2017) Eichler, M., Dahlhaus, R., and Dueck, J. (2017).  Graphical modeling for multivariate Hawkes processes with nonparametric link functions.  Journal of Time Series Analysis, 38(2):225–242.
  * Ertekin et al., (2015) Ertekin, Ş., Rudin, C., and McCormick, T. H. (2015).  Reactive point processes: A new approach to predicting power failures in underground electrical systems.  The Annals of Applied Statistics, 9(1):122–144.
  * Gelfand, (1941) Gelfand, I. (1941).  Normierte ringe.  Rech. Math. [Mat. Sbornik], 9(1):3–24.
  * Guo et al., (2018) Guo, X., Hu, A., Xu, R., and Zhang, J. (2018).  Consistency and computation of regularized MLEs for multivariate Hawkes processes.  arXiv:1810.02955 [math, stat].
  * Hawkes, (1971) Hawkes, A. G. (1971).  Spectra of some self-exciting and mutually exciting point processes.  Biometrika, 58(1):83–90.
  * Hawkes and Oakes, (1974) Hawkes, A. G. and Oakes, D. (1974).  A cluster process representation of a self-exciting process.  Journal of Applied Probability, 11(3):493–503.
  * Helmstetter and Sornette, (2002) Helmstetter, A. and Sornette, D. (2002).  Subcritical and supercritical regimes in epidemic models of earthquake aftershocks.  Journal of Geophysical Research: Solid Earth, 107(B10):ESE 10–1–ESE 10–21.
  * Jovanović et al., (2015) Jovanović, S., Hertz, J., and Rotter, S. (2015).  Cumulants of Hawkes point processes.  Physical Review E, 91(4):042802.
  * Kalair et al., (2021) Kalair, K., Connaughton, C., and Alaimo Di Loro, P. (2021).  A non-parametric Hawkes process model of primary and secondary accidents on a UK smart motorway.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 70(1):80–97.
  * Kolev and Ross, (2019) Kolev, A. A. and Ross, G. J. (2019).  Inference for etas models with non-poissonian mainshock arrival times.  Statistics and Computing, 29(5):915–931.
  * Lai et al., (2016) Lai, E. L., Moyer, D., Yuan, B., Fox, E., Hunter, B., Bertozzi, A. L., and Brantingham, P. J. (2016).  Topic time series analysis of microblogs.  IMA Journal of Applied Mathematics, 81(3):409–431.
  * Lemonnier and Vayatis, (2014) Lemonnier, R. and Vayatis, N. (2014).  Nonparametric markovian learning of triggering kernels for mutually exciting and mutually inhibiting multivariate hawkes processes.  In Proceedings of the 2014th European Conference on Machine Learning and Knowledge Discovery in Databases - Volume Part II, page 161–176, Berlin, Heidelberg. Springer-Verlag.
  * Lu and Abergel, (2018) Lu, X. and Abergel, F. (2018).  High-dimensional hawkes processes for limit order books: modelling, empirical analysis and numerical calibration.  Quantitative Finance, 18(2):249–264.
  * Mei and Eisner, (2017) Mei, H. and Eisner, J. (2017).  The neural Hawkes process: A neurally self-modulating multivariate point process.  In Advances in Neural Information Processing Systems, volume 30, page 6757–6767. Curran Associates Inc.
  * Miscouridou et al., (2018) Miscouridou, X., Caron, F., and Teh, Y. W. (2018).  Modelling sparsity, heterogeneity, reciprocity and community structure in temporal interaction data.  In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc.
  * Mohler, (2013) Mohler, G. (2013).  Modeling and estimation of multi-source clustering in crime and security data.  The Annals of Applied Statistics, 7(3):1525–1539.
  * Ogata, (1981) Ogata, Y. (1981).  On Lewis’ simulation method for point processes.  IEEE Transactions on Information Theory, 27(1):23–31.
  * Ogata, (1988) Ogata, Y. (1988).  Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association, 83(401):9–27.
  * Rambaldi et al., (2017) Rambaldi, M., Bacry, E., and Lillo, F. (2017).  The role of volume in order book dynamics: a multivariate Hawkes process analysis.  Quantitative Finance, 17(7):999–1020.
  * Rasmussen, (2013) Rasmussen, J. G. (2013).  Bayesian inference for Hawkes processes.  Methodology and Computing in Applied Probability, 15(3):623–642.
  * Serre, (2002) Serre, D. (2002).  Matrices: Theory and Applications.  Springer Science & Business Media, New York, NY, 2 edition.
  * Shelton et al., (2018) Shelton, C., Qin, Z., and Shetty, C. (2018).  Hawkes process inference with missing data.  Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence.
  * Stan Development Team, (2019) Stan Development Team (2019).  RStan: the R interface to Stan.  R package version 2.19.2.
  * Sulem et al., (2021) Sulem, D., Rivoirard, V., and Rousseau, J. (2021).  Bayesian estimation of nonlinear Hawkes process.  arXiv:2103.17164 [math, stat].
  * Tetereva, (2018) Tetereva, A. (2018).  Do financial companies communicate to one another in the news? (application of multivariate hawkes graphs to uncover granger causality of financial news).  Econometric Modeling: Capital Markets - Asset Pricing eJournal.
  * Tucker et al., (2019) Tucker, D. J., Shand, L., and Lewis, J. R. (2019).  Handling missing data in self-exciting point process models.  Spatial Statistics, 29:160–176.
  * Veen and Schoenberg, (2008) Veen, A. and Schoenberg, F. P. (2008).  Estimation of space-time branching process models in seismology using an em-type algorithm.  Journal of the American Statistical Association, 103(482):614–624.
  * Xu et al., (2016) Xu, H., Farajtabar, M., and Zha, H. (2016).  Learning Granger causality for Hawkes processes.  In International Conference on Machine Learning, pages 1717–1726. PMLR.
  * Yu et al., (2020) Yu, X., Shanmugam, K., Bhattacharjya, D., Gao, T., Subramanian, D., and Xue, L. (2020).  Hawkesian Graphical Event Models.  In International Conference on Probabilistic Graphical Models, pages 569–580. PMLR.

## Appendices

##  Appendix A Roots of the Intensity Function

This section shows how the roots of the intensity function can be found when
using the exponential kernel. As described in Section 4.2, these can be used
to calculate the integral of the intensity exactly.

###  A.1 Exponential Kernel with All $\beta_{ij}$ Equal

We start with a special case where all $\beta_{ij}$ are equal and examine the
roots both in one dimension and $M$ dimensions.

####  A.1.1 One Dimension

Assume that we have data $Y=(t_{1}\dots t_{N})$ from one dimension and the
intensity function is defined with the exponential kernel:

$\lambda(t)=\mu+\sum_{i:t>t_{i}}K\,\beta\,exp\left(-\beta(t-t_{i})\right)$
(19)

The intensity can only drop below zero when an event happens. Therefore it is
sufficient to check only intervals after an event at which the intensity is
negative to see if the intensity becomes positive again before the next event
happens.

Then we can find the root $t$ between observation $t_{n}$ and $t_{n+1}$ in the
following way:

$\displaystyle\mu+\sum_{i=1}^{n}K\,\beta\,exp\left(-\beta(t-t_{i})\right)=0$
$\displaystyle\mu+\sum_{i=1}^{n}K\,\beta\,exp\left(-\beta
t\right)\,exp\left(\beta t_{i}\right)=0$ $\displaystyle exp\left(-\beta
t\right)\left[\sum_{i=1}^{n}K\,\beta exp\left(\beta t_{i}\right)\right]=-\mu$
$\displaystyle t=\frac{log\left(\frac{-\mu}{\sum_{i=1}^{n}K\,\beta
exp\left(\beta t_{i}\right)}\right)}{-\beta}$

if $t_{n}<t<t_{n+1}$.

####  A.1.2 $M$ Dimensions

Now assume we have $M$ dimensional data $Y_{1}=\left(t_{1\,1}\dots
t_{1\,N_{1}}\right)\dots Y_{M}=\left(t_{M\,1}\dots t_{M\,N_{M}}\right)$ with
intensity function

$\lambda_{i}(t)=\mu_{i}+\sum_{j=1}^{M}\sum_{q:t>t_{q\,l}}K_{ji}\,g_{ji}(t-t_{j\,l})$
(20)

Again, we only intervals after an event at which the intensity would be
negative to see if the intensity becomes positive again before the next event
happens. However, we need to check in each interval in each dimension.

To check in dimension $m$ after an observation $t_{n}$ from arbitrary
dimension we can start with

$\displaystyle\mu_{m}+\sum_{q=1}^{M}\sum_{i:t_{qi}<t_{n}}K_{q}\,\beta\,exp\left(-\beta(t-t_{qi})\right)=0$

which leads to the following expression for the root

$\displaystyle
t=\frac{log\left(\frac{-\mu_{m}}{\sum_{q=1}^{M}\sum_{i:t_{q\,i}<t_{n}}K_{q}\,\beta
exp\left(\beta t_{q\,i}\right)}\right)}{-\beta}$

if $t_{n}<t<t_{n+1}$.

###  A.2 General Exponential Kernel

Now we examine the general case when all betas are different. Here, finding
the roots is a polynomial problem:

$\displaystyle\mu_{m}+\sum_{q=1}^{M}\sum_{i:t_{qi}<t_{n}}K_{q}\,\beta_{q}\,exp\left(-\beta_{q}(t-t_{qi})\right)$
$\displaystyle=0$
$\displaystyle\mu_{m}+\sum_{q=1}^{M}\sum_{i:t_{qi}<t_{n}}K_{q}\,\beta_{q}\,u^{-\beta_{q}}\,exp\left(\beta_{q}t_{qi}\right)$
$\displaystyle=0\text{ where }x=exp(t)$
$\displaystyle\mu_{m}+\sum_{q=1}^{M}x^{-\beta_{q}}\underbrace{\sum_{i:t_{qi}<t_{n}}K_{q}\,\beta_{q}\,exp\left(\beta_{q}t_{qi}\right)}_{v_{q}}$
$\displaystyle=0$ $\displaystyle\mu_{m}+\sum_{q=1}^{M}x^{-\beta_{q}}v_{q}$
$\displaystyle=0$

This is now a polynomial in $x$ that needs to be solved.

##  Appendix B Approximating $\Lambda$ using the Simpson’s Rule

This section describes the approximation procedure for
$\Lambda=\sum_{m=1}^{M}\int_{0}^{T_{max}}\lambda_{m}(x)\,dx$ as mentioned in
Section 4.2. We use the Cubic Simpson’s Rule where

$\int_{a}^{b}f(x)\,dx\approx\frac{(b-a)}{8}\left[f(a)+3f\left(\frac{2a+b}{3}\right)+3f\left(\frac{a+2b}{3}\right)+f(b)\right]$
(21)

Assume we have $M$ dimensional data $Y_{1}=\left(t_{1\,1}\dots
t_{1\,N_{1}}\right)\dots Y_{M}=\left(t_{M\,1}\dots t_{M\,N_{M}}\right)$ and
intensity function $\lambda_{m}(\cdot)$ for each dimension $m=1\dots M$.
Algorithm 1 outlines the approximation procedure used.

Algorithm 1 Approximating $\Lambda$ using the Simpson’s Rule

1:Set $Y_{total}=(0,t_{1\,1}\dots t_{1\,N_{1}}\dots t_{M\,1}\dots
t_{M\,N_{M}},T_{max})$

2:Order all entries of $Y_{total}$ and call the result $X=(x_{1}<\dots<x_{P})$
where $P=2+\sum_{m}^{M}N_{m}$

3:Define $X^{t}=\\{x_{i}:x_{i}\leq t\\}$

4:Set $res=0$

5:for $i$ in $1:(P-1)$ do

6: Set $a=x_{i}$

7: Set $b=x_{i+1}$

8: for $m$ in $1:M$ do

9: Set res = res +
$\frac{(b-a)}{8}\left[\lambda_{m}(a|X^{a})+3\lambda_{m}\left(\frac{2a+b}{3}|X^{a}\right)+3\lambda_{m}\left(\frac{a+2b}{3}|X^{a}\right)+\lambda_{m}(b|X^{a})\right]$

10: end for

11:end for

12:Return $res$ as the approximation of $\Lambda$

##  Appendix C Branching Structure with Inhibition

This section outlines a sampling procedure that preserves the branching
structure (see Section 2.3) of a Hawkes process, even if inhibition is
present.

We define the direct parent of an event as follows. If $t_{j}$ is a direct
offspring of $t_{i}$, then $t_{i}$ is the parents of $t_{j}$. For each event
$t_{i}$ we define $B(t_{i})=0$ if $t_{i}$ was an immigrant event. Otherwise
$B(t_{i})$ takes the index the direct parent of $t_{i}$. For example if
$t_{i}$ has one offspring $t_{j}$ then $B(t_{j})=i$. Moreover, let $C(t_{i})$
contain all direct and indirect (‘total’) offsprings of $t_{i}$.

We assume that $\phi(x)=\max(0,x)$, as it ensures that
$\lambda(t\,|\,\mathbf{K}^{+})>\lambda(t\,|\,\mathbf{K})$ for all $t$.

First, we review the classic thinning approach.

######  Lemma 7.

Let $X=(x_{1}\dots x_{P})$ be event times from a point process with intensity
$\bar{\lambda}(\cdot)$. Suppose there is another point process with intensity
function $\lambda(\cdot)$ such that $\bar{\lambda}(x)\geq\lambda(x)$ for all
$x\geq 0$. To obtain samples from $\lambda(\cdot)$, take $X=(x_{1}\dots
x_{P})$ as candidates. Let $Y$ be the empty set. Each $x_{i}$, $i=1\dots P$,
is added to $Y$ with probability $\lambda(x_{i})/\bar{\lambda}(x_{i})$. The
event times in $Y$ are samples from the process with intensity
$\lambda(\cdot)$.

This is the popular ‘thinning’ sampler for point processes, as introduced by
Ogata, (1981). Algorithm 2 outlines this approach, where each event is thinned
independently, regardless of the branching structure. In particular, in Line 6
we calculate $b_{original}$ using all generated events in $Y_{original}$,
regardless how many of those are still in $Y$.

Now we look at a different sampling procedure as outlines in Algorithm 3. It
differs from the algorithm above in two aspects. Here, $b$ is calculated with
only the retained events (Line 8), and offsprings are removed with its parents
(Line 10). Note that our sampling approach under inhibition is inherently
different to the one introduced by Apostolopoulou et al., (2019) as they work
with a different model set up that specifies inhibition through a
multiplicative factor.

######  Theorem 8.

Algorithm 3 generates samples from the correct target intensity
$\lambda(\cdot\,|\,\mathbf{K})$.

###### Proof.

Before any event is removed, the algorithm mirrors the classic thinning
approach.

Let us now assume that event $t_{i}$ is removed (Line 10), along with all of
its offsprings, of which $k$ lie before event $t_{j}$, which remains.

$\displaystyle Y_{original}$ $\displaystyle=(t_{1}\dots
t_{i-1},t_{i},t_{i+1}\dots t_{i+k},t_{j},\dots)$ $\displaystyle Y$
$\displaystyle=(t_{1}\dots t_{i-1},t_{j},\dots)$

Hence, at this point we could discard all events after $t_{i-1}$ and instead
sample from any proposal (that is larger than the target) and thin again
according to the classic thinning approach. We choose to sample this new
proposal event from $\lambda(\cdot\,|\,\mathbf{K^{+}},\mathcal{H}=Y)$ using
the branching structure. However, we have already generated such an event,
namely $t_{j}$. Since $t_{j}$ has not been removed, it was either an immigrant
event or an offspring of an event in $Y$, therefore it can be viewed as
sampled from the desired $\lambda(\cdot\,|\,\mathbf{K^{+}},\mathcal{H}=Y)$.
Hence, the ration in Line 9 follows the usual thinning procedure and Algorithm
3 describes a sampling procedure that samples from the target intensity
function $\lambda(\cdot\,|\,\mathbf{K})$. ∎

Crucially, all parents of events in $Y$ are also in $Y$ by construction. Hence
the branching structure is retained, even under inhibition.

Algorithm 2 Classic Thinning Sampling without branching structure

1:Sample $Y_{original}=(t_{1}<\dots<t_{N})$ from
$\lambda(\cdot\,|\,\mathbf{K}^{+})$

2:Set $Y=Y_{original}$

3:for $i$ in $1:N$ do

4: sample $u\sim\mathcal{U}[0,1]$

5: set $a=\lambda(t_{i}\,|\,\mathbf{K},\mathcal{H}=Y)$

6: set
$b_{original}=\lambda(t_{i}\,|\,\mathbf{K^{+}},\mathcal{H}=Y_{original})$

7: if $u>a/b_{original}$ then

8: remove $t_{i}$ from $Y$

9: end if

10:end for

11:Return $Y$ as a sample from the multivariate Hawkes Process

Algorithm 3 Sampling with inhibition while retaining the branching structure

1:Sample $Y_{original}=(t_{1}<\dots<t_{N})$ from
$\lambda(\cdot\,|\,\mathbf{K}^{+})$ using branching

2:Keep the $B(t_{i})$ for all $i=1\dots N$

3:Set $Y=Y_{original}$

4:for $i$ in $1:N$ do

5: if $t_{i}\not\in Y$ then next

6: else sample $u\sim\mathcal{U}[0,1]$

7: set $a=\lambda(t_{i}\,|\,\mathbf{K},\mathcal{H}=Y)$

8: set $b=\lambda(t_{i}\,|\,\mathbf{K^{+}},\mathcal{H}=Y)$

9: if $u>a/b$ then

10: remove $t_{i}$ and $C(t_{i})$ from $Y$

11: end if

12: end if

13:end for

14:Return $Y$ as a sample from the multivariate Hawkes Process

##  Appendix D Total Number of Offsprings

This section contains the proof of the total number of offsprings from Section
5.1.

We write $K_{ij}^{*}$ is the total number of offsprings in dimension $j$ that
are produced by an event cascade started by an immigrant event in $i$. This
can be interpreted as the marginal influence from $i$ onto $j$ as both direct
and indirect influences are taken into account. We write $\mathbf{K}^{*}$
where $K_{ij}^{*}$ is the entry in row $i$ and column $j$.

Here, $N=(N_{1}\dots N_{M})^{T}=$ is the total number of events in each
dimension, which can be calculated by $N=(I-\mathbf{K}^{T})^{-1}\mathbf{\mu}$,
as described by Hawkes, (1971); Jovanović et al., (2015).

######  Lemma 9.

The number of events in each dimension is

$N=(N_{1}\dots N_{M})^{T}=\left(\mathbf{K}^{*}+I\right)^{T}\mathbf{\mu}$

where $I$ is the identity matrix of appropriate dimension $M$.

###### Proof.

The Hawkes process can be written as the superposition of Poisson processes as
outlined by the branching structure interpretation (see Section 2.3). Hence,
each $N_{i}$ consist of two parts:

  1. 1.

The number of immigrants events from the background process with rate
$\mu_{i}$

  2. 2.

The number of offsprings in $i$ from an immigrant event in each dimension $j$

Therefore, we can write

$N_{i}=\mu_{i}T_{max}+\sum_{j=1}^{M}\mu_{j}T_{max}K_{ji}^{*}$ (22)

In matrix notation this is

$N=(N_{1}\dots N_{M})^{T}=\left(\mathbf{K}^{*}+I\right)^{T}\mathbf{\mu}$ (23)

where $I$ is the identity matrix of appropriate dimension $M$. ∎

######  Theorem 10.

The total number of offsprings is

$\mathbf{K}^{*}=\left(I-\mathbf{K}\right)^{-1}-I$

###### Proof.

We use the definition of $N$ from the literature

$N=(I-\mathbf{K}^{T})^{-1}\mathbf{\mu}$ (24)

and relate it to Lemma 9:

$(I-\mathbf{K^{T}})^{-1}\mathbf{\mu}=\left(\mathbf{K}^{*}+I\right)^{T}\mathbf{\mu}$
(25)

Rearranging, this allows us to write $\mathbf{K}^{*}$ as

$\displaystyle\mathbf{K}^{*}$ $\displaystyle=\left(I-\mathbf{K}\right)^{-1}-I$
(26)

∎


# Quantifying Noise of Dynamic Vision Sensor

Evgeny V. Votyakov  DeepCamera MRG, CYENS Centre of Excellence, Nicosia,
Cyprus  Alessandro Artusi  DeepCamera MRG, CYENS Centre of Excellence,
Nicosia, Cyprus

###### Abstract

Dynamic visual sensors (DVS) are characterized by a large amount of background
activity (BA) noise, which it is mixed with the original (cleaned) sensor
signal. The dynamic nature of the signal and the absence in practical
application of the ground truth, it clearly makes difficult to distinguish
between noise and the cleaned sensor signals using standard image processing
techniques. In this letter, a new technique is presented to characterise BA
noise derived from the Detrended Fluctuation Analysis (DFA). The proposed
technique can be used to address an existing DVS issues, which is how to
quantitatively characterised noise and signal without ground truth, and how to
derive an optimal denoising filter parameters. The solution of the latter
problem is demonstrated for the popular real moving-car dataset.

###### Index Terms:

event cameras, background activity noise, optimal filter parameters, detrended
fluctuation analysis.

11(2.5, 0.4)

This work has been submitted to the IEEE Signal Processing Letters. Copyright
may be transferred without notice, after which this version may no longer be
accessible. ©IEEE

##  I Introduction

Dynamic vision sensor (DVS), also known as silicon retina, is a
microelectronic chip for events camera detecting changes in brightness or
intensity on a per-pixel basis and storing these changes asynchronously as
events [1]. One of the major issues with this type of sensors, is the fact
that large amount of background activity (BA) noise can be present in the
clean signal, which it may limit the usage of the signal in practical
applications. Properly setting of the sensor may partially solve the issue by
means of improving the cleanness of the received events stream and this has
been demonstrated in various previous works [2], [3], [4], [5], [6]. However,
large areas of the signal may still be affected by the BA noise, due to the
difficulties to recognise and remove noise induced by the background activity
(BA) of the DVS matrix.

A key principle to distinguish BA noise from a signal has been firstly
introduced by Delbruck [7],[8]. It is based on the observation that a random
noise event is spatially and temporally isolated. The noise is identified by
checking eight neighboring pixels and using a single $\Delta T$ parameter for
the time interval between the checked event and its correlated events. Also a
variety of modifications of the original approach have been presented [9]: (1)
accelerate search of neighboring pixels in the vicinity of the central event
without change of the original principle (k-noise [10] and DWF [11] filters);
(2) introduce additional parameters $\Delta S_{x},\Delta S_{y}$ to expand the
allowed space vicinity as well as the minimal number of captured events
$\rho_{\text{min}}$ in a 3D $xyt$ box as a threshold (y-noise filter [12]);
(3) make projections of the events, which are captured in the 3D box $\Delta
S_{x},\Delta S_{y},\Delta T$ onto central plane of the box with various
weights (TS [13], IETS [14], MLPF [11] filter). Also, deep-learning based
approaches have been presented: PUGM [15], EDnCNN [16], AEDNet [17], and
partly MLPF [11].

It is commonly known that it is hard to assess quality of denoised DVS data
for two reasons. First, the evaluation may be influenced by the so called
accumulation time, i.e., too small only few events are captured, too big the
clean signal may result blurred and so it will be identified as noise. Second,
any BA filter has tuning parameters, and they need to be carefully selected
and based on the adopted accumulation time same parameters values may provide
different results.

To avoid an evaluation based mainly on subjective taste, one needs an
objective criterion. We may have two approaches, using a ground truth or using
image recognition applications performances to evaluate the quality of the
noise filtering output. The former approach is available only for synthetic
datasets with an artificial noise injected [11], and the latter one exploits a
concept of warped events [18], where one needs (a) to perform a clustering of
events streams, (b) to determine a optical flow for each cluster by means of
optimisation of a contrast maximization goal function, and (c) to use the
results of the optimisation to compare different denoising approaches.
Obviously, the latter is much more complex. To this scope TS [13] and IETS
[14] filters were proposed.

Attempts to label each event as noise or signal have been proposed[16], [17],
however, due to the heavy data load of DVS it is completely impractical
approach [19].

The above motivations are clearly showing how traditional approaches are
limiting themselves in working with each event by checking if it is noise or
signal. In this letter we introduce a novel objective criterion, to evaluate
noise filtering performances, where a single parameter is derived from the
statistical properties of the DVS time series. We quantify filtered DVS noise
with detrended fluctuation analysis (DFA). After a brief introduction of the
DFA principles, we are showing how DFA might be adapted for DVS data, and then
demonstrating its efficiency by separating noise from the cleaned signal. This
is shown by testing it on a popular moving-car dataset to find optimal filter
characteristics. To our best knowledge, this is the first attempt to use DFA
approach to DVS data.

##  II Dynamic Visual Sensor (DVS)

Dynamic Visual Sensor (DVS) [1] is a 2D matrix of asynchronously pixels. Each
pixel has static and dynamic variables. The static one is for previous
reference illumination of the pixel and the dynamic one is for the
instantaneous illumination. At every tact, the instantaneous and reference
illuminations are compared. If their difference is more than a threshold, the
pixel updates its reference and triggers an event $e=e(t,x,y,p)$, being $t$ an
instance of time, when the event is triggered, $(x,y)$ pixel’s 2D coordinates,
and $p=\pm 1$ is polarity of the triggered event (positive, if illumination
increases, and negative otherwise). As a result, DVS data is a stream of
events $\\{e_{1},e_{2},\dots,e_{N}\\}$.

##  III Background Activity (BA) filter

A simple BA filter works as follows [7], [8], [10], [20], [11]. Let
$\sigma_{ij}$ be a Boolean variable for a correlation between $e_{i}$ and
$e_{j}$ events in the 3D box $(\Delta S_{x},\Delta S_{y},\Delta T)$:

$\displaystyle\sigma_{ij}$ $\displaystyle=|x_{i}-x_{j}|\leq\Delta S_{x}$
$\displaystyle\land|y_{i}-y_{j}|\leq\Delta S_{y}\land|t_{i}-t_{j}|\leq\Delta
T.$ (1)

When $e_{j}$ is vicinity of $e_{i}$, then $\sigma_{ij}$ is true. With
$\llbracket\cdot\rrbracket$ as Iverson bracket
($\llbracket\text{true}\rrbracket=1$,$\llbracket\text{false}\rrbracket=0$),
the total number of correlated events $\rho_{i}$ in $e_{i}$’s vicinity is
given by

$\displaystyle\rho_{i}(\Delta S_{x},\Delta S_{y},\Delta T)$
$\displaystyle=\sum_{j\neq i}\llbracket\sigma_{ij}\rrbracket,$ (2)

Then

$\displaystyle e_{i}(t_{i}x_{i},y_{i},p_{i})=$
$\displaystyle\begin{cases}\text{signal},\text{ if
}\rho_{i}\geq\rho_{\text{min}},\\\ \text{noise, otherwise.}\end{cases}$ (3)

We use $\Delta S_{x}=\Delta S_{y}=1$, $\rho_{\text{min}}=1$, and $\Delta T$ is
a single parameter to optimize the BA filter.

##  IV Formulation of the mathematical problem

Take a time series (TS) of DVS events and split it into clean and noise parts
with a BA fliter, Fig. 1a. The clean signal, $(TS)_{clean}$, is composed of
correlated events, i.e. it contains long term correlations as much as
possible. The noise signal, $(TS)_{noise}$, must be as random as possible,
i.e. it contains long term correlations as less as possible. To quantify
correlations, DFA scaling exponents $\alpha$ extracted from corresponding TS
is used, Fig. 1b.

##  V Detrended Fluctuation Analysis

Detrended Fluctuation Analysis (DFA) [21] is a statistical method to discover
long-term correlations in a time series. DFA is numerously exemplified:
anaerobic threshold derived from heart rate variabilities [22], earthquakes
[23], lightning thunderstorms flash sequences in atmosphere science [24],
financial markets [25]. In image processing, Ramirez [26] explored DFA at
describing roughness of images. Technical aspects are given by Hu [27] and
Kantelhardt [28]. For DVS applications, the DFA has not been yet published.

Figure 1: (a) clean $(TS)_{\text{clean}}$ and noise $(TS)_{\text{noise}}$ are
produced by a BA filter (1) from an initial $(TS)_{\text{full}}$. (b) DFA
results. Slope of the line is DFA scaling exponent $\alpha$ characterizing an
impact of long term correlations.

Let $x$ be a random variable, and there is a finite series composed of $x$:

$\displaystyle\\{x_{1},x_{2},\dots,x_{N}\\}.$ (4)

In the case of DVS, there is a time series of events, and as a random variable
$x$ stands a time interval $\delta t_{i}=t_{i}-t_{i-1}$ between adjacent
events.

Let $X$ be a new random variable equal to a cumulative sum, $X_{1}=x_{1}$,
$X_{2}=x_{1}+x_{2}$, $\dots$, $X_{N}=x_{1}+x_{2}+\dots+x_{N}$. Notice that for
a DVS stream, the cumulative sums are exactly timestamps: $X_{i}=t_{i}$.
Compose a cumulative series:

$\displaystyle\\{X_{1},X_{2},\dots,X_{N}\\}.$ (5)

Series (5) is partitioned in local segments of equal length $n$, and an
independent least-square fit (local trend) is performed inside the $k$-th
segment, $1\leq k\leq\lfloor N/n\rfloor$. Let $\hat{X}_{i}$ be the values
predicted by a local trend, then, $(X_{i}-\hat{X}_{i})$ are detrended
fluctuations. Mean local fluctuation of the $k$-th segment is:

$\displaystyle
F_{k}(n)=\left\\{\frac{1}{n}\sum_{l=in+1}^{in+n}\left(X_{l}-\hat{X}_{l}\right)^{2}\right\\}^{1/2},$
(6)

and the mean total fluctuation depending on $n$ is averaged over all mean $k$
local fluctuations:

$\displaystyle F(n)=\left\\{\frac{1}{\lfloor N/n\rfloor}\sum_{k=1}^{\lfloor
N/n\rfloor}F_{k}(n)^{2}\right\\}^{1/2}.$ (7)

The computation of eq. 6 and 7 run for segments of regularly increasing length
$n\in\\{m_{1},m_{2},\dots,m_{M}$}, where integers $m_{k}$ must obey the
condition $m_{k}\simeq qm_{k-1}$, $q>1$ is a constant, i.e. $m_{k}$ must form
a geometric progression, being $m_{1}=4$ and $m_{M}\simeq N$. As a result, a
fluctuation function $F(n)$ is produced. If the series (4) is of long-range
power-law correlations, $F(n)$ increases by a power-law:

$\displaystyle F(n)\simeq Cn^{\alpha}.$ (8)

The scaling exponent $\alpha$ gives a metric on self-correlations inside the
initial series (4). To find $\alpha$, it is convenient to make double $\log
n-\log F(n)$ plot as shown in Fig. 1b: $\alpha\simeq 0.5$ means uncorrelated
data, $\alpha>0.5$ positive long-range correlations, $\alpha<0.5$ \-
anticorrelation, and $\alpha>1$ the data non-stationary [29], [30].

##  VI Results

To demonstrate how to assess objectively the quality of BA denoising without
having a ground truth, we have used as noise input data a traditional dataset
(slot-car dataset), and as filter we have adopted a commonly used basic filter
in the field (Delbruck’s filter [7]). Even do our approach is not sensitive to
hot pixels, to simplify the demonstration of our results without loosing in
generality, hot pixels were removed form the input dataset.

Time bins from 4e6 till 16e6 with starting time 1e6 (all time units are in
$\mu s$) are used. The time bin duration has no effect on the results, except
for the longest (32e6) are better for the analysis (they provide more data for
statistics), and short time bins are better for 3D $xyt$ visualisation (the
long-time bins produce dirty 3D figures).

Fig.(2(a) \- 2(h)) show BA denoising for different $\Delta T$ from 1000 till
16000. All the noise is removed well, but signal is removed too, as one can
see easily in these 3D figures. It is easy to see signal traces in the
filtered noise data because we are dealing with a simple dataset.

(a) $\Delta T=1000$ clean

(b) $\Delta T=1000$ noise

(c) $\Delta T=4000$ clean

(d) $\Delta T=4000$ noise

(e) $\Delta T=8000$ clean

(f) $\Delta T=8000$ noise

(g) $\Delta T=16000$ clean

(h) $\Delta T=16000$ noise

Figure 2:  BAF applied at various $\Delta T$.

Now we show how the signal captured by BA filter and manifested in the noise
time series, disturbs noise statistical properties. Fig. 3 presents DFA
results for the clean and noise time series at various $\Delta T$ shown in
Fig.(2(a) \- 2(h)). The DFA exponent $\alpha$ corresponds to the slope of
crossover lines, see legend in Fig. 3b. At $\Delta T=1000$, $\alpha=0.74$ what
is an evidence that noise is correlated, while $\alpha=0.50$ at $\Delta
T=16000$ proves the filtered BA noise is random. Indeed, by comparing Fig.
2(b) with Fig. 2(h), one can see that in the first sample there are
significant traces of spatially ordered signal in the noise time series, while
in the second case, the signal traces are almost absent.

Figure 3: DFA applied to the noise time series.  Figure 4: SNR and $\alpha$ as
a function of $\Delta T$.

The Signal-Noise-Ratio (SNR) and DFA exponent $\alpha$ as a function of BA
filter parameter $\Delta T$ are summarized in Fig. 4 to draw a conclusion that
the BA filter $\Delta T$ for the slot-car data must be taken as large as
possible since it provides the higher SNR with the $\alpha$ confirming random
noise, and these both are asymptotically converging to their stationary
values. This is trivial, since the highest SNR is of course better, however,
such triviality is valid only for simple datasets, when BA noise is obviously
random, i.e. we can easily distinguish in separate snapshots signal and noise
as in the given example, and there is no doubt that suspected noise is random.
Otherwise, noise long term correlations must be checked.

Now a question arise about how to select the optimal $\Delta T$. At larger
$\Delta T$, BA filter works slower since the filter window captures more
events. Each event in the window must be checked for correlations with the
central event, i.e. a naive BA filter scales $O(N^{2})$ (there are technical
implementations of BA filters free of this constrain, but then there appear
other complexities). Then, we have two natural limits for BAF $\Delta T$: the
lower $\Delta T$ is because of most events are treated as noisy but the filter
is fast, and the higher limit is because of the finite computational
resources. This is true already for simple scenes; for complex scenes one can
not distinguish by eye all the details and separate suspected noise from the
signal. Then, one can check statistical noise properties.

##  VII Conclusion

Detrended functional analysis (DFA) is presented to characterise quality of
denoised DVS data without the availability of the ground truth. It is derived
from statistical properties of filtered random noise extracted from a full DVS
time series. The DFA is adapted to DVS data and tested with a popular slot-car
dataset refined by a simple BA filter. It is demonstrated that DFA might be an
useful tool since scaling exponents capture DVS signal in the filtered BA
noise. As future works we would like to compare more complex BA filters, e.g.,
ynoise or knoise on complex datasets, such as data characterised by
independently moving objects (IMO), where each object produces events with own
statistical properties. Other possible future directions are for example, to
explore how DFA exponents are linked with Poisson distribution of time
intervals, investigate the possible outcome if instead of time intervals taken
as random variable $x=t_{i}-t_{i-1}$, Eq.(4), one can consider the number of
events $x=N$ coming inside some a priori fixed interval $\Delta
t_{\text{ref}}$.

##  VIII Acknowledgment

This project has received funding from the European Union’s Horizon 2020
Research and Innovation Programme under Grant Agreement No 739578 and the
Government of the Republic of Cyprus through the Deputy Ministry of Research,
Innovation and Digital Policy.

## References

  * [1] G. Gallego, T. Delbrück, G. Orchard, C. Bartolozzi, B. Taba, A. Censi, S. Leutenegger, A. J. Davison, J. Conradt, K. Daniilidis, and D. Scaramuzza, “Event-based vision: A survey,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 44, no. 1, pp. 154–180, 2022.
  * [2] G. Gallego, M. Gehrig, and D. Scaramuzza, “Focus is all you need: Loss functions for event-based vision,” _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 12 272–12 281, 2019.
  * [3] M. Muglikar, M. Gehrig, D. Gehrig, and D. Scaramuzza, “How to calibrate your event camera,” _2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , pp. 1403–1409, 2021.
  * [4] S. Shiba, Y. Aoki, and G. Gallego, “Secrets of event-based optical flow,” in _European Conference on Computer Vision_ , 2022.
  * [5] R. Graça, B. Mcreynolds, and T. Delbrück, “Optimal biasing and physical limits of dvs event noise,” _ArXiv_ , vol. abs/2304.04019, 2023.
  * [6] ——, “Shining light on the dvs pixel: A tutorial and discussion about biasing and optimization,” _ArXiv_ , vol. abs/2304.04706, 2023.
  * [7] T. Delbruck, “Frame-free dynamic digital vision,” in _International Symposium on Secure-Life Electronics_ , vol. 1, no. 1. University of Tokyo, March 2008, pp. 21–26, in: Proceedings of International Symposium on Secure-Life Electronics, Advanced Electronics for Quality Life and Society, Univ. of Tokyo, Mar. 6-7, 2008. [Online]. Available: https://doi.org/10.5167/uzh-17620
  * [8] H. Liu, C. Brandli, C. Li, S.-C. Liu, and T. Delbrück, “Design of a spatiotemporal correlation filter for event-based sensors,” _2015 IEEE International Symposium on Circuits and Systems (ISCAS)_ , pp. 722–725, 2015.
  * [9] S. J. Ding, J. Chen, Y. Wang, Y. Kang, W. Song, J. Cheng, and Y. Cao, “E-mlb: Multilevel benchmark for event-based camera denoising,” _ArXiv_ , vol. abs/2303.11997, 2023. [Online]. Available: https://arxiv.org/abs/2303.11997
  * [10] A. Khodamoradi and R. Kastner, “O(n)-space spatiotemporal filter for reducing noise in neuromorphic vision sensors,” _IEEE Transactions on Emerging Topics in Computing_ , 2018.
  * [11] S. Guo and T. Delbruck, “Low cost and latency event camera background activity denoising,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2022.
  * [12] Y. Feng, H. Lv, H. Liu, Y. Zhang, Y. Xiao, and C. Han, “Event density based denoising method for dynamic vision sensor,” _Applied Sciences_ , 2020.
  * [13] X. Lagorce, G. Orchard, F. Galluppi, B. E. Shi, and R. B. Benosman, “Hots: a hierarchy of event-based time-surfaces for pattern recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , pp. 1346–1359, 2016.
  * [14] R. W. Baldwin, M. Almatrafi, J. R. Kaufman, V. Asari, and K. Hirakawa, “Inceptive event time-surfaces for object classification using neuromorphic cameras,” in _Image Analysis and Recognition_ , F. Karray, A. Campilho, and A. Yu, Eds. Cham: Springer International Publishing, 2019, pp. 395–403.
  * [15] J. Wu, C. Ma, L. Li, W. Dong, and G. Shi, “Probabilistic undirected graph based denoising method for dynamic vision sensor,” _IEEE Transactions on Multimedia_ , vol. 23, pp. 1148–1159, 2021.
  * [16] R. W. Baldwin, M. Almatrafi, V. Asari, and K. Hirakawa, “Event probability mask (epm) and event denoising convolutional neural network (edncnn) for neuromorphic cameras,” 2020.
  * [17] H. Fang, J. Wu, L. Li, J. Hou, W. Dong, and G. Shi, “Aednet: Asynchronous event denoising with spatial-temporal correlation among irregular data,” in _Proceedings of the 30th ACM International Conference on Multimedia_ , ser. MM ’22. New York, NY, USA: Association for Computing Machinery, 2022, p. 1427–1435. [Online]. Available: https://doi.org/10.1145/3503161.3548048
  * [18] G. Gallego, H. Rebecq, and D. Scaramuzza, “A unifying contrast maximization framework for event cameras, with applications to motion, depth, and optical flow estimation,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 3867–3876.
  * [19] N. Xu, L. Wang, J. Zhao, and Z. Yao, “Denoising for dynamic vision sensor based on augmented spatiotemporal correlation,” _IEEE Transactions on Circuits and Systems for Video Technology_ , vol. 33, no. 9, pp. 4812–4824, 2023.
  * [20] C. Yang, H. Feng, Z. hai Xu, Q. Li, and Y. ting Chen, “The spatial correlation problem of noise in imaging deblurring and its solution,” _J. Vis. Commun. Image Represent._ , vol. 56, pp. 167–176, 2018.
  * [21] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger, “Mosaic organization of dna nucleotides,” _Phys. Rev. E_ , vol. 49, pp. 1685–1689, Feb 1994. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.49.1685
  * [22] B. Rogers, D. Giles, N. Draper, O. Hoos, and T. Gronwald, “A new detection method defining the aerobic threshold for endurance exercise and training prescription based on fractal correlation properties of heart rate variability,” _Front. Physiol._ , vol. 11, p. 596567, 2020.
  * [23] T. Kataoka, T. Miyaguchi, and T. Akimoto, “Detrended fluctuation analysis of earthquake data,” _Phys. Rev. Res._ , vol. 3, p. 033081, Jul 2021. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevResearch.3.033081
  * [24] X. Gou, M. Chen, and G. Zhang, “Time correlations of lightning flash sequences in thunderstorms revealed by fractal analysis,” _Journal of Geophysical Research: Atmospheres_ , vol. 123, no. 2, pp. 1351–1362, 2018. [Online]. Available: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2017JD027206
  * [25] K. Shrestha, “Multifractal detrended fluctuation analysis of return on bitcoin*,” _International Review of Finance_ , vol. 21, no. 1, pp. 312–323, 2021. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1111/irfi.12256
  * [26] J. Alvarez-Ramirez, E. Rodriguez, I. Cervantes, and J. Carlos Echeverria, “Scaling properties of image textures: A detrending fluctuation analysis approach,” _Physica A: Statistical Mechanics and its Applications_ , vol. 361, no. 2, pp. 677–698, 2006. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378437105007193
  * [27] K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, and H. Eugene Stanley, “Effect of trends on detrended fluctuation analysis,” _Physical Review E_ , vol. 64, no. 1, Jun. 2001. [Online]. Available: http://dx.doi.org/10.1103/PhysRevE.64.011114
  * [28] J. W. Kantelhardt, E. Koscielny-Bunde, H. H. Rego, S. Havlin, and A. Bunde, “Detecting long-range correlations with detrended fluctuation analysis,” _Physica A: Statistical Mechanics and its Applications_ , vol. 295, no. 3, pp. 441–454, 2001. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0378437101001443
  * [29] M. S. Taqqu, V. Teverovsky, and W. Willinger, “Estimators for long-range dependence: An empirical study,” _Fractals_ , vol. 03, no. 04, pp. 785–798, 1995. [Online]. Available: https://doi.org/10.1142/S0218348X95000692
  * [30] M. Höll, K. Kiyono, and H. Kantz, “Theoretical foundation of detrending methods for fluctuation analysis such as detrended fluctuation analysis and detrending moving average,” _Phys. Rev. E_ , vol. 99, p. 033305, Mar 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.99.033305
  * [31] L. Kristoufek, “Detrending moving-average cross-correlation coefficient: Measuring cross-correlations between non-stationary series,” _Physica A: Statistical Mechanics and its Applications_ , vol. 406, pp. 169–175, 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S037843711400209X