∎

11institutetext:  K.Sakuraba, S.Mori 22institutetext: Department of
Mathematics and Physics, Graduate School of Science and Technology, Hirosaki
University
Bunkyo-cho 3, Hirosaki, Aomori 036-8561, Japan
W. Kurebayashi 33institutetext: Department of Mechanical Science and
Engineering, Graduate School of Science and Technology, Hirosaki University
Bunkyo-cho 3, Hirosaki, Aomori 036-8561, Japan
M. Hisakado 44institutetext: Nomura Holdings, Inc., Otemachi 2-2-2, Chiyoda-
ku, Tokyo 100-8130, Japan

# Self-exciting negative binomial distribution process and critical properties
of intensity distribution

Kotaro Sakuraba  Wataru Kurebayashi  Masato Hisakado  Shintaro Mori

(Received: date / Accepted: date)

###### Abstract

We study the continuous time limit of a self-exciting negative binomial
process and discuss the critical properties of its intensity distribution. In
this limit, the process transforms into a marked Hawkes process. The
probability mass function of the marks has a parameter $\omega$, and the
process reduces to a ”pure” Hawkes process in the limit $\omega\to 0$. We
investigate the Lagrange–Charpit equations for the master equations of the
marked Hawkes process in the Laplace representation close to its critical
point and extend the previous findings on the power-law scaling of the
probability density function (PDF) of intensities in the intermediate
asymptotic regime to the case where the memory kernel is the superposition of
an arbitrary finite number of exponentials. We develop an efficient sampling
method for the marked Hawkes process based on the time-rescaling theorem and
verify the power-law exponents.

##  1 Introduction

Recently, high-frequency financial data have become available, and many
studies have been devoted to calibrating models of market micro-
structure(Bacry et al. 2015). Initially, many studies adopted discrete time
models that captured trade dynamics at regular time intervals or through a
sequence of discrete time events, such as trades. Subsequently, the framework
of the continuous time model—the point process— has been gradually applied to
model financial data at the transaction level (Hasbrouck 1991, Engle and
Russell 1998, Bowsher 2007, Hawkes 1971a, b, Hawkes and Oakes 1974, Hawkes
2018, Filimonov and Sornette 2012, 2015, Wheatley et al. 2019).

A point or counting process is characterized by its ”intensity” function,
which represents the conditional probability of an event occurring in the
immediate future. The Hawkes process, introduced by A. G. Hawkes (Hawkes
1971a, b, Hawkes and Oakes 1974), is a type of point process originally
developed to model the occurrence of seismic events. Owing to its simplicity
and flexibility, this model is being increasingly used in high-frequency
finance (Bacry et al. 2015, Kirchner 2017, Blanc et al. 2017, Errais et al.
2010). It can easily capture the interactions between different types of
events, incorporate the influence of intensive factors through marks, and
accommodate non-stationarities.

The non-stationarity of the Hawkes process arise from the reproduction of an
event induced by the interaction terms described by the kernel function that
estimate the effect of the past events based on the elapsed time from them.
The intensity of the process is determined by summing the effects of all past
events through the kernel function. Additionally, by introducing a real number
”mark” that represents the strength of event effects, both the marks and the
kernel values can be considered in the intensity estimation. If the average
number of events triggered by one event, known as the branching ratio, exceeds
one, the system becomes non-stationary. If the branching ratio approaches one
from below, the system becomes critical, leading to power-law scaling in the
PDF of the intensities in the intermediate asymptotic region. In the sub-
critical case, the distribution decays exponentially beyond the characteristic
intensity, which diverges as the branching ratio approaches one (Kanazawa and
Sornette 2020b, a).

The discrete time self-exciting negative binomial distribution (DT-SE-NBD) was
proposed for the modeling of time-series data related to defaults (Hisakado
and Mori 2020). This model incorporates a parameter $\omega$ that captures the
correlation of events within the same term, and the process converges to a
discrete-time self-exciting Poisson process in the limit $\omega\to 0$. In the
context of time-series default data, defaults are typically recorded on a
yearly or quarterly basis, and defaults within the same period exhibit high
correlation. Consequently, the probability distribution of the number of
defaults displays overdispersion, where the variance is significantly larger
than the mean value. To accurately capture this overdispersion, the negative
binomial distribution (NBD), which has two parameters to control both the mean
and variance, is more suitable than the Poisson distribution (Hisakado et al.
2022a, b). In the continuous time limit, the DT-SE-NBD process transforms into
a point process with marks. As $\omega$ tends to zero, the marked point
process becomes a ”pure” Hawkes process as the value of the marks equals one.
The PDF of the intensities also exhibits power-law scaling and exponentially
decays beyond the characteristic intensity at the critical point and in the
subcritical case, respectively.

This study extends the results of the SE-NBD process with a single exponential
kernel to the case where the memory kernel is the sum of an arbitrary finite
number of exponentials. Empirically, power-law memory kernels have been
observed in studies on financial transaction data(Bacry et al. 2015). The
theoretical analysis of Hawkes processes with power-law memory kernels
requires decomposing the power-law memory kernel into a superposition of
exponential kernels with weights following the inverse gamma distribution. By
considering the case of multiple exponential kernels, we can estimate the
power-law exponent of the intensity distribution for the marked Hawkes process
with a power-law memory kernel. Additionally, we validate the theoretical
predictions for the power-law behavior of the intensities through numerical
simulations. The remaining sections of this paper are organized as follows: In
Section 2, we introduce an SE-NBD process with an arbitrary finite number of
exponentials and review related results on the process as well as on the
Hawkes process. Section 3 focuses on studying the multivariate stochastic
differential equation for the process, providing the solution to the two-
variate master equations with a double exponential memory kernel, and deriving
the PDF of the intensities. Furthermore, we discuss the PDF of the intensities
for the case where the memory kernel is a sum of finite $K$ exponential
functions. In Section 4, we verify the theoretical predictions for the power-
law exponents of the PDF of the intensities at the critical point through
numerical simulations. Finally, we present our conclusions in Section 5.

##  2 Model

We consider a DT-SE-NBD process $\\{X_{t}\\},t=1,\cdots$. The variable
$X_{t}\in\\{0,1,\cdots\\}$ represents the size of the event at time $t$. In
the context of modeling of time-series default data, $X_{t}$ specifically
represents the number of defaults that occur in the $t$-th period(Hisakado et
al. 2022a, b). $X_{t}$ obeys NBD for the condition
$\hat{\lambda}_{t}=\lambda_{t}$ as

$\displaystyle X_{t+1}$ $\displaystyle\sim$
$\displaystyle\mbox{NBD}\left(\alpha=\frac{\hat{\lambda}_{t}}{\omega},p=\frac{1}{\omega+1}\right),t\geq
0$ (1) $\displaystyle\hat{\lambda}_{t}$ $\displaystyle=$
$\displaystyle\nu_{0}+n\sum_{s=1}^{t}h_{t-s}X_{s},t\geq
1\,\,,\,\,\hat{\lambda}_{0}=\nu_{0}$ (2) $\displaystyle h_{t}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{k=1}^{K}n_{k}h^{k}_{t}=\frac{1}{n}\sum_{k=1}^{K}n_{k}(1-e^{-1/\tau_{k}})e^{-t/\tau_{k}},n=\sum_{k=1}^{K}n_{k}.$
(3)

Here, $\alpha,p$ are the parameters of the NBD. $\alpha>0$ represents the
number of successes until the experiment is terminated, and $p\in(0,1]$ is the
probability of success in each individual experiment. The sequence
$\\{h_{t}\\}_{t=0,1,\cdots}$ refers to the discount factor (memory kernel),
which is assumed to be normal; specifically, $\sum_{t=0}^{\infty}h_{t}=1$. The
prefactor $(1-e^{-/\tau_{k}})$ is introduced to ensure the normalization of
the $k$-th exponential term $h^{k}_{t}=(1-e^{-1/\tau_{k}})e^{-t/\tau_{k}}$,
where $\tau_{k}$ represents the memory length associated with the $k$-th term.
Consequently, $\sum_{t=0}^{\infty}(1-e^{-1/\tau_{k}})e^{-t/\tau_{k}}=1$. The
set of coefficients $\\{n_{k}\\}_{k=1,\cdots,K}$ quantifies the contribution
of the $k$-th exponential term to the memory kernel $h_{t},t=0,1\cdots$, with
$n$ representing their sum.

