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

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.

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

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:

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.

Fig.1 illustrates the transformation from the multi-term Pòlya urn process to the DT-SE Negative binomial and DT-SE 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

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:

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

We then obtain,

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:

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

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

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

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

$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

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

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

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

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

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

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

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

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

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.

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

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

The SDE(8) can be interpreted as

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

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:

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

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

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