By considering a double-scaling limit, the DT-SE-NBD process can be
constructed as an extension of the multi-term Pólya urn process with a memory
kernel as follows:

$X_{t}\sim\lim_{\begin{subarray}{c}N,n_{0}\to\infty\\\
n_{0}/N=1\end{subarray}}\mbox{BBD}\left(N,\alpha=\frac{\hat{\lambda}_{t}}{\omega},\beta=\frac{n_{0}-\hat{\lambda}_{t}}{\omega}\right)=\mbox{NBD}\left(\alpha=\frac{\hat{\lambda}_{t}}{\omega},p=\frac{1}{\omega+1}\right).$

Here, BBD represents the beta-binomial distribution, and the variables
$N,\hat{\lambda}_{t},n_{0}-\hat{\lambda}_{t}$ represent the number of trials,
the number of red balls, and the number of blue balls in an urn in the $t$-th
term, respectively. In each term, we sequentially remove a ball and return it
along with $\omega$ additional balls of the same color. Consequently, the
total number of balls in the urn increases by $\omega$. The process is
repeated $N$ times within one term. The term $1/(1+n_{0}/\omega)$ represents
Pearson’s correlation coefficient, , which reflects the correlation between
the color choices of the balls in the process (Hisakado et al. 2006). The
feedback term $h_{t}$ induces an intertemporal correlation of the ball
choices. Furthermore, as $\omega$ approaches zero, $X_{t+1}$ under the
condition $\hat{\lambda}_{t}=\lambda_{t}$ follows a Poisson random variable
with a mean of $\lambda_{t}$, thereby transforming the DT-SE-NBD process into
a DT-SE Poisson process.

$X_{t+1}\sim\mbox{Po}(\hat{\lambda}_{t}),t\geq 0$

Fig.1 illustrates the transformation from the multi-term Pòlya urn process to
the DT-SE Negative binomial and DT-SE Poisson process.

Figure 1: Transformation of multi-term Pòlya urn process (beta-binomial) to
discrete-time Negative binomial and Poisson process.

We study the (unconditional) expected value of $X_{t+1}$. As
$X_{t+1}\sim\mbox{NBD}\left(\frac{\lambda_{t}}{\omega},\frac{1}{\omega+1}\right)$
for $\hat{\lambda}_{t}=\lambda_{t}$, we have

$E[X_{t+1}]=E_{\lambda_{t}}[E[X_{t+1}|\hat{\lambda}_{t}=\lambda_{t}]]=E_{\lambda_{t}}[\hat{\lambda}_{t}].$

where $E_{\lambda_{t}}[\,\,\,\,]$ is the ensemble average of the stochastic
process $\\{X_{s},s\leq t\\}$ [i.e., Filtration $F_{t}$]. We are now
interested in the steady state of the process and write the unconditional
expected value of $X_{t}$ as $\mu$. We have the following relationship:

$\mu\equiv\lim_{t\to\infty}E_{\lambda_{t}}[\hat{\lambda}_{t}].$

With $\lambda_{t}=\nu_{0}+n\sum_{s=1}^{t}X_{s}h_{t-s}$,$E[X_{s}]=\mu$ and the
normalization of the discount factor $h_{t}$, we have

$\mu=\nu_{0}+n\mu.$

We then obtain,

$\mu=\frac{\nu_{0}}{1-n}$

The phase transition between the steady and non-steady state occurs at
$n=n_{c}=1$, which is the critical point. This is the same as that in the
Hawkes process. In this sense, $n$ corresponds to the branching ratio of the
Hawkes process, and the steady-state condition is $n<1$.

We proceed to the continuous time SE-NBD process. We begin with the DT-SE-NBD
process in (1), (2), and (3). The stochastic process $\\{X_{t}\\},t=1,\cdots$
is non-Markovian, so we focus on the time evolution of
$\hat{z}_{t}\equiv\hat{\lambda}_{t}-\nu_{0}$. We decompose $\hat{z}_{t}$ as
follows:

$\displaystyle\hat{z}_{t}$ $\displaystyle\equiv$
$\displaystyle\hat{\lambda}_{t}-\nu_{0}=\sum_{k}\hat{z}^{k}_{t}$
$\displaystyle\hat{z}^{k}_{t}$ $\displaystyle\equiv$ $\displaystyle
n_{k}\sum_{s=1}^{t}h^{k}_{t-s}X_{s}\,\,,\,\,\
\hat{z}_{0}^{k}=0\color[rgb]{0,0,0}.$

$\hat{z}^{k}_{t}$ satisfies the following recursive equation:

$\hat{z}^{k}_{t+1}=e^{-1/\tau_{k}}\hat{z}^{k}_{t}+n_{k}h^{k}_{0}X_{t+1}.$

Here, we use the relation
$\sum_{s=1}^{t+1}X_{s}h^{k}_{t+1-s}=X_{t+1}h^{k}_{0}+e^{-1/\tau_{k}}\sum_{s=1}^{t}X_{s}h^{k}_{t-s}$.
The stochastic difference equation for $\hat{z}^{k}_{t}$ is

$\Delta\hat{z}^{k}_{t}\equiv\hat{z}^{k}_{t+1}-\hat{z}^{k}_{t}=(e^{-1/\tau_{k}}-1)\hat{z}^{k}_{t}+n_{k}h^{k}_{0}X_{t+1},k=1,\cdots,K.$
(4)

The DT-SE-NBD process is now represented by a multivariate stochastic
difference equation.

Now, we consider the continuous time limit. We divide the unit time interval
by infinitesimal time intervals with width $dt$. The decreasing factor
$e^{-1/\tau_{k}}$ for unit time period is replaced by $e^{-dt/\tau_{k}}\simeq
1-dt/\tau_{k}+o(dt/\tau_{k})$ for $dt$. In addition, we adopt the notation
$h_{k}(t)=\frac{1}{\tau_{k}}e^{-t/\tau_{k}}$. Here, we replace the
normalization factor $(1-e^{-1/\tau_{k}})$ of $h^{k}_{t}$ with $1/\tau_{k}$ in
$h_{k}(t)$ to ensure that
$\int_{0}^{\infty}h_{k}(t)dt=\frac{1}{\tau_{k}}\int_{0}^{\infty}e^{-t/\tau_{k}}=1$.
We introduce the continuous time memory kernel $h(t)$ as

$nh(t)=\sum_{k=1}^{K}n_{k}h_{k}(t)=\sum_{k=1}^{K}n_{k}\frac{1}{\tau_{k}}e^{-t/\tau_{k}}.$
(5)

We denote the continuous time limits of $\hat{z}^{k}_{t}$ and
$\hat{\lambda}_{t}$ as $\hat{z}_{k}(t)$ and $\hat{\lambda}(t)$, respectively.
The timing, size of events, and counting process are denoted by $\hat{t}_{i}$,
$m_{i}$ and $\hat{N}(t)$, respectively. $\hat{z}_{k}(t)$ and
$\hat{\lambda}(t)$ are defined as

$\displaystyle\hat{z}_{k}(t)$ $\displaystyle=$ $\displaystyle
n_{k}\sum_{i=1}^{\hat{N}(t)}h_{k}(t-\hat{t}_{i})\cdot m_{i}$
$\displaystyle\hat{\lambda}(t)$ $\displaystyle=$
$\displaystyle\nu_{0}+\sum_{k=1}^{K}\hat{z}_{k}(t)$

$X_{t}$ is the common noise for the time interval $[t,t+1)$ in the DT SE-NBD
process. To take the continuous time limit, it is necessary to define the
noise $d\xi(t)$ for the infinitesimal interval $[t,t+dt)$. There is freedom in
defining $d\xi(t)$, and we adopt the following, which is based on the
reproduction properties of NBD. When $X_{1}\sim
NBD(\alpha_{1},p),X_{2}\sim\mbox{NBD}(\alpha_{2},p)$,$X_{1}+X_{2}\sim\mbox{NBD}(\alpha_{1}+\alpha_{2},p)$.
We define $d\xi^{NBD}_{(\alpha,p)}(t)$ for the interval $[t,t+dt)$ as

$d\xi^{NBD}_{(\alpha,p)}(t)\sim\mbox{NBD}(\alpha dt,p).$

If we consider the above-mentioned noise, by the reproduction property, the
noise for the interval $[t,t+1)$ becomes

$\int_{t}^{t+1}d\xi^{NBD}_{(\alpha,p)}(s)\sim\mbox{NBD}(\alpha,p).$

In the self-exciting process, we replace $\alpha$ and $p$ by
$\frac{\hat{\lambda}(t)}{\omega}$ and $\frac{1}{\omega+1}$, respectively. The
conditional expected value and variance of
$d\xi^{NBD}_{(\alpha=\frac{\hat{\lambda}(t)}{\omega}dt,p=\frac{1}{\omega+1})}$
are

$\displaystyle
E\left[d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega}dt,\frac{1}{\omega+1})}\right|\left.\hat{\lambda}(t)=\lambda(t)\right]=\lambda(t)dt$
$\displaystyle
V\left[d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1})}\right|\left.\hat{\lambda}(t)=\lambda(t)\right]=\lambda(t)(\omega+1)dt.$

The conditional probability of occurrence of an event with size $m$ during the
time interval $[t,t+dt)$ is given as

$P\left(d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1})}=m\right|\left.\hat{\lambda}(t)=\lambda(t)\right)=\left\\{\begin{array}[]{cc}1-\frac{\lambda(t)}{\omega}\log(\omega+1)dt&m=0\\\
\frac{1}{m}\left(\frac{\lambda(t)}{\omega}\right)\left(\frac{\color[rgb]{1,0,0}\omega\color[rgb]{0,0,0}}{\omega+1}\right)^{m}dt&m\geq
1.\end{array}\right.$

In the limit $\omega\to 0$, the probabilities converge to

$\lim_{\omega\to
0}P\left(d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1})}=m\right|\left.\hat{\lambda}(t)=\lambda(t)\right)=\left\\{\begin{array}[]{cc}1-\lambda(t)dt&m=0\\\
\lambda(t)dt&m=1\\\ o(dt)&m\geq 2,\end{array}\right.$

and the event size $m$ is restricted to one.

The NBD noise $\xi^{NBD}(t)$ is then written as

$\xi^{NBD}_{\left(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1}\right)}(t)=\sum_{i=1}^{\hat{N}(t)}m_{i}\delta(t-\hat{t}_{i}).$

The probability mass function (PMF) for event size $m=1,\cdots$ is given by

$\rho(m)=\frac{P\left(d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1})}=m\right)}{P(d\xi^{NBD}>0)}=\frac{1}{m\log(\omega+1)}\left(\frac{\color[rgb]{1,0,0}\omega\color[rgb]{0,0,0}}{\omega+1}\right)^{m}.$
(6)

The intensity of the process, $P(d\xi^{NBD}>0)/dt$, is then given as

$\hat{\nu}(t)\equiv\frac{\log(\omega+1)}{\omega}\hat{\lambda}(t)\,\,,\,\,\hat{\lambda}(t)=\nu_{0}+\sum_{i=1}^{\hat{N}(t)}\sum_{k=1}^{K}n_{k}h_{k}(t-\hat{t}_{i})\cdot
m_{i}.$ (7)

This process is known as the marked Hawkes process. Each event has a distinct
influential power that is encoded in the marks
$\\{m_{i}\\}_{i=1,\cdots,\hat{N}(t)}$. $\\{m_{i}\\}_{i=1,\cdots,\hat{N}(t)}$
are IID random numbers and obey the PMF of (6). When the parameter $\omega$
tends to zero ($\omega\to 0$), the PMF $\rho(m)$ approaches the delta function
$\delta(m-1)$, and the marked Hawkes process reduces to the “pure” Hawkes
process. In the following analysis, we will focus on studying the PDF of the
intensities in the marked Hawkes process.

##  3 Solution

The multivariate difference equation (4) can be transformed into a
multivariate SDE as follows:

$\displaystyle\hat{z}_{k}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau_{k}}\hat{z}_{k}(t)dt+\frac{n_{k}}{\tau_{k}}d\xi^{NBD}_{(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1})},k=1,\cdots,K.$
(8) $\displaystyle\hat{\lambda}(t)$ $\displaystyle=$
$\displaystyle\nu_{0}+\sum_{k=1}^{K}\hat{z}_{k}(t)$ (9)
$\displaystyle\xi^{NBD}_{\left(\frac{\hat{\lambda}(t)}{\omega},\frac{1}{\omega+1}\right)}(t)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{\hat{N}(t)}m_{i}\delta(t-\hat{t}_{i}).$ (10)

Note that the same state-dependent NBD noise
$\xi^{\mathrm{NBD}}_{(\frac{\hat{\lambda}}{\omega},\frac{1}{\omega+1})}$
affects every component of the multivariate SDE
$\\{\hat{z}_{k}\\}_{k=1,...K}$. In other words, each shock event
simultaneously affects the trajectories for all excess intensities
$\\{\hat{z}_{k}\\}_{k=1,...K}$.

The formal solution of the SDE is

$\hat{z}_{k}(t)=\frac{n_{k}}{\tau_{k}}\int^{t}_{0}dt^{\prime}e^{-(t-t^{\prime})/\tau_{k}}\xi^{\mathrm{NBD}}_{(\frac{\hat{\lambda}}{\omega},\frac{1}{\omega+1})}(t^{\prime}).$

The SDE(8) can be interpreted as

$\hat{z}_{k}(t+dt)-\hat{z}_{k}(t)=\left\\{\begin{aligned}
-\frac{\hat{z}_{k}(t)}{\tau_{k}}dt\qquad&\mathrm{Prob.}=1-\frac{\hat{\lambda}(t)}{\omega}\log(\omega+1)dt\\\
\frac{n_{k}m}{\tau_{k}}\hskip
25.0pt&\mathrm{Prob.}=\frac{1}{m}\left(\frac{\lambda(t)}{\omega}\right)\left(\frac{\color[rgb]{1,0,0}\omega\color[rgb]{0,0,0}}{\omega+1}\right)^{m}dt,m\geq
1.\end{aligned}\right.$

We adopted the same procedure to derive the master equation for the PDF of
$\hat{z}_{k}$ in Ref.(Kanazawa and Sornette 2020b, a). As the SDEs for
$\boldsymbol{z}:=(z_{1},...,z_{K})$ are standard Markovian stochastic
processes, we obtain the corresponding master equation

$\displaystyle\frac{\partial}{\partial t}P_{t}$
$\displaystyle(\boldsymbol{z})=\sum^{K}_{k=1}\frac{\partial}{\partial
z_{k}}\frac{z_{k}}{\tau_{k}}P_{t}(\boldsymbol{z})+\sum_{m=1}^{\infty}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}$
$\displaystyle\times\left\\{\left[\nu_{0}+\sum^{K}_{k=1}\left(z_{k}-\frac{n_{k}m}{\tau_{k}}\right)\right]P_{t}(\boldsymbol{z-h})-\left[\nu_{0}+\sum^{K}_{k=1}z_{k}\right]P_{t}(\boldsymbol{z})\right\\}$
(11)

The jump-size vector is given by
$\boldsymbol{h}:=(n_{1}m/\tau_{1},...,n_{K}m/\tau_{K})$.

The master equation (3) takes a simplified form under the Laplace
representation:

$\tilde{P}_{t}(\boldsymbol{s}):=\mathcal{L}_{\color[rgb]{1,0,0}K\color[rgb]{0,0,0}}[P_{t}(\boldsymbol{z});\boldsymbol{s}],$
(12)

where the Laplace transformation in the $K$-dimensional space is defined as

$\mathcal{L}_{K}[P_{t}(\boldsymbol{z});\boldsymbol{s}]:=\int^{\infty}_{0}d\boldsymbol{z}e^{-\boldsymbol{s\cdot
z}}P_{t}(\boldsymbol{z})$ (13)

with the volume element $d\boldsymbol{z}:=\prod^{K}_{k=1}dz_{k}$. The wave
vector $\boldsymbol{s}:=(s_{1},...,s_{K})$ is the conjugate of the excess
intensity vector $\boldsymbol{z}:=(z_{1},...,z_{K})$.

The Laplace representation of the master equation (3) is given by

$\frac{\partial\tilde{P}_{t}(\boldsymbol{s})}{\partial
t}=-\sum^{K}_{k=1}\frac{s_{k}}{\tau_{k}}\frac{\partial\tilde{P}_{t}(\boldsymbol{s})}{\partial
s_{k}}+\sum^{\infty}_{m=1}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}\left(\nu_{0}-\sum^{K}_{k=1}\frac{\partial}{\partial
s_{k}}\right)(e^{-\boldsymbol{h\cdot s}}-1)\tilde{P}_{t}(\boldsymbol{s})$ (14)

Then, the Laplace representation (12) of $P_{t}(\boldsymbol{z})$, which is the
solution of (14), enables us to obtain the (one-dimensional) Laplace
representation $\tilde{Q}_{t}(s)$ of the intensity $\mathrm{PDF}\ P_{t}(\nu)$
according to

$\tilde{Q}_{t}(s):=\mathcal{L}[P_{t}(\nu);s]=e^{-\nu_{0}s}\tilde{P}_{t}(\boldsymbol{s}=(s,s,...,)).$
(15)

### A. Single exponential kernel

We now consider the case in which the memory function (5) consists of $K=1$
exponential functions(Hisakado et al. 2022a). The Laplace representation of
the master equation (3) is given by

$\frac{d\tilde{P}_{t}(s)}{dt}=-\frac{s}{\tau}\frac{d\tilde{P}_{t}(s)}{ds}+\sum^{\infty}_{m=1}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}\left(\lambda_{0}-\frac{d}{ds}\right)\left(e^{-\frac{nk}{\tau}s}-1\right)\tilde{P}_{t}(s).$

The steady-state PDF of the intensity $\hat{\lambda}(t)$ is

$P_{SS}(\lambda)\propto\lambda^{-1+2\frac{\nu_{0}\tau}{\omega+1}}e^{-\frac{2\tau\epsilon}{\omega+1}\lambda}\,\,,\,\,\epsilon=1-n.$
(16)

The power-law exponent of the PDF of the excess intensity is
$1-\frac{2\nu_{0}\tau}{\omega+1}$ and depends on $\omega$. In the limit
$\omega\rightarrow 0$, the result coincides with that in Ref.(Kanazawa and
Sornette 2020b, a). The power-law exponent increases with $\omega$ and
converges to $1$ in the limit $\omega\rightarrow\infty$. In addition, the
length scale beyond which the intensity shows exponential decay for the off-
critical case is $(\omega+1)/(2\tau\epsilon)$, and it diverges in the limit
$\epsilon=1-n\to 0$. This is also an increasing function of $\omega$.

### B. Double exponential kernel

We consider the case where the memory function (5) consists of $K=2$
exponential functions. To derive the solution for the Laplace representation
of the master equation (14), we apply the method of characteristics as in
Ref.(Kanazawa and Sornette 2020a). One can find a brief explanation in
Refs.(Gardiner 2009) and (Kanazawa and Sornette 2020a). In appendix A, we
provide a brief review.

We start from the Lagrange–Charpit equations, which are given by

$\displaystyle\frac{ds_{1}}{dl}$ $\displaystyle=$
$\displaystyle-\sum^{\infty}_{m=1}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}\left(e^{-\boldsymbol{h\cdot
s}}-1\right)-\frac{s_{1}}{\tau_{1}},$ (17) $\displaystyle\frac{ds_{2}}{dl}$
$\displaystyle=$
$\displaystyle-\sum^{\infty}_{m=1}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}\left(e^{-\boldsymbol{h\cdot
s}}-1\right)-\frac{s_{2}}{\tau_{2}},$ (18)
$\displaystyle\frac{d}{dl}\log\tilde{P}_{ss}$ $\displaystyle=$
$\displaystyle-\sum^{\infty}_{m=1}\frac{1}{m\omega}\left(\frac{\omega}{\omega+1}\right)^{m}\nu_{0}\left(e^{-\boldsymbol{h\cdot
s}}-1\right)$ (19)

and $l$ is the auxiliary “time” parameterizing the position on the
characteristic curve. Let us develop the stability analysis around $s=0$
(i.e., for large $\lambda^{\prime}s$) for this pseudo-dynamical system.

#### a. Sub-critical case $n<1$

Assuming $n:=n_{1}+n_{2}<1$, we first expand $e^{-\boldsymbol{h\cdot s}}$ to
compute the sum of $m$:

$e^{-\left(\frac{n_{1}s_{1}}{\tau_{1}}+\frac{n_{2}s_{2}}{\tau_{2}}\right)m}\simeq
1-\left(\frac{n_{1}s_{1}}{\tau_{1}}+\frac{n_{2}s_{2}}{\tau_{2}}\right)m+\frac{1}{2}\left(\frac{n_{1}s_{1}}{\tau_{1}}+\frac{n_{2}s_{2}}{\tau_{2}}\right)^{2}m^{2}+\cdots$
(21)

We obtain the linearized dynamics of the system (17),(18),(19) as follows:

$\frac{d\boldsymbol{s}}{dl}\simeq-\boldsymbol{H}\boldsymbol{s},\
\frac{d}{dl}\log\tilde{P}_{ss}\simeq\nu_{0}\boldsymbol{K}\boldsymbol{s}$ (22)

with

$\boldsymbol{H}:=\begin{pmatrix}\frac{1-n_{1}}{\tau_{1}}&\frac{-n_{2}}{\tau_{2}}\\\
\frac{-n_{1}}{\tau_{1}}&\frac{1-n_{2}}{\tau_{2}}\end{pmatrix},\
\boldsymbol{K}:=\left(\frac{n_{1}}{\tau_{1}},\frac{n_{2}}{\tau_{2}}\right).$

We introduce the eigenvalues $\beta_{1},\beta_{2}$ and eigenvectors
$\boldsymbol{e}_{1},\boldsymbol{e}_{2}$ of $\boldsymbol{H}$ such that

$\boldsymbol{P}:=(e_{1},e_{2}),\
\boldsymbol{P}^{-1}\boldsymbol{H}\boldsymbol{P}=\begin{pmatrix}\beta_{1}&0\\\
0&\beta_{2}\end{pmatrix}.$

The matrix $\boldsymbol{H}$ is the same with the results of Ref.(Kanazawa and
Sornette 2020b, a) and all eigenvalues are real. We denote them as
$\beta_{1}\geq\beta_{2}$. The determinant of $\boldsymbol{H}$ is given by

$\mathrm{det}\boldsymbol{H}=\frac{1-n}{\tau_{1}\tau_{2}}.$

This implies that the zero eigenvalue $\beta_{1}=0$ appears at the critical
point $n=1$. Below the critical point $n<1$, all eigenvalues are positive
($\beta_{1},\beta_{2}>0$). For $n<1$, the dynamics can be rewritten as

$\frac{d}{dl}\boldsymbol{P}^{-1}\boldsymbol{s}=-\begin{pmatrix}\beta_{1}&0\\\
0&\beta_{2}\end{pmatrix}\boldsymbol{P}^{-1}\boldsymbol{s}\Longrightarrow\boldsymbol{s}(l)=\boldsymbol{P}\begin{pmatrix}e^{-\beta_{1}(l-l_{0})}\\\
e^{-\beta_{2}(l-l_{0})}/C_{1}\end{pmatrix}.$

Here, $l_{0}$ and $C_{1}$ are integration constants. We can assume $l_{0}=0$
as the initial point of the characteristic curve without loss of generality.
Integrating the second equation in (22), we obtain

$\displaystyle\log\tilde{P}_{ss}$ $\displaystyle=$
$\displaystyle\nu_{0}\boldsymbol{K}\int\boldsymbol{s}(l)dl+C_{2}=-\nu_{0}\boldsymbol{KP}\begin{pmatrix}1/\beta_{1}&0\\\
0&1/\beta_{2}\end{pmatrix}\boldsymbol{P}^{-1}\boldsymbol{s}+C_{2}$
$\displaystyle=$ $\displaystyle-\nu\boldsymbol{KH}^{-1}\boldsymbol{s}+C_{2}.$

The general solution is given by

$\mathcal{H}(C_{1})=C_{2}$ (23)

with function $\mathcal{H}$ determined by the initial condition of the
characteristic curve. Let us introduce

$\bar{\boldsymbol{s}}:=\boldsymbol{P}^{-1}\boldsymbol{s}=\begin{pmatrix}\bar{s}_{1}\\\
\bar{s}_{2}\end{pmatrix}\Longrightarrow
C_{1}=(\bar{s}_{1})^{\beta_{2}/\beta_{1}}(\bar{s}_{2})^{-1}.$

This implies that the solution is given in the following form:

$\log\tilde{P}_{ss}(\boldsymbol{s})=-\nu\boldsymbol{KH}^{-1}\boldsymbol{s}+\mathcal{H}((\bar{s}_{1})^{\beta_{2}/\beta_{1}}(\bar{s}_{2})^{-1}).$

Owing to the renormalization of the PDF, the relation

$\lim_{\boldsymbol{s\rightarrow 0}}\log\tilde{P}_{ss}(\boldsymbol{s})=0$

must hold for any path in the ($s_{1},s_{2}$) space ending at the origin
($\mathrm{limit}\ \boldsymbol{s\rightarrow 0}$). Let us consider a specific
limit such that $\bar{s}_{2}=x^{-1}(\bar{s}_{1})^{\beta_{2}/\beta_{1}}$ and
$\bar{s}_{1}\rightarrow 0$ for an arbitrary positive $x$.

$\lim_{\bar{s}_{1}\rightarrow
0}\log\tilde{P}_{ss}(\boldsymbol{s})=\mathcal{H}(x)$

Because the left-hand side is zero for any $x$, the function
$\mathcal{H}(\cdot)$ must be exactly zero. Thus, this leads to

$\log\tilde{P}_{ss}(\boldsymbol{s})=-\nu\boldsymbol{KH}^{-1}\boldsymbol{s}.$

By substituting $\boldsymbol{s}=(s_{1}=s,s_{2}=s)$, we obtain

$\log\tilde{Q}_{ss}(\boldsymbol{s})=-\nu_{0}s+\log\tilde{P}_{t=\infty}(s(1,1))\simeq-\frac{\nu_{0}}{1-n}s.$

This is consistent with the asymptotic mean intensity in the steady state
(Kanazawa and Sornette 2020b, a).

#### b. Critical case $n=1$

In this case, the eigenvalues and eigenvectors of $\boldsymbol{H}$ are given
by

$\beta_{1}=0,\
\beta_{2}=\frac{n_{1}\tau_{1}+n_{2}\tau_{2}}{\tau_{1}\tau_{2}},\
\boldsymbol{e}_{1}=\begin{pmatrix}\tau_{1}\\\ \tau_{2}\end{pmatrix},\
\boldsymbol{e}_{2}=\begin{pmatrix}-n_{2}\\\ n_{1}\end{pmatrix}.$

This means that the eigenvalue matrix and its inverse matrix are given by

$\boldsymbol{P}=\begin{pmatrix}\tau_{1}&-n_{2}\\\
\tau_{2}&n_{1}\end{pmatrix},\
\boldsymbol{P}^{-1}=\frac{1}{\alpha}\begin{pmatrix}n_{1}&n_{2}\\\
-\tau_{2}&\tau_{1}\end{pmatrix},\
\alpha:=\det\boldsymbol{P}=\tau_{1}n_{1}+\tau_{2}n_{2}.$

Here, let us introduce

$\boldsymbol{X}=\begin{pmatrix}X\\\
Y\end{pmatrix}=\boldsymbol{P}^{-1}\boldsymbol{s},\Longleftrightarrow
X=\frac{n_{1}s_{1}+n_{2}s_{2}}{\alpha},\
Y=\frac{-\tau_{2}s_{1}+\tau_{1}s_{2}}{\alpha}.$ (24)

We then obtain

$\frac{dX}{dl}=0,\frac{dY}{dl}=-\beta_{2}Y$

at the leading linear order in the expansions of the powers of $X$ and $Y$.
Because the first linear term is zero in the dynamics of $X$, corresponding to
a transcritical bifurcation for the Lagrange–Charpit equations (16), we need
to consider the second-order term in $X$, namely,

$e^{-\left(\frac{n_{1}}{\tau_{1}}s_{1}+\frac{n_{2}}{\tau_{2}}s_{2}\right)m}\simeq
1-Xm+\frac{X^{2}m^{2}}{2}+n_{1}n_{2}\left(\frac{1}{\tau_{1}}-\frac{1}{\tau_{2}}\right)mY+\mathcal{O}(XY,X^{2}Y,Y^{2})$

where we dropped the terms of the order $Y^{2}$, $XY$, and $X^{2}Y$. We then
obtain the dynamic equations at the transcritical bifurcation to the leading
order:

$\frac{dY}{dl}\simeq-\beta_{2}Y,\
\frac{dX}{dl}\simeq-\frac{\omega+1}{2\alpha}X^{2}$

whose solutions are given by

$X(l)=\frac{2\alpha}{\omega+1}\frac{1}{l-l_{0}},\
Y(l)=C_{1}e^{-\beta_{2}(l-l_{0})}$

with constants of integration $l_{0}$ and $C_{1}$. We can assume $l_{0}=0$ is
the initial point on the characteristic curve. Remarkably, only the
contribution along the $X$ axis is dominant for the large $l$ limit (i.e.,
$|X|\gg|Y|$ for $l\rightarrow\infty$), which corresponds to the asymptotic
$\boldsymbol{s}\rightarrow 0$. We then obtain

$\displaystyle\log\tilde{P}_{ss}$ $\displaystyle\simeq$
$\displaystyle\nu_{0}\int
dl\left(\frac{n_{1}s_{1}(l)}{\tau_{1}}+\frac{n_{2}s_{2}(l)}{\tau_{2}}\right)$
$\displaystyle\simeq$ $\displaystyle-\frac{2\nu_{0}\alpha}{\omega+1}\log
X+\frac{\nu_{0}n_{1}n_{2}}{\beta_{2}}\left(\frac{1}{\tau_{1}}-\frac{1}{\tau_{2}}\right)Y+C_{2}$

with integration constant $C_{2}$. $C_{2}$ is a divergent constant because it
has to compensate the diverging logarithm $\log X$ to ensure that $\log
P_{ss}(s=0)=0$. This divergent constant appears as a result of neglecting the
ultraviolet (UV) cutoff for small $s$ (which corresponds to neglecting the
exponential tail of the PDF of intensities) (Kanazawa and Sornette 2020a).

Therefore, we obtain the steady solution

$\log\tilde{P}_{ss}(\boldsymbol{s})=-\frac{2\nu_{0}\alpha}{\omega+1}\log
X+\frac{\nu_{0}n_{1}n_{2}}{\beta_{2}}\left(\frac{1}{\tau_{1}}-\frac{1}{\tau_{2}}\right)Y$

for small $X$ and $Y$ by ignoring the UV cutoff and constant contribution.
This recovers the power-law formula of the intermediate asymptotic of the PDF
of the Hawkes intensities:

$\displaystyle\log\tilde{Q}_{ss}(s):$ $\displaystyle=$
$\displaystyle-\nu_{0}s+\log\tilde{P}_{ss}(s,s)\simeq-\frac{2\nu_{0}\alpha}{\omega+1}\log
s\quad(s\sim 0)$ (25) $\displaystyle\Longleftrightarrow$ $\displaystyle
P(\lambda)\sim\lambda^{-1+\frac{2\nu_{0}\alpha}{\omega+1}}\quad(\lambda\rightarrow\infty),$

where $\alpha=\tau_{1}n_{1}+\tau_{2}n_{2}$, as defined in (24). The power-law
exponent of the PDF of the excess intensity is
$1-\frac{2\nu_{0}\alpha}{\omega+1}$ and depends on $\omega$. In the limit
$\omega\rightarrow 0$, the result coincides with that in Ref.(Kanazawa and
Sornette 2020b, a).

### C. Discrete superposition of exponential kernels

We now study the case in which the memory kernel is the sum of $K$
exponentials for an arbitrary finite number $K$. We obtain the steady-state
PDF of the intensity $\hat{\lambda}(t)$ as

$\log\tilde{Q}_{ss}(s)\simeq-\frac{2\nu_{0}\alpha}{\omega+1}\log
s\Longleftrightarrow
P(\lambda)\sim\lambda^{-1+\frac{2\nu_{0}\alpha}{\omega+1}},\quad\mathrm{with}\
\alpha:=\sum^{K}_{k=1}n_{k}\tau_{k}.$ (26)

The power-law exponent of the PDF of the excess intensity is
$1-\frac{2\nu_{0}\alpha}{\omega+1}$ and depends on $\omega$. In the limit
$\omega\rightarrow 0$, the result coincides with that in Ref.(Kanazawa and
Sornette 2020b, a).

### D. General case

We study the case where the memory kernel is a continuous superposition of
exponential functions. We decompose the kernel as

$h(t)=\frac{1}{n}\int^{\infty}_{0}n(\tau)\frac{1}{\tau}e^{-t/\tau}d\tau\ ,\
n=\int^{\infty}_{0}n(\tau)d\tau.$

This decomposition satisfies the normalization condition
$\int^{\infty}_{0}h(t)=1$. The function $n(\tau)$ quantifies the contribution
of the exponential kernel $\frac{1}{\tau}e^{-t/\tau}$ to the branching ratio.
The steady-state PDF of the intensity $\hat{\lambda}_{t}$ is

$\log\tilde{Q}_{ss}(s)\simeq-\frac{2\nu_{0}\alpha}{\omega+1}\log
s\Longleftrightarrow
P(\lambda)\sim\lambda^{-1+\frac{2\nu_{0}\alpha}{\omega+1}},\quad\mathrm{with}\
\alpha:=\int^{\infty}_{0}d\tau n(\tau)\tau.$ (27)

The power-law exponent of the PDF of the excess intensity is
$1-\frac{2\nu_{0}\alpha}{\omega+1}$ and depends on $\omega$. In the limit
$\omega\rightarrow 0$, the result coincides with that in Ref.(Kanazawa and
Sornette 2020b, a).

We now consider the case in which the memory kernel exhibits power-law decay.
We express $h(t)$ as the superposition of $e^{-rt}$ with the weight of the
gamma distribution $f_{\gamma}(r)=r^{\gamma-1}e^{-r}/\Gamma(\gamma)$ as

$h(t)=\gamma(1+t)^{-(\gamma+1)}=\int^{\infty}_{0}re^{-rt}\cdot
f_{\gamma}(r)dr,\ \gamma>0.$

By the change of variable $r\rightarrow 1/\tau$, we rewrite $h(t)$ as the
superposition of $e^{-t/\tau}$ with the weight of the inverse gamma
distribution
$f^{\prime}_{\gamma}(\tau)=\tau^{-\gamma-1}e^{-1/\tau}/\Gamma(\gamma)$ as

$h(t)=\int^{\infty}_{0}\frac{1}{\tau}e^{-t/\tau}\cdot
f_{\gamma}(1/\tau)\cdot\frac{1}{\tau^{2}}d\tau=\int^{\infty}_{0}n(\tau)\frac{1}{\tau}e^{-t/\tau}d\tau\
,n(\tau)=f^{\prime}_{\gamma}(\tau).$ (28)

$\alpha$ of Eq.(27) is obtained as the expected value of the inverse gamma
distribution,

$\alpha=\int^{\infty}_{0}\tau
f^{\prime}_{\gamma}(\tau)d\tau=\frac{1}{\gamma-1},\ \gamma>1.$ (29)

The power-law exponent is $1-\frac{2\nu_{0}}{(\omega+1)(\gamma-1)}$.

##  4 Numerical verification

We conducted numerical simulations to study the steady-state PDF of the
intensity in the marked Hawkes process. In particular, our goal was to verify
the theoretical predictions of Eqs. (16), (25), (26), and (29). We employed
the time-rescaling theorem, which enables us to sample data more efficiently
than the rejection method. The time-rescaling theorem states that any point
process with an integrable conditional intensity function can be transformed
into a Poisson process with a unit rate (Brown et al. 2002). That is, it is
possible to convert a marked Hawkes process into a marked point process with a
constant intensity of 1 by performing time-dependent rescaling. A point
process with intensity 1 can easily generate the time of the next event
occurrence because the event interval follows an exponential distribution. The
time-rescaling theorem introduces the following time transformations:

$\Lambda(t)\equiv\int^{t}_{0}\hat{\nu}(s)ds.$ (30)

When event $\boldsymbol{t}_{n}=\\{t_{1},t_{2},...,t_{n}\\}$ follows a marked
Hawkes process with intensity function $\hat{\nu}(t)$ in the observation
interval $[0,T]$, event
$\boldsymbol{t}^{\prime}=\\{\Lambda(t_{1}),\Lambda(t_{2}),...,\Lambda(t_{n})\\}$
obtained by the time transformation (30) follows a marked point process with
intensity $1$ in the observation interval $[0,\Lambda(T)]$. The event interval
$\Lambda(t_{i})-\Lambda(t_{i-1})$ is given by

$\Lambda^{\prime}\equiv\Lambda(t_{i})-\Lambda(t_{i-1})=\int^{t_{i}}_{t_{i-1}}\hat{\nu}(s)ds\sim\mbox{Ex}(1).$

Hence, to numerically realize the marked Hawkes process, it is sufficient to
generate a random number $\Lambda^{\prime}$ following an exponential
distribution with mean 1 and successively find $t_{i}$ that satisfies the
above relation.

We rewrite the relation more concretely. For $t_{i-1}\leq s\leq t_{i}$,
$z_{k}(s)$ is written as

$z_{k}(s)=n_{k}\sum_{j=1}^{i-1}\frac{1}{\tau_{k}}e^{-(s-t_{j})/\tau_{k}}m_{j}=z_{k}(t_{i-1})e^{-(s-t_{i-1}/\tau_{k}}.$

The integral of $\nu(s)$ is

$\int_{t_{i-1}}^{t_{i}}\nu(s)ds=\frac{\ln(\omega+1)}{\omega}\left(\nu_{0}(t_{i}-t_{i-1})+\sum_{k=1}^{K}z_{k}(t_{i-1})\int_{t_{i-1}}^{t_{i}}e^{-(s-t_{i-1})/\tau_{k}}ds\right).$

We have to solve the next relation,

$\displaystyle\Lambda^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\ln(\omega+1)}{\omega}\left(\nu_{0}(t_{i}-t_{i-1})+\sum_{k=1}^{K}z_{k}(t_{i-1})\tau_{k}(1-e^{-(t_{i}-t_{i-1})/\tau_{k}}\right)$
(31) $\displaystyle\Longrightarrow$
$\displaystyle\nu_{0}(t_{i}-t_{i-1})+\sum_{k=1}^{K}z_{k}(t_{i-1})\tau_{k}\left(1-\mbox{exp}\left(-\frac{t_{i}-t_{i-1}}{\tau_{k}}\right)\right)-\frac{\omega\Lambda^{\prime}}{\ln(\omega+1)}=0$

Note that the jump size of this process is determined according to the
distribution of the mark $\rho(m)$ in Eq.(6) when the event occurrence time
$t_{i}$ is determined. The numerical procedures are given below.

Algorithm Sampling process of marked Hawkes process
We denote $z_{k,i}=z_{k}(t_{i})$.

Step 1 $i\leftarrow 1$.

Step 2 Get a random number $\Lambda^{\prime}\sim\mbox{Ex}(1)$.

Step 3 $t_{1}\leftarrow\frac{\Lambda^{\prime}}{\nu_{0}}$.

Step 4 Get a random number $m_{i}\sim\rho(m)$ in Eq.(6).

Step 5 $\hat{z}_{k,1}\leftarrow\frac{n_{k}m_{1}}{\tau_{k}},k=1,\cdots,K$.

Step 6 The following steps are repeated.

Step 6.1 $i\leftarrow i+1$.

Step 6.2 Get a random number $\Lambda^{\prime}\sim\mbox{Ex}(1)$.

Step 6.3 Solve Eq.(31) to get $t_{i}$.

Step 6.4 Get a random number $m_{i}\sim\rho(m)$ in Eq.(6).

Step 6.5
$\hat{z}_{k,i}\leftarrow\hat{z}_{k,i-1}e^{-\frac{t_{i}-t_{i-1}}{\tau_{k}}}+\frac{n_{k}m_{i}}{\tau_{k}},k=1,\cdots,K$.

We performed numerical simulations using Julia 1.7.3. The simulation time for
$t\in[0,T=10^{7}]$ is approximately 0.5 h for background intensity
$\nu_{0}=0.01$ and $K=1$. The execution environment is as follows: OS: Ubuntu
20.045 LTS, Memory: 64 GB, CPU: 4.7 GHz. The code is available on
github(Sakuraba 2023). We sampled the point process $\\{t_{i}\\},i=1,...,N(T)$
with $T=10^{7}$. The common settings of the sampling processes are
$\nu_{0}\in\\{0.01,0.2,1.0\\},\omega\in\\{0.01,1.0,10.0\\},\tau=1.0$ and
$n\in\\{0.999,0.99,0.9\\}$.

#### a. Single exponential kernel

Here, we test the theoretical prediction of Eq. (16),

$P_{SS}(\lambda)\propto\lambda^{-1+2\frac{\nu_{0}\tau}{\omega+1}}e^{-\frac{2\tau\epsilon}{\omega+1}\lambda}.$

Fig.2 presents the PDFs of $\hat{\lambda}(t)$ for the cases.

|  |
---|---|---
|  |
|  |
Figure 2: Plot of the steady state PDF $P_{ss}(\lambda)$ for the single
exponential kernel case $K=1$ near the critical point $n=n_{c}=1$. We set
$n=0.9(\bullet)$, $n=0.99)(\circ)$, and $n=0.999)(\bullet)$ for
$\tau=1.0,\omega\in\\{0.01,1.0,10.0\\}$ and $\nu_{0}\in\\{0.01,0.2,1.0\\}$.
(a-c) Background intensity $\nu_{0}=0.01$ and power-law exponents are
$0.9802$, $0.99$, and $0.9982$ for $\omega=0.01$, $1.0$, and $10$,
respectively. (d-f) Background intensity $\nu_{0}=0.2$ and power-law exponents
are $0.6039$, $0.8$, and $0.9636$ for $\omega=0.01$, $1.0$, and $10$,
respectively. (g-i) Background intensity $\nu_{0}=1.0$ and power-law exponents
are $-0.9802$, $0.0$, and $0.8182$ for $\omega=0.01$, $1.0$, and $10$,
respectively.

For a small background intensity $\nu_{0}=0.01\ll 1$, the power-law exponents
are $0.9802$ and $0.9982$ for $\omega=0.01$ and $\omega=10.0$, respectively.
These exponents are close to 1 and show little dependence on $\omega$. For
large $\nu_{0}$, the power-law exponent becomes more dependent on $\omega$.
When $\nu_{0}=1.0$, the power-law exponents are $-0.9802$ and $0.8182$ for
$\omega=0.01$ and $\omega=10.0$, respectively. When $\omega\gg 2\nu_{0}\tau$,
The power-law exponent becomes universally 1, as the expression
$1-2\nu_{0}\tau/(\omega+1)\simeq 1$ holds. Furthermore, the length scale of
the power-law region, denoted as $(\omega+1)/(2\tau\epsilon)$, increases as
$\omega$ increases. The numerical results confirm these findings by
demonstrating agreement between the slopes of the PDFs and the theoretical
values, and the exponents of the power-law were verified. In addition, by
observing the $x$-axis, it can be seen that the range of the straight power-
law region becomes wider as $\omega$ increases.

#### b. Double exponential kernel

We verified the theoretical predictions of Eq. (25).

$P_{SS}(\lambda)\propto\lambda^{-1+2\frac{\nu_{0}(\tau_{1}n_{1}+\tau_{2}n_{2})}{\omega+1}}.$

To realize $n=n_{1}+n_{2}\in\\{0.9,0.99,0.999\\}$, we set
$(n_{1},n_{2})=(0.5,0.4),(0.5,0.49)$ and $(0.5,0.499)$. Fig.3 presents the
double logarithmic plot of PDFs of $\hat{\lambda}_{t}$. One can observe the
agreement between the slopes of the PDFs and the theoretical values and the
theoretical prediction is verified.

|  |
---|---|---
|  |
|  |
Figure 3:  Plot of the steady-state PDF $P_{ss}(\lambda)$ for the double
exponential kernel case $K=2$ near the critical point $n=n_{c}=1$. We set
$n=0.0(\bullet)$, $n=0.99(\circ)$, and $n=0.999(\bullet)$ for
$(\tau_{1},\tau_{2})=(1.0,3.0),\omega\in\\{0.01,1.0,10.0\\}$, and
$\nu_{0}\in\\{0.01,0.2,1.0\\}$. (a-c) Background intensity $\nu_{0}=0.01$ and
power-law exponents are $0.9604$, $0.98$, and $0.9964$ for $\omega=0.01$,
$1.0$, and $10$, respectively. (d-f) Background intensity $\nu_{0}=0.2$ and
power-law exponents are $0.2079$, $0.6$, and $0.9273$ for $\omega=0.01$,
$1.0$, and $10$, respectively. (g-i) Background intensity $\nu_{0}=1.0$ and
power-law exponents are $-2.96$, $-1.0$, and $0.6363$ for $\omega=0.01$,
$1.0$, and $10$, respectively.

#### c. Triple exponential kernel

We verified the theoretical predictions of Eq. (26) for $K=3$.

$P_{SS}(\lambda)\propto\lambda^{-1+2\frac{\nu_{0}(\tau_{1}n_{1}+\tau_{2}n_{2}+\tau_{3}n_{3})}{\omega+1}}.$

To achieve $n=n_{1}+n_{2}+n_{3}\in\\{0.9,0.99,0.999\\}$, we set
$(n_{1},n_{2},n_{3})=(0.3,0.2,0.4),(0.3,0.2,0.49)$, and $(0.3,0.2,0.499)$.
Fig.4 shows the resulting PDFs of $\hat{\lambda}(t)$. One can observe the
agreement between the slopes of the PDFs and the theoretical values and the
theoretical prediction is verified.

|  |
---|---|---
|  |
|  |
Figure 4:  Plot of the numerical results for the steady-state PDF
$P_{ss}(\lambda)$ for the double exponential case $K=3$. $\lambda$ near the
critical point $n=n_{c}=1$. We set $n=0.9(\bullet)$, $n=0.99(\circ)$, and
$n=0.999(\bullet)$ for
$(\tau_{1},\tau_{2},\tau_{3})=(1.0,2.0,3.0),\omega\in\\{0.01,1.0,10.0\\}$, and
$\nu_{0}\in\\{0.01,0.2,1.0\\}$. (a-c) Background intensity $\nu_{0}=0.01$ and
power-law exponents are $0.956$, $0.978$, and $0.996$ for $\omega=0.01$,
$1.0$, and $10$, respectively. (d-f) Background intensity $\nu_{0}=0.2$ and
power-law exponents are $0.129$, $0.56$, and $0.92$ for $\omega=0.01$, $1.0$,
and $10$, respectively. (g-i) Background intensity $\nu_{0}=1.0$ and power-law
exponents are $-3.356$, $-1.2$, and $0.6$ for $\omega=0.01$, $1.0$, and $10$,
respectively.

#### d. Power-law memory kernel

We verified the theoretical predictions of Eq.(29) for $\gamma=11$.

$P(\lambda)\sim\lambda^{-1+\frac{2\nu_{0}}{10(\omega+1)}}.$

As $\tau\sim\mathrm{InvGamma}(\gamma,1)$ in (28), we set $K=100$ and
$\\{\tau_{i}\\}_{i=1,\cdots,K}$ as the $(i-1/2)$% point of the distribution.
To achieve $n=n_{1}+...+n_{K}\in\\{0.9,0.99,0.999\\}$, we set
$n_{i}=n/K,i=1,\cdots,K$. Fig.5 illustrates the resulting PDFs of
$\hat{\lambda}(t)$, confirming the validity of the theoretical prediction for
$\gamma=11$.

|  |
---|---|---
|  |
|  |
Figure 5:  Plot of the numerical results for the steady-state PDF
$P_{ss}(\lambda)$ for the general case $K=100$ and $\gamma=11$. We set
$n=0.9(\bullet)$, $n=0.99(\circ)$, and $n=0.999(\bullet)$ for
$(\tau_{1},...,\tau_{100})\sim\mathrm{InvGamma}(11,1),\omega\in\\{0.01,1.0,10.0\\}$,
and $\nu_{0}\in\\{0.01,0.2,1.0\\}$. (a-c) Background intensity $\nu_{0}=0.01$
and power-law exponents are $0.99802$, $0.999$, and $0.99982$ for
$\omega=0.01$, $1.0$, and $10$, respectively. (d-f) Background intensity
$\nu_{0}=0.2$ and power-law exponents are $0.9604$, $0.98$, and $0.9964$ for
$\omega=0.01$, $1.0$, and $10$, respectively. (g-i) Background intensity
$\nu_{0}=1.0$ and power-law exponents are $0.802$, $0.9$, and $0.982$ for
$\omega=0.01$, $1.0$, and $10$, respectively.

##  5 Conclusion

Here, we investigated the continuous time self-exciting negative binomial
process with a memory kernel consisting of a sum of finite $K$ exponential
functions. This process is equivalent to the marked Hawkes process. Our aim
was to extend the previous findings on the intermediate power-law behavior of
the intensities near the critical point for the steady and non-steady state
phase transition to the case of $K>1$. To achieve this, we developed an
efficient sampling method for the marked Hawkes process based on the time-
rescaling theorem. We conducted extensive simulations, focusing specifically
on the scenario where $K=100$ to capture the process with a power-law memory
kernel. Through these simulations, we were able to verify the theoretical
predictions and confirm their accuracy.

##  Appendix A Method of Characteristics

The method of characteristics is a standard method to solve first-order
partial differential equations (Gardiner 2009, Kanazawa and Sornette 2020b,
a). The equation for the characteristic function

$\phi(s)=\int^{\infty}_{-\infty}\mathrm{e^{isx}}p(x,t|x_{0},0)dx$

is

$\partial_{t}\phi+ks\partial_{s}\phi=-\frac{1}{2}Ds^{2}\phi.$ (A.1)

We consider the corresponding Lagrange–Charpit equations

$\frac{dt}{dl}=-1,\quad\frac{ds}{dl}=-ks,\quad\frac{d\phi}{dl}=\frac{1}{2}Ds^{2}\phi$

with the parameter $l$ encoding the position along the characteristic curves.
These equations are equivalent to an invariant form in terms of $l$

$\frac{dt}{1}=\frac{ds}{ks}=-\frac{d\phi}{\frac{1}{2}Ds^{2}\phi}.$

The method of characteristics can be used to solve this equation. Namely, if

$u(s,t,\phi)=a,\qquad\mathrm{and}\qquad v(s,t,\phi)=b$

are two integrals of the subsidiary equation (with $a$ and $b$ arbitrary
constants), then a general solution of (A.1) is given by

$f(u,v)=0.$

##  Appendix B Addendum

In a recent publication by K.Kanazawa and D.Sornette (Kanazawa and Sornette
2023), the power-law exponent of the intensity distribution of general marked
Hawkes process was given. We explain the correspondence for the reader’s
convenience.

In (Kanazawa and Sornette 2023), the power-law exponent is given using the
second moment of the mark’s PDF E$[m^{2}]$ as

$P_{SS}(\nu)\propto\nu^{-1-a},a=\frac{2\tau\nu_{0}}{\mbox{E}[m^{2}]}.$

The normalization of $\rho(m)$ with E$[m]=1$ is adopted. In our model, $a$ is
given as

$a=\frac{2\tau\nu_{0}}{\mbox{E}[m^{2}]/\mbox{E}[m]}.$

We use $\rho(m)$ in Eq.(6) to estimate the power-law exponent. As
E$[m]=\omega/\ln(\omega+1)$, E$[m^{2}]=\omega(\omega+1)/\ln(\omega+1)$, we
obtain

$P_{SS}(\nu)\propto\nu^{-1-a},a=\frac{2\tau\nu_{0}}{\mbox{E}[m^{2}]/\mbox{E}[m]}=2\tau\nu_{0}/(\omega+1).$

The result is consistent with ours in (Hisakado et al. 2022a).

###### Acknowledgements.

This work was supported by JPSJ KAKENHI [Grant No. 22K03445]. We would like to
thank Editage (www.editage.com) for English language editing.

## References

  * Bacry et al. (2015) Bacry E, Mastromatteo I, Muzy J (2015) Hawkes processes in finance. market microstructure and liquidity, 1, 1550005
  * Blanc et al. (2017) Blanc P, Donier J, Bouchaud JP (2017) Quadratic hawkes processes for financial prices. Quantitative Finance 17(2):171–188
  * Bowsher (2007) Bowsher CG (2007) Modelling security market events in continuous time: Intensity based, multivariate point process models. Journal of Econometrics 141(2):876–912
  * Brown et al. (2002) Brown EN, Barbieri R, Ventura V, Kass RE, Frank LM (2002) The time-rescaling theorem and its application to neural spike train data analysis. Neural computation 14(2):325–346
  * Engle and Russell (1998) Engle RF, Russell JR (1998) Autoregressive conditional duration: a new model for irregularly spaced transaction data. Econometrica pp 1127–1162
  * Errais et al. (2010) Errais E, Giesecke K, Goldberg LR (2010) Affine point processes and portfolio credit risk. SIAM Journal on Financial Mathematics 1(1):642–665
  * Filimonov and Sornette (2012) Filimonov V, Sornette D (2012) Quantifying reflexivity in financial markets: Toward a prediction of flash crashes. Physical Review E 85(5):056108
  * Filimonov and Sornette (2015) Filimonov V, Sornette D (2015) Apparent criticality and calibration issues in the hawkes self-excited point process model: application to high-frequency financial data. Quantitative Finance 15(8):1293–1314
  * Gardiner (2009) Gardiner C (2009) Stochastic Methods A Handbook for the Natural and Social Sciences. Springer, Berlin
  * Hasbrouck (1991) Hasbrouck J (1991) Measuring the information content of stock trades. The Journal of Finance 46(1):179–207
  * Hawkes (1971a) Hawkes A (1971a) Point spectra of some mutually exciting point processes. Journal of the Royal Statistical Society Series B 33, DOI 10.1111/j.2517-6161.1971.tb01530.x
  * Hawkes (1971b) Hawkes AG (1971b) Spectra of some self-exciting and mutually exciting point processes. Biometrika 58(1):83–90
  * Hawkes (2018) Hawkes AG (2018) Hawkes processes and their applications to finance: a review. Quantitative Finance 18(2):193–198
  * Hawkes and Oakes (1974) Hawkes AG, Oakes D (1974) A cluster process representation of a self-exciting process. Journal of Applied Probability 11(3):493–503
  * Hisakado and Mori (2020) Hisakado M, Mori S (2020) Phase transition in the bayesian estimation of the default portfolio. Physica A: Statistical Mechanics and its Applications 544:123480
  * Hisakado et al. (2006) Hisakado M, Kitsukawa K, Mori S (2006) Correlated binomial models and correlation structures. Journal of Physics A: Mathematical and General 39(50):15365
  * Hisakado et al. (2022a) Hisakado M, Hattori K, Mori S (2022a) From the multiterm urn model to the self-exciting negative binomial distribution and hawkes processes. Physical Review E 106(3):034106
  * Hisakado et al. (2022b) Hisakado M, Hattori K, Mori S (2022b) Multi-dimensional self-exciting nbd process and default portfolios. The Review of Socionetwork Strategies 16(2):493–512
  * Kanazawa and Sornette (2020a) Kanazawa K, Sornette D (2020a) Field master equation theory of the self-excited hawkes process. Physical Review Research 2(3):033442
  * Kanazawa and Sornette (2020b) Kanazawa K, Sornette D (2020b) Nonuniversal power law distribution of intensities of the self-excited hawkes process: A field-theoretical approach. Physical Review Letters 125(13):138301
  * Kanazawa and Sornette (2023) Kanazawa K, Sornette D (2023) Asymptotic solutions to nonlinear hawkes processes: A systematic classification of the steady-state solutions. Physical Review Research 5(1):013067
  * Kirchner (2017) Kirchner M (2017) An estimation procedure for the hawkes process. Quantitative Finance 17(4):571–595
  * Sakuraba (2023) Sakuraba K (2023) Self-exciting negative binomial distribution process and critical properties of intensity distribution. URL https://github.com/Kotaro-Sakuraba/SE-NBD_Process.git
  * Wheatley et al. (2019) Wheatley S, Wehrli A, Sornette D (2019) The endo–exo problem in high frequency financial price fluctuations and rejecting criticality. Quantitative Finance 19(7):1165–1178