2.1. Preliminary models and concepts

Introduced and pioneered through the series of papers [Reference Hawkes35–Reference Hawkes and Oakes37], the Hawkes process is a stochastic-intensity point process in which the current rate of arrivals is dependent on the history of the arrival process itself. Formally, this is defined as follows: let $(\lambda_t, N_{t,\lambda})$ be an intensity and counting process pair such that

\begin{align*}\mathbb{P}\left( N_{t+\Delta,\lambda} - N_{t,\lambda} = 1 \mid \mathcal{F}^N_t \right) &= \lambda_t \Delta + o(\Delta), \\\mathbb{P}\left( N_{t+\Delta,\lambda} - N_{t,\lambda} > 1 \mid \mathcal{F}^N_t \right) &= o(\Delta), \\\mathbb{P}\left( N_{t+\Delta,\lambda} - N_{t,\lambda} = 0 \mid \mathcal{F}^N_t \right) &= 1 - \lambda_t \Delta + o(\Delta),\end{align*}

where $\mathcal{F}^N_t$ is the filtration of $N_{t,\lambda}$ up to time t and $\lambda_t$ is given by

\begin{align*}\lambda_t = \lambda^* + \int_{-\infty}^t g(t-u) \textrm{d}N_{u,\lambda},\end{align*}

where $\lambda^* > 0$ and $g\, :\, \mathbb{R}^+ \to \mathbb{R}^+$ is such that $\int_0^\infty g(x) \textrm{d}x < 1$ . Through this definition, the intensity $\lambda_t$ captures the history of the arrival process up to time t. Thus, $\lambda_t$ encapsulates the sequence of past events and uses it to determine the rate of future occurrences. We refer to $\lambda^*$ as the baseline intensity and $g({\cdot})$ as the excitation kernel. The baseline intensity represents an underlying stationary arrival rate, and the excitation kernel governs the effect that the history of the process has on the current intensity. A common modeling choice is to set $g(x) = \alpha e^{-\beta x}$ , where $\beta > \alpha > 0$ . This is often referred to as the ‘exponential’ kernel, and it is perhaps the most widely used form of the Hawkes process. In this case, $(\lambda_t, N_{t,\lambda})$ is a Markov process obeying the stochastic differential equation

\begin{align*}\textrm{d}\lambda_t = \beta(\lambda^* - \lambda_t)\textrm{d}t + \alpha \textrm{d}N_{t,\lambda}.\end{align*}

That is, at arrival epochs $\lambda_t$ jumps upward by the amount $\alpha$ and the $N_{t,\lambda}$ increases by 1; between arrivals $\lambda_t$ decays exponentially at rate $\beta$ towards the baseline intensity $\lambda^*$ . Thus, each arrival increases the likelihood of additional arrivals occurring soon afterwards—hence, it self-excites. This form of the Hawkes process is also often alternatively stated with an initial value for $\lambda_t$ , say $\lambda_0 \geq \lambda^*$ . In this case, the intensity can be expressed as

\begin{align*}\lambda_t=\lambda^*+(\lambda_0 - \lambda^*)e^{-\beta t}+\alpha\int_{0}^t e^{-\beta (t-u)} \textrm{d}N_{u,\lambda}.\end{align*}

Additional overview of the Hawkes process with the exponential kernel can be found in Section 2 of [Reference Daw and Pender20]. Another common choice for excitation kernel is the ‘power-law’ kernel $g(x) = \frac{k}{(c+x)^p}$ , where $k > 0$ , $c > 0$ , and $p > 0$ . This kernel was originally popularized in seismology [Reference Ogata51].

2.2. Defining the ephemerally self-exciting process

As we have discussed in the introduction, a plethora of natural phenomena exhibit self-exciting features but only for a finite amount of time. This prompts the notion of ephemeral self-excitement. By comparison to the traditional Hawkes process we have reviewed in Subsection 2.1, we seek a model in which a new occurrence increases the arrival rate only so long as the newly entered entity remains active in the system. Thus, we now define the ephemerally self-exciting process (ESEP), which trades the Hawkes process’s eternal decay for randomly drawn expiration times. Moreover, in the following Markovian model, exponential decay is replaced with exponentially distributed durations. In Section 4, we extend these concepts to generally distributed service. As another generalization, in Appendix B we consider a Markovian model with both decay and down-jumps. For now, we explore the effects of ephemerality through the ESEP model in Definition 1.

Definition 1. (Ephemerally self-exciting process.) For times $t \geq 0$ , a baseline intensity $\eta^* > 0$ , intensity jump size $\alpha > 0$ , and expiration rate $\beta > 0$ , let $N_t$ be a counting process with stochastic intensity $\eta_t$ such that

(1) \begin{align}\eta_t = \eta^* + \alpha Q_t,\end{align}

where $Q_t$ is incremented with $N_t$ and then is depleted at unit down-jumps according to the rate $\beta Q_t$ . We say that $(\eta_t, N_t)$ is an ephemerally self-exciting process (ESEP).

We will assume that $\eta_0$ and $Q_0$ are known initial values such that $\eta_0 = \eta^* + \alpha Q_0$ . In addition to this definition, one could also describe the ESEP through its dynamics. In particular, the behavior of this process can be summarily cast through the life cycle of its arrivals:

1. i. At each arrival, the arrival rate $\eta_t$ increases by $\alpha$ .

2. ii. Each arrival remains active for an activity duration drawn from an independent and identically distributed (i.i.d.) sequence of exponential random variables with rate $\beta$ .

3. iii. At the expiration of a given activity duration, $\eta_t$ decreases by $\alpha$ .

The ephemerality of the ESEP is embodied by this cycle. Because arrivals only contribute to the intensity for the length of their activity duration, their effect on the process’s excitation vanishes when this clock expires. Furthermore, there is an affine relationship between the number of active ‘exciters’—meaning unexpired arrivals still causing excitation—and the intensity, i.e. $\eta_t = \eta^* + \alpha Q_t$ . Thus, we could also track the arrival rate through $Q_t$ in place of $\eta_t$ and still have full understanding of this process. This also means that results are readily transferrable between these two processes; we will often make use of this fact.

Because the ESEP is quite parsimonious, there are many alternative perspectives we could take to gain additional understanding of it. For example, one could consider $Q_t$ a Markovian queueing system with infinitely many servers and a state-dependent arrival rate. Equivalently, one could also describe the ESEP as a Markov chain on the nonnegative integers where transitions at state i are to $i+1$ at rate $\eta^* + \alpha i$ and to $i-1$ at rate $\mu i$ , with the counting process then defined as the epochs of the upward jumps in this chain. This Markov chain perspective certifies the existence and uniqueness of Definition 1. A visualization of this linear birth–death–immigration process is given in Figure 1. Stability for this chain occurs when $\beta > \alpha$ ; we will assume this henceforth, although of course it is not necessary for transient results. One could also view the ESEP as a generalization of Hawkes’s original definition where the excitation kernel function $g({\cdot})$ is replaced with a randomly drawn indicator function that is different for each arrival. Each indicator function compares time to an independently drawn exponential random variable, and this perspective is closely aligned with our analysis in Section 3.

Figure 1. The transition diagram of the Markov chain for $Q_t$ .

In the remainder of this subsection, let us now develop a few fundamental quantities for this stochastic process, particularly its intensity and active number in system, as these capture the self-exciting behavior of the process. First, in Proposition 1 we compute the transient moment generating function for the intensity $\eta_t$ . As we have noted, this can also be used to immediately derive the same transform for $Q_t$ , and the proof makes use of this fact.

Proposition 1. Let $\eta_t = \eta^* + \alpha Q_t$ be the intensity of an ESEP with baseline intensity $\eta^* > 0$ , intensity jump $\alpha > 0$ , and expiration rate $\beta > \alpha$ . Then the moment generating function for $\eta_t$ is given by

\begin{align*}{\mathbb{E}\left[e^{\theta \eta_t}\right]}&=\left(\frac{\beta - \alpha e^{\alpha\theta} - \beta \big(1-e^{\alpha\theta}\big)e^{-(\beta - \alpha)t}}{\beta - \alpha e^{\alpha\theta} - \alpha \big(1-e^{\alpha\theta}\big)e^{-(\beta - \alpha)t}}\right)^{\!\frac{\eta_0 - \eta^*}{\alpha}}\\&\qquad \cdot \left(\frac{\beta e^{\alpha \theta}}{\beta - \alpha e^{\alpha\theta}}-\frac{\alpha e^{\alpha \theta}}{\beta - \alpha e^{\alpha\theta}}\left(\frac{\beta - \alpha e^{\alpha\theta} - \beta \big(1-e^{\alpha\theta}\big)e^{-(\beta - \alpha)t}}{\beta - \alpha e^{\alpha\theta} - \alpha \big(1-e^{\alpha\theta}\big)e^{-(\beta - \alpha)t}}\right)\right)^{\!\frac{\eta^*}{\alpha}},\end{align*}

for all $t \geq 0$ and $\theta < \frac{1}{\alpha}\log\left(\frac{\beta}{\alpha}\right)$ .

Proof. We will approach this through the perspective of the active number in system, $Q_t$ . Using Lemma 2, we have that the probability generating function for $Q_t$ , say $\mathcal{P}(z,t) = {\mathbb{E}\!\left[z^{Q_t}\right]}$ for $z \in [0,1]$ , is given by the solution to the following partial differential equation (PDE):

\begin{align*}\frac{\partial}{\partial t}{\mathbb{E}\big[z^{Q_t}\big]}&={\mathbb{E}\!\left[\left(\eta^* + \alpha Q_t\right)\left(z^2 - z\right)z^{Q_t - 1} + \beta Q_t \left(1 - z\right)z^{Q_t - 1}\right]},\end{align*}

which is equivalently expressed by

\begin{align*}\frac{\partial}{\partial t}\mathcal{P}(z,t)&=\eta^* \left(z - 1\right)\mathcal{P}(z,t)+\left(\alpha \left(z^2 - z\right)+\beta\! \left(1 - z\right)\right)\frac{\partial}{\partial z}\mathcal{P}(z,t),\end{align*}

with initial condition $\mathcal{P}(z,0) = z^{Q_0}$ . The solution to this initial value problem is given by

\begin{align*}\mathcal{P}(z,t)&=\left(\frac{\beta - \alpha z - \beta (1-z)e^{-(\beta - \alpha)t}}{\beta - \alpha z - \alpha (1-z)e^{-(\beta - \alpha)t}}\right)^{Q_0}\\&\qquad \cdot \left(\frac{\beta}{\beta - \alpha z}-\frac{\alpha}{\beta - \alpha z}\left(\frac{\beta - \alpha z - \beta (1-z)e^{-(\beta - \alpha)t}}{\beta - \alpha z - \alpha (1-z)e^{-(\beta - \alpha)t}}\right)\right)^{\!\frac{\eta^*}{\alpha}},\end{align*}

yielding the probability generating function for $Q_t$ . By setting $z = e^{\theta}$ we obtain the moment generating function. Finally, using the affine relationship $\eta_t = \eta^* + \alpha Q_t$ , we have that

\begin{align*}{\mathbb{E}\!\left[e^{\theta \eta_t}\right]}={\mathbb{E}\!\left[e^{\theta \left(\eta^* + \alpha Q_t\right)}\right]}=e^{\theta \eta^*}{\mathbb{E}\!\left[e^{\alpha \theta Q_t}\right]},\end{align*}

with $\eta_0 = \eta^* + \alpha Q_0$ .

As we have mentioned, this Markov chain can be shown to be stable for $\beta > \alpha$ through standard techniques. Thus, using the moment generating function from Proposition 1, we can find the steady-state distributions of the intensity and the active number in system by taking the limit of t. We can quickly observe that this leads to a negative binomial distribution, as we state in Theorem 1. Because of the varying definitions of the negative binomial distribution, we state the probability mass function explicitly.

Theorem 1. Let $\eta_t = \eta^* + \alpha Q_t$ be an ESEP with baseline intensity $\eta^* > 0$ , intensity jump $\alpha > 0$ , and expiration rate $\beta > \alpha$ . Then the active number in system in steady state follows a negative binomial distribution with parameters $\frac{\alpha}{\beta}$ and $\frac{\eta^*}{\alpha}$ , which is to say that the steady-state probability mass function is

(2) \begin{align}\mathbb{P}\left( Q_\infty = k \right)&=\frac{\Gamma\left(k + \frac{\eta^*}{\alpha}\right)}{\Gamma\left(\frac{\eta^*}{\alpha}\right) k! }\left(\frac{\beta - \alpha}{\beta}\right)^{\!\frac{\eta^*}{\alpha}}\left(\frac{\alpha}{\beta}\right)^k.\end{align}

Consequently, the steady-state distribution of the intensity is given by a shifted and scaled negative binomial with parameters $\frac{\alpha}{\beta}$ and $\frac{\eta^*}{\alpha}$ , shifted by $\eta^*$ and scaled by $\alpha$ .

Proof. Using Proposition 1, we can see that the steady-state moment generating function of $Q_t$ is given by

\begin{align*}\lim_{t\to\infty} {\mathbb{E}\!\left[e^{\theta Q_t}\right]}=\left(\frac{\beta-\alpha}{\beta - \alpha e^{\theta}}\right)^{\!\frac{\eta^*}{\alpha}}.\end{align*}

We can observe that this steady-state moment generating function is equivalent to that of a negative binomial. By the affine transformation $\eta_t = \eta^* + \alpha Q_t$ , we find the steady-state distribution for the intensity.

The negative binomial distribution presented in Theorem 1 allows for $\frac{\eta^*}{\alpha}$ to take on any positive real value; integrality is not required. If $\frac{\eta^*}{\alpha}$ is in fact an integer, then the gamma functions will match the corresponding factorial functions and this ratio of factorials will simplify to a binomial coefficient, reproducing the most familiar form of the negative binomial distribution. In this case, $\frac{\eta^*}{\alpha}$ can be interpreted as the number of failures upon which an experiment is stopped, where the success probability is $\frac{\alpha}{\beta}$ and $k + \frac{\eta^*}{\alpha}$ is the total number of trials.

Let us pause to note that this explicit characterization of the steady-state intensity is already an advantage of the ESEP over the traditional Markovian Hawkes process, for which there is not a closed-form intensity stationary distribution available. As a consequence of Theorem 1, we can observe that the steady-state mean of the intensity is $\eta_\infty \,:\!=\, \frac{\beta \eta^*}{\beta - \alpha}$ . Interestingly, this would also be the steady-state mean of the Hawkes process when given the same baseline intensity, the same intensity jump size, and an exponential decay rate equal to the rate of expiration. This leads us to ponder how the processes would otherwise compare when given equivalent parameters. In Proposition 2 we find that although this equivalence of means extends to transient settings, for all higher moments the ESEP dominates the Hawkes process in terms of both the intensity and the counting process.

Proposition 2. Let $(\eta_t, N_{t,\eta})$ be an ESEP intensity and counting process pair with jump size $\alpha > 0$ , expiration rate $\beta > \alpha$ , and baseline intensity $\eta^* > 0$ . Similarly, let $(\lambda_t, N_{t,\lambda})$ be a Hawkes process intensity and counting process pair with jump size $\alpha > 0$ , decay rate $\beta > 0$ , and baseline intensity $\eta^* > 0$ . Then, if the two processes have equal initial values, their means will satisfy

(3) \begin{align}&{\mathbb{E}\!\left[\lambda_t\right]} = {\mathbb{E}\!\left[\eta_t\right]},&&{\mathbb{E}\!\left[N_{t,\lambda}\right]} = {\mathbb{E}\!\left[N_{t,\eta}\right]},\end{align}

and for $m \geq 2$ their mth moments are ordered so that

(4) \begin{align}&{\mathbb{E}\!\left[\lambda_t^m\right]} \leq {\mathbb{E}\!\left[\eta_t^m\right]},&&{\mathbb{E}\!\left[N_{t,\lambda}^m\right]} \leq {\mathbb{E}\!\left[N_{t,\eta}^m\right]},\end{align}

for all time $t \geq 0$ .

Proof. Let us start with the means. For the intensities, we can note that these are given by the solutions to

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\lambda_t\right]}=\alpha {\mathbb{E}\!\left[\lambda_t\right]} - \beta({\mathbb{E}\!\left[\lambda_t\right]} - \eta^*)\qquad \text{and}\qquad \frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\eta_t\right]}=\alpha{\mathbb{E}\!\left[\eta_t\right]} - \beta \alpha\left( \frac{{\mathbb{E}\!\left[\eta_t\right]} - \eta^*}{\alpha} \right),\end{align*}

and through simplification one can quickly observe that these two ordinary differential equations (ODEs) are equivalent. Thus, because we have assumed that the processes have the same initial values, we find that ${\mathbb{E}\!\left[\lambda_t\right]} = {\mathbb{E}\!\left[\eta_t\right]}$ . Since

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[N_{t,\lambda}\right]} = {\mathbb{E}\!\left[\lambda_t\right]}\qquad \text{and}\qquad \frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[N_{t,\eta}\right]} = {\mathbb{E}\!\left[\eta_t\right]},\end{align*}

this equality immediately extends to the means of the counting processes as well. We will now use these equations for the means as the base cases for inductive arguments, beginning again with the intensities. For the inductive step, we will assume that the intensity moment ordering holds for moments 1 to $m-1$ . The mth moment of the Hawkes process intensity is thus given by the solution to

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[{\lambda_t}^m\right]}=\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[{\lambda_t}^{k+1}\right]} \alpha^{m-k}-m \beta {\mathbb{E}\!\left[{\lambda_t}^m\right]}+m \beta \eta^{*} {\mathbb{E}\!\left[{\lambda_{t}}^{m-1}\right]}\,:\!=\,f_{\lambda}\left(t, {\mathbb{E}\!\left[{\lambda_t}^{m}\right]}\right),\end{align*}

where $f_\lambda\!\left(t, {\mathbb{E}\!\left[\lambda_t^{m}\right]}\right)$ is meant to capture that this ODE depends on the value of the mth moment and of the lower moments, which by the inductive hypothesis we take as known functions of the time t. Then, the mth moment of the ESEP intensity will satisfy

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[{\eta_t}^m\right]}=&\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[{\eta_t}^{k+1}\right]} \alpha^{m-k}+\frac{\beta}{\alpha}\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[{\eta_t}^{k+1}\right]} ({-}\alpha)^{m-k}\\&\quad -\frac{\beta \eta^*}{\alpha}\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[{\eta_t}^{k}\right]} ({-}\alpha)^{m-k}.\end{align*}

By pulling the $k = m-1$ terms off the top of each summation, we can re-express this ODE as

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\eta_t^m\right]}&=\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[\eta_t^{k+1}\right]} \alpha^{m-k}-m\beta {\mathbb{E}\!\left[\eta_t^{m}\right]}+m\beta \eta^* {\mathbb{E}\!\left[\eta_t^{m-1}\right]}\\&\quad+\frac{\beta}{\alpha}\sum_{k=0}^{m-2} {\left(\begin{array}{c}m\\ k\end{array}\right)} ({-}\alpha)^{m-k}\left({\mathbb{E}\!\left[\eta_t^{k+1}\right]}-\eta^* {\mathbb{E}\!\left[\eta_t^{k}\right]}\right),\end{align*}

and through the definition of $f_\lambda({\cdot})$ and the inductive hypothesis, we can find the following lower bound:

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\eta_t^m\right]}&\geq f_\lambda\!\left(t, {\mathbb{E}\!\left[\eta_t^m\right]}\right)+\frac{\beta}{\alpha}\sum_{k=0}^{m-2} {\left(\begin{array}{c}m\\ k\end{array}\right)} ({-}\alpha)^{m-k}\left({\mathbb{E}\!\left[\eta_t^{k+1}\right]}-\eta^* {\mathbb{E}\!\left[\eta_t^{k}\right]}\right)\!.\end{align*}

This rightmost term can then be expressed as

\begin{align*}&\frac{\beta}{\alpha}\sum_{k=0}^{m-2} {\left(\begin{array}{c}m\\ k\end{array}\right)} ({-}\alpha)^{m-k}\left({\mathbb{E}\!\left[\eta_t^{k+1}\right]}-\eta^* {\mathbb{E}\!\left[\eta_t^{k}\right]}\right)\\&\qquad =\frac{\beta}{\alpha}{\mathbb{E}\!\left[(\eta_t - \eta^*)\left((\eta_t - \alpha)^m-\eta_t^m+m \alpha \eta_t^{m-1}\right)\right]},\end{align*}

and we can now reason about the quantity inside the expectation. By definition, we have that $\eta_t \geq \eta^*$ with probability one, and furthermore we can observe that if $\eta_t - \eta^* > 0$ , then $\eta_t \geq \eta^* + \alpha > \alpha$ . Thus, let us consider $(\eta_t - \alpha)^m-\eta_t^m+m \alpha \eta_t^{m-1}$ assuming $\eta_t > \alpha$ . Dividing through by $\eta_t^m$ , we have the expression

(5) \begin{align}\left(1 - \frac{\alpha}{\eta_t}\right)^m-1+\frac{m \alpha}{ \eta_t}.\end{align}

Since $(1-x)^m - 1 + mx$ is equal to 0 at $x = 0$ and is non-decreasing on $x \in [0,1)$ via a first-derivative check, we can note that (5) is nonnegative for all $\eta_t > \eta^*$ . Thus, we have that

\begin{align*}{\mathbb{E}\!\left[(\eta_t - \eta^*)\left((\eta_t - \alpha)^m-\eta_t^m+m \alpha \eta_t^{m-1}\right)\right]}\geq 0,\end{align*}

and by consequence,

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[\eta_t^m\right]}\geq f_\lambda\!\left(t, {\mathbb{E}\!\left[\eta_t^m\right]}\right),\end{align*}

completing the proof of the intensity moment ordering via Lemma 3. For the counting processes, let us again assume as an inductive hypothesis that the moment ordering holds for moments 1 through $m-1$ , with the mean equality serving as the base case. Then, the mth moment of the ESEP counting process will satisfy

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[N_{t,\eta}^m\right]}&=\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[\eta_t N_{t,\eta}^k\right]},\end{align*}

and the ODE for the mth Hawkes counting process moment is analogous. By the Fortuin–Kasteleyn–Ginibre inequality, we can observe that

\begin{align*}\sum_{k=0}^{m-1} {\left(\begin{array}{c}m \\ k\end{array}\right)} {\mathbb{E}\!\left[\eta_t N_{t,\eta}^k\right]}\geq \sum_{k=0}^{m-1} {\left(\begin{array}{c}m \\ k\end{array}\right)} {\mathbb{E}\!\left[\eta_t\right]}\, {\mathbb{E}\!\left[N_{t,\eta}^k\right]}.\end{align*}

By the inductive hypothesis, we have that ${\mathbb{E}\!\left[N_{t,\eta}^k\right]} \geq {\mathbb{E}\!\left[N_{t,\lambda}^k\right]}$ for each $k \leq m - 1$ , and thus we can observe that

\begin{align*}\sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[\eta_t\right]}\, {\mathbb{E}\!\left[N_{t,\eta}^k\right]}\geq \sum_{k=0}^{m-1} {\left(\begin{array}{c}m\\ k\end{array}\right)} {\mathbb{E}\!\left[\lambda_t\right]}\, {\mathbb{E}\!\left[N_{t,\lambda}^k\right]}.\end{align*}

Finally, by another application of Lemma 3, we have ${\mathbb{E}\!\left[N_{t,\lambda}^m\right]} \leq {\mathbb{E}\!\left[N_{t,\eta}^m\right]}$ .

The fact that the ESEP variance dominates the Hawkes variance should not be surprising, since the presence of both up- and down-jumps means that the ESEP sample paths should be subject to more abrupt changes. Nevertheless, this also shows that the ESEP is more over-dispersed than the Hawkes process is. This may be an attractive feature for data modeling. It is worth noting that matrix computations are available for all moments of these intensities via [Reference Daw and Pender21], through which one could use the method of moments to fit the processes to data.

2.3. The ephemerally self-exciting counting process

Thus far we have studied the intensity of the ESEP, as this process is by definition tracking the self-excitement. However, this excitation is manifested in the actual arrivals from the process, which are counted in $N_t$ . We now turn our attention to developing fundamental quantities for this counting process. To begin, we give the probability generating function of the counting process in closed form below in Proposition 3. One can note that by comparison, the generating functions of the Hawkes process are instead only expressible as functions of ODEs with no known closed-form solutions; see for example Subsection 3.5 of [Reference Daw and Pender20].

Proposition 3. Let $N_t$ be the number of arrivals by time $t \geq 0$ in an ESEP with baseline intensity $\eta^* > 0$ , intensity jump $\alpha > 0$ , and expiration rate $\beta > \alpha$ . Then, for $z \in [0,1]$ , the probability generating function of $N_t$ is given by

(6) \begin{align}&{\mathbb{E}\!\left[z^{N_t}\right]}=e^{\frac{\eta^*(\beta-\alpha)}{2\alpha}t}\left(\frac{2e^{\frac{t}{2} \sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}}{1-\frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}+\left(1 + \frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}\right)e^{t \sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}}\right)^{\!\frac{\eta^*}{\alpha}}\nonumber \\&\quad\cdot \left(\!\frac{\beta + \alpha}{2\alpha}+\frac{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}{2\alpha}\!\left(\!\frac{1-\frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}-\left(1 + \frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}\right)e^{t \sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}}{1-\frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}\!+\!\left(1 + \frac{ \beta+\alpha-2\alpha z}{\sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}\right)e^{t \sqrt{(\beta+\alpha)^2 - 4\alpha\beta z}}}\right)\!\right)^{Q_0},\end{align}

where $Q_0$ is the active number in system at time 0.

In addition to calculating the probability generating function, we can also find a matrix calculation for the transient probability mass function of the counting process. To do so, we recognize that the time until the next arrival occurs can be treated as the time to absorption in a continuous-time Markov chain. By building from this idea to construct a transition matrix for several successive arrivals, we find the form for the distribution given in Proposition 4.

Proposition 4. Let $N_t$ be the number of arrivals by time t in an ESEP with baseline intensity $\eta^* > 0$ , intensity jump $\alpha > 0$ , and expiration rate $\beta > \alpha$ . Further, let $Q_0 = k$ be the initial active number in system. Then for $i \in \mathbb{N}$ , define the matrices $\textbf{D}_i \in \mathbb{R}^{k + i + 1 \times k + i + 1}$ and $\textbf{S}_i \in \mathbb{R}^{k+i+1 \times k + i+2}$ as

\begin{align*}\textbf{D}_i&=\begin{bmatrix}-(\eta^* + (k+i)(\alpha + \beta)) &\quad (k+i)\beta & & &\\ & \quad -(\eta^*+ (k+i-1)(\alpha + \beta)) & & &\\ & & \ddots &&\\ &&&\quad -(\eta^*+ \alpha+\beta)&\beta\\ &&&&\quad -\eta^*\end{bmatrix}\end{align*}

and

\begin{align*}\textbf{S}_i&=\begin{bmatrix}\eta^* + \alpha(k+i)\quad & & & & & \quad 0\\ & \quad\eta^* + \alpha(k+i-1) & & && \quad 0\\ & & \quad\ddots &&& \quad\vdots\\ &&&\quad\eta^* + \alpha&& \quad 0\\ &&&&\quad\eta^* &\quad 0\end{bmatrix}.\end{align*}

Further, let $\textbf{Z}_n \in \mathbb{R}^{\hat d_n \times \hat d_n}$ for $\hat d_n = \frac{n(n+1)}{2} + (n+1)(k+1)$ be a matrix such that

\begin{align*}\textbf{Z}_n=\begin{bmatrix}\textbf{D}_0 & \textbf{S}_0\quad &&&&\\&\quad \textbf{D}_1 & \textbf{S}_1 &&&\\&&\quad \ddots & \ddots &&\\&&& \quad\ddots & \quad \textbf{S}_{n-2} &\\&&&&\quad \textbf{D}_{n-1} & \quad \textbf{S}_{n-1}\\&&&&&\quad \textbf{D}_n\end{bmatrix}.\end{align*}

Then the probability that $N_t = n$ is given by

(7) \begin{align}\mathbb{P}\left( N_t = n \right)={{\textbf{v}}_{{1}}}^{\textrm{T}} e^{\textbf{Z}_n t} {{\textbf{v}}_{{:}}}\end{align}

where ${{\textbf{v}}_{{j}}} \in \mathbb{R}^{\hat d_n}$ is the unit column vector for the jth coordinate, and ${{\textbf{v}}_{{:}}} = \sum_{j = 0}^{k+n} {{\textbf{v}}_{{\hat d_n - j}}}$ .

With these fundamental quantities in hand, let us now turn to explore more nuanced connections between the ESEP and other stochastic processes in the following section. Doing so will provide further comparison between the Hawkes process and the ESEP, and moreover will formally connect the notion of self-excitement to similar concepts such as contagion, virality, and rich-get-richer effects.

3.1. Discrete-time perspectives through branching processes

Let us first view these processes through a discrete-time lens as branching processes. In this subsection we will interpret classical branching process results in application to these self-exciting processes. Taking the immigration–birth representation as inspiration, we start by considering the distribution of the total number of offspring of a single arrival. That is, we want to calculate the probability mass function for the number of arrivals that are generated directly from the excitement caused by the initial arrival. To constitute the total number of offspring, we will consider all the children of this initial entity across all time. For the ESEP, this equals the number of arrivals generated by the entity throughout its duration in the system; in the Hawkes process this counts the number of arrivals spurred by the entity as time goes to infinity. Given that the stability conditions are satisfied throughout, in Proposition 5 we calculate these distributions by way of inhomogeneous Poisson processes, yielding a Poisson mixture form for each.

Proposition 5. Let $\beta > \alpha > 0$ . Let $X^\eta$ be the number of new arrivals generated by the excitement caused by an arbitrary initial arrival throughout its duration in the system in an ESEP with jump size $\alpha$ and expiration rate $\beta$ . Then this offspring distribution is geometrically distributed with probability mass function

(8) \begin{align}\mathbb{P}\left( X^\eta = k \right) &= \left(\frac{\beta}{\alpha+\beta}\right)\left(\frac{\alpha}{\alpha+\beta}\right)^k.\end{align}

Similarly, let $X^\lambda$ be the number of new arrivals generated by the excitement caused by an arbitrary initial arrival in a Hawkes process with jump size $\alpha$ and decay rate $\beta$ . This offspring distribution is then Poisson distributed with probability mass function

(9) \begin{align}\mathbb{P}\left( X^\lambda = k \right) &= \frac{e^{-\frac{\alpha}\beta}}{k!}\left(\frac{\alpha}\beta\right)^k,\end{align}

where all $k \in \mathbb{N}$ .

Proof. Without loss of generality, we assume that the initial arrival in each process occurred at time 0. Then, at time $t \geq 0$ , the excitement generated by these initial arrivals has intensities given by $\alpha e^{-\beta t}$ and $\alpha \textbf{1}\{t < S\}$ for the Hawkes and ESEP processes, respectively, where $S \sim \textrm{Exp}(\beta)$ . Using [Reference Daley and Vere-Jones16], one can note that the offspring distributions across all time can then be expressed as

\begin{align*}X^\lambda \sim \textrm{Pois}\left(\alpha \int_0^\infty e^{-\beta t} \textrm{d}t\right)\qquad \text{and}\qquad X^\eta \sim \textrm{Pois}\left(\alpha \int_0^\infty \textbf{1}\{t < S\} \textrm{d}t\right),\end{align*}

which are equivalently stated as $X^\lambda \sim \textrm{Pois}\left(\alpha/\beta\right)$ and $X^\eta \sim \textrm{Pois}\left(\alpha S\right)$ . This now immediately yields the stated distributions for $X^\lambda$ and $X^\eta$ , as the Poisson–exponential mixture is known to yield a geometric distribution; see for example the overview of Poisson mixtures in [Reference Karlis and Xekalaki39].

We now move towards considering the total progeny of an initial arrival, meaning the total number of arrivals generated by the excitement of an initial arrival and the excitement of its offspring, and of their offspring, and so on across all time. It is important to note that in contrast to the number of offspring, the progeny includes the initial arrival itself. As we will see, the stability of the self-exciting processes implies that this total number of descendants is almost surely finite. This demonstrates the necessity of immigration for these processes to survive. From the offspring distributions in Proposition 5, the Hawkes descendant process is a Poisson branching process, and similarly the ESEP is a geometric branching process. These are well-studied models in branching processes, so we have many results available to us. In fact, we now use a result for random walks with potentially multiple simultaneous steps forward to derive the progeny distributions for these two processes. This is through the well-known hitting time theorem, stated below in Lemma 1.

Lemma 1. (Hitting time theorem.) The total progeny Z of a branching process with descendant distribution equivalent to $X_1$ is

\begin{align*}\mathbb{P}\left( Z = k \right) = \frac{1}{k}\mathbb{P}\left( X_1 + X_2 + \dots + X_k = k -1 \right),\end{align*}

where $X_1, \dots, X_k$ are i.i.d. for all $k \in \mathbb{Z}^+$ .

We now use the hitting time theorem to give the progeny distributions for the Hawkes process and the ESEP in Proposition 6. This is a common technique for branching processes, and it now yields valuable insight into these two self-exciting models.

Proposition 6. Let $\beta > \alpha > 0$ . Let $Z^\eta$ be a random variable for the total progeny of an arbitrary arrival in an ESEP with intensity jump $\alpha$ and expiration rate $\beta$ . Likewise, let $Z^\lambda$ be a random variable for the total progeny of an arbitrary arrival in a Hawkes process with intensity jump $\alpha$ and decay rate $\beta$ . Then the probability mass functions for $Z^\eta$ and $Z^\lambda$ are given by

(10) \begin{align}&\mathbb{P}\left( Z^\eta = k \right)=\frac{1}{k} \left(\begin{array}{c}2k-2\\k-1\end{array}\right) \left(\frac{\beta}{\beta + \alpha}\right)^k\left(\frac{\alpha}{\beta + \alpha}\right)^{k-1}\text{and}&\mathbb{P}\left( Z^\lambda = k \right)=\frac{e^{-\frac{\alpha}{\beta} k}}{k!}\left(\frac{\alpha k}{\beta}\right)^{k-1},\end{align}

where $k \in \mathbb{Z}^+$ .

Proof. This follows from applying Lemma 1 to Proposition 5. Because the sum of independent Poisson random variables is Poisson distributed with the sum of the rates, we have that

\begin{align*}\frac{1}{k}\mathbb{P}\left( X^\lambda_{1} + X^\lambda_{2} + \dots X^\lambda_{k} = k -1 \right)=\frac{1}k \mathbb{P}\left( K_1 = k -1 \right),\end{align*}

where $K_1 \sim \textrm{Pois}\left(\frac{\alpha k}{\beta}\right)$ . This now yields the expression for the probability mass function for $Z^\lambda$ . Similarly for $Z_\eta$ we note that the sum of independent geometric random variables has a negative binomial distribution, which implies that

\begin{align*}\frac{1}{k}\mathbb{P}\left( X^\eta_{1} + X^\eta_{2} + \dots X^\eta_{k} = k -1 \right)=\frac{1}k \mathbb{P}\left( K_2 = k -1 \right),\end{align*}

where $K_2 \sim \textrm{NegBin}\left(k,\frac{\alpha}{\beta + \alpha}\right)$ , and this completes the proof.

For a visual comparison of the descendants in the ESEP and the Hawkes process, we plot these two progeny distributions for equivalent parameters in Figure 2. As suggested by the variance ordering in Proposition 2, the tail of the ESEP progeny distribution is heavier than that of the Hawkes process.

Figure 2. Progeny distributions for the ESEP and the Hawkes process with matching parameters.

We can note that while one can calculate the mean of each progeny distribution via the probability mass functions in Proposition 6, it can also easily be found using Wald’s identity. Through standard infinitesimal generator approaches, one can calculate that the expected number of arrivals (including by immigration) in the ESEP is

\begin{align*}{\mathbb{E}\!\left[N_t\right]} = \frac{\beta \eta^* t}{\beta - \alpha} + \frac{\eta_0 - \eta_\infty}{\beta - \alpha}(1-e^{-(\beta - \alpha)t}).\end{align*}

However, using these branching process representations, we can also express this as

\begin{align*}{\mathbb{E}\!\left[N_t\right]} = {\mathbb{E}\!\left[\sum_{i=1}^{M_t} Z_{i}(t)\right]},\end{align*}

where $M_t$ is a Poisson process with rate $\eta^*$ and $Z_i(t)$ are the total progeny up to time $t \geq 0$ that are descended from the ith immigrant arrival. Now, by applying Wald’s identity to the limit of $\frac{1}{t}{\mathbb{E}\!\left[N_t\right]}$ as $t \to \infty$ , we see that

\begin{align*}\frac{\beta \eta^* }{\beta - \alpha} = \lim_{t \to \infty} \frac{{\mathbb{E}\!\left[N_t\right]}}{t} = \lim_{t\to\infty}\frac{1}{t}{\mathbb{E}\!\left[\sum_{i=1}^{M_t} Z_{i}(t)\right]} = \eta^* {\mathbb{E}\!\left[Z^\eta\right]} ,\end{align*}

and so ${\mathbb{E}\!\left[Z^\eta\right]} = \frac{\beta}{\beta - \alpha}$ . By analogous arguments for the Hawkes process, we see that ${\mathbb{E}\!\left[Z^\lambda\right]} = \frac{\beta}{\beta-\alpha}$ .

As a final branching process comparison between these two processes, we calculate the distribution of the total number of generations of descendants of an initial arrival in the ESEP and the Hawkes process. That is, let the first entity be the first generation, its offspring the second generation, their offspring the third, and so on. In Proposition 7 we find the probability mass function for the ESEP in closed form and a recurrence relation for the cumulative distribution function for the Hawkes process.

Proposition 7. Let $\beta > \alpha > 0$ . Let $\mathcal{G}^\eta$ be the number of distinct arrival generations across the full progeny of an initial arrival in an ESEP with intensity jump $\alpha$ and service rate $\beta$ . Then the probability mass function for $\mathcal{G}^\eta$ is given by

(11) \begin{align}\mathbb{P}\left( \mathcal{G}^\eta = k \right)=\frac{\alpha^{k-1}(\beta-\alpha)}{\beta^k-\alpha^k}-\frac{\alpha^k(\beta-\alpha)}{\beta^{k+1} - \alpha^{k+1}}.\end{align}

Likewise, let $\mathcal{G}^\lambda$ be the number of distinct arrival generations in the full progeny of an initial arrival for a Hawkes process with intensity jump $\alpha > 0$ and decay rate $\beta$ . Then $\mathcal{G}^\lambda$ has cumulative distribution function $F_{\mathcal{G^\lambda}}(k) = \mathbb{P}\left( \mathcal{G}^\lambda \leq k \right)$ satisfying the recursion

(12) \begin{align}F_{\mathcal{G^\lambda}}(k) = e^{-\frac{\alpha}{\beta}\left(1 - F_{\mathcal{G^\lambda}}(k-1)\right)} ,\end{align}

where $F_{\mathcal{G^\lambda}}(0) = 0$ and all $k \in \mathbb{Z}^+$ .

Proof. Let $Y_k^\lambda$ and $Y_k^\eta$ be Galton–Watson branching processes defined as

(13) \begin{align}&Y_k^\lambda = \sum_{i=1}^{Y_{k-1}^\lambda} X_{\lambda, i}^{(k)},&Y_k^\eta = \sum_{i=1}^{Y_{k-1}^\eta} X_{\eta, i}^{(k)},\end{align}

with $X_{\lambda, i}^{(k)} \stackrel{i.i.d.}{\sim} \textrm{Pois}\left(\frac{\alpha}{\beta}\right)$ , $X_{\eta, i}^{(k)} \stackrel{i.i.d.}{\sim} \textrm{Geo}\left(\frac{\alpha}{\alpha+\beta}\right)$ , and $Y_0^\lambda = Y_0^\eta = 1$ . These processes then have probability generating functions

\begin{align*}\mathcal{P}_k^\lambda(z)=\sum_{j=0}^\infty z^j\, \mathbb{P}\left( Y_k^\lambda = j \right)\quad\text{ and }\quad\mathcal{P}_k^\eta(z)=\sum_{j=0}^\infty z^j\, \mathbb{P}\left( Y_k^\eta = j \right)\end{align*}

that are given by the recursions $\mathcal{P}_{k+1}^\lambda(z) = \mathcal{P}_{X^\lambda}\left(\mathcal{P}_{k}^\lambda(z)\right)$ and $\mathcal{P}_{k+1}^\lambda(z) = \mathcal{P}_{X^\eta}\left(\mathcal{P}_{k}^\eta(z)\right)$ with $\mathcal{P}_{1}^\lambda(z) = \mathcal{P}_{X^\lambda}(z)$ and $\mathcal{P}_{1}^\eta(z) = \mathcal{P}_{X^\eta}(z)$ , where $\mathcal{P}_{X^\lambda}(z)$ and $\mathcal{P}_{X^\eta}(z)$ are the probability generating functions of $X_{\lambda,1}^{(1)}$ and $X_{\eta,1}^{(1)}$ , respectively; see e.g. Section XII.5 of [Reference Feller27]. One can then use induction to observe that

\begin{align*}\mathcal{P}_k^\eta(z)=1 - \frac{\alpha^k(1-z)}{\beta^k + \sum_{j=1}^k\alpha^j \beta^{k-j}(1-z)},\end{align*}

whereas $\mathcal{P}_k^\lambda(z) = e^{-\frac{\alpha}{\beta}\left(1 - \mathcal{P}^\lambda_{k-1}(z)\right)}$ , with $\mathcal{P}^\lambda_1(z) = e^{-\frac{\alpha}{\beta}(1-z)}$ . Because of their shared offspring distribution constructions, the number of the progeny in the kth arrival generations of the Hawkes process and the ESEP are equivalent in distribution to $Y_k^\lambda$ and $Y_k^\eta$ , respectively. In this way, we can express $\mathcal{G}^\lambda$ and $\mathcal{G}^\eta$ as

\begin{align*}\mathcal{G}^\lambda=\inf\{ k \in \mathbb{Z}^+ \mid Y_k^\lambda = 0 \}\quad\text{ and }\quad\mathcal{G}^\eta=\inf\{ k \in \mathbb{Z}^+ \mid Y_k^\eta = 0 \}.\end{align*}

This leads us to observe that the events $\{\mathcal{G}^\lambda = j\}$ and $\{Y_j^\lambda = 0, Y_{j-1}^\lambda > 0\}$ are equivalent, as are $\{\mathcal{G}^\eta = j\}$ and $\{Y_j^\eta = 0, Y_{j-1}^\eta > 0\}$ . Focusing for now on $\mathcal{G}^\lambda$ , we have that

\begin{align*}\mathbb{P}\left( Y_j^\lambda = 0, Y_{j-1}^\lambda > 0 \right)&=\sum_{i=1}^\infty \mathbb{P}\left( X_{\lambda,1}^{(1)} = 0 \right)^i\mathbb{P}\left( Y_{j-1}^\lambda = i \right)\\&=\mathcal{P}^\lambda_{j-1}\left(\mathbb{P}\left( X_{\lambda,1}^{(1)} = 0 \right)\right)-\mathbb{P}\left( Y_{j-1}^\lambda = 0 \right),\end{align*}

and since $\mathbb{P}\left( K = 0 \right) = \mathcal{P}(0)$ for any nonnegative discrete random variable K with probability generating function $\mathcal{P}(z)$ , this yields

\begin{align*}\mathbb{P}\left( \mathcal{G}^\lambda = j \right) = \mathcal{P}^\lambda_{j}(0) - \mathcal{P}^\lambda_{j-1}(0) .\end{align*}

Using $\mathcal{P}^\lambda_{0}(0) = 0$ , this telescoping sum now produces the stated form of the cumulative distribution function for $\mathcal{G}^\lambda$ . By analogous arguments for $\mathcal{G}^\eta$ , we complete the proof.

In the following subsection we focus on the ESEP, using the insight we have now gained from branching processes to connect this process to stochastic models for preferential attachment that are popular in the Bayesian nonparametric and machine learning literatures.

3.2. Similarities with preferential attachment and Bayesian statistics models

In the branching process perspective of the ESEP, consider the total number of active families at one point in time. That is, across all the entities present in the system at a given time, we are interested in the number of distinct families to which these entities belong. As each arrival occurs, the new entity either belongs to one of the existing families, meaning that the entity is a descendant, or it forms a new family, which is to say that it is an immigrant. If the entity is joining an existing family, it is more likely to join families that have more presently active family members.

We can note that these dynamics are quite similar to the definition of the Chinese restaurant process (CRP); see Chapter 11 in [Reference Aldous3]. The CRP models the successive arrivals of customers to a restaurant with infinitely many tables that each have infinitely many seats. Each arriving customer chooses which table to join based on the decisions of those before. Specifically, the nth customer to arrive joins table i with probability $\frac{c_i}{n-1+\lambda}$ or otherwise starts a new table with probability $\frac{\lambda}{n-1+\lambda}$ , where $c_i$ is the number of customers at table i and $\lambda > 0$ . As the number seated at table i grows larger, it is increasingly likely that the next customer will choose to sit at table i. In the ESEP, a new arrival at time $t \geq 0$ is generated as part of active excitement family i with probability $\frac{\alpha Q_{t,i}}{\alpha Q_t + \eta^*}$ and otherwise is an externally generated arrival with probability $\frac{\eta^*}{\alpha Q_t + \eta^*}$ , where $Q_{t,i}$ is the number of active exciters in the system at time t in the ith excitement family with $Q_t = \sum_{i} Q_{t,i}$ . By normalizing the numerator and denominator of these probabilities by $\frac{1}\alpha$ , we see that these dynamics match the CRP almost exactly. The difference is hardly a novel idea for restaurants—in the ESEP diners eventually leave. This departure then decreases the number of customers at the table, making it less attractive to the next person to arrive.

In addition to being an intriguing stochastic model, the CRP is also of interest for Bayesian statistics and machine learning through its connection to Bayesian nonparametric mixture models, specifically Dirichlet process mixtures. By consequence, the CRP then also has commonality with urn models and models for preferential attachment; see e.g. [Reference Blackwell and MacQueen10]. The CRP is also established enough to have its own generalizations, such as the distance-dependent CRP in [Reference Blei and Frazier12], in which the probability that a customer joins a table is dependent on a distance metric, and the recurrent CRP in [Reference Ahmed and Xing1], in which the restaurant closes at the end of each day, forcing all of that day’s customers to simultaneously depart. Drawing inspiration from the CRP and from the branching process perspectives of the ESEP, we investigate the distribution of the number of active families in the ESEP. This is equivalently stated as the number of active tables in a continuous-time CRP in which customers leave after their exponentially distributed meal durations. To begin, we first find the expected amount of time until a newly formed table becomes empty.

Proposition 8. Suppose that an ESEP receives an initial arrival at time 0. Let $X_t$ be the number of entities in the system at time $t \geq 0$ that are progeny of the initial arrival, and let $\tau$ be a stopping time such that $\tau = \inf \{t \geq 0 \mid X_t = 0\}$ . Then the expected value of $\tau$ is

(14) \begin{align}{\mathbb{E}\!\left[\tau\right]}=\frac{1}{\alpha}\log\left(\frac{\beta}{\beta - \alpha}\right),\end{align}

where $\alpha > 0$ is the intensity jump size and $\beta > \alpha$ is the expiration rate.

Proof. To observe this, we note that $X_t$ can be viewed as the state of an absorbing continuous-time Markov chain on the nonnegative integers. State 0 is the single absorbing state; in any other state j, the two possible transitions are to $j+1$ at rate $\alpha j$ and to $j-1$ at rate $\mu j$ , as visualized below.

Then, $\tau$ is the time of absorption into state 0 when starting in state 1, and so ${\mathbb{E}\!\left[\tau\right]}$ can be calculated by standard first-step analysis approaches, yielding

\begin{align*}{\mathbb{E}\!\left[\tau\right]}=\sum_{i=1}^\infty \frac{1}{\alpha i} \prod_{j=1}^{i} \frac{\alpha j}{\beta j} = \frac{1}{\alpha}\sum_{i=1}^\infty \frac{1}{i}\left(\frac{\alpha}{\beta}\right)^i = \frac{1}{\alpha}\log\left(\frac{1}{1-\frac{\alpha}{\beta}}\right),\end{align*}

and this simplifies to the stated result.

Proposition 8 gives the expectation of the total time an excitement family is active in the system. Using this, in Proposition 9 we now employ a classical queueing theory result to find the exact distribution of the number of active families simultaneously in the system in steady state.

An interesting consequence of the number of active families being Poisson distributed and the total number in system being negative-binomially distributed is that it suggests that the number of simultaneously active family members is logarithmically distributed. We observe this via the known compound Poisson representation of the negative binomial distribution [Reference Willmot60]. For $B \sim \textrm{Pois}\left(\frac{\eta^*}{\alpha}\log \left(\frac{\beta}{\beta - \alpha}\right)\right)$ , $Q \sim \textrm{NegBin}\left(\frac{\alpha}{\beta}, \frac{\eta^*}{\alpha}\right)$ , and $L_i \stackrel{\textrm{i.i.d.}}{\sim} \textrm{Log}\left(\frac{\alpha}{\beta} \right)$ , one can observe that

\begin{align*}Q \stackrel{D}{=} \sum_{i=1}^B L_i,\end{align*}

where $\mathbb{P}\left( L_1 = k \right) = \left(\frac{\alpha}{\beta}\right)^k\left(k \log\left(\frac{\beta}{\beta-\alpha}\right)\right)^{-1}$ for all $k \in \mathbb{Z}^+$ . Thus, the idea that the number of active members of each family is logarithmically distributed follows from the fact that this is a sum of positive integer-valued random variables, of which there are as many as there are active families, and this sum is equal to the total number in system.

3.3. Connections to epidemic models

As a final observation regarding the ESEP and its connections to other stochastic models, consider disease spread. As we discussed in the introduction to this paper, when a person becomes sick with a contagious disease she increases the rate of new infections through her contact with others. Furthermore, when a person recovers from a disease such as the flu, she is no longer contagious and thus no longer contributes to the rate of disease spread. While we have discussed in the introduction that this scenario has the hallmarks of self-excitement qualitatively, a classic model for studying this phenomenon is the susceptible–infected–susceptible (SIS) process.

In the SIS model there is a finite population of $N \in \mathbb{Z}^+$ individuals. Each individual takes on one of two states, either infected or susceptible. Let $I_t$ be the number infected at time $t \geq 0$ , and let $S_t$ be the number susceptible. In the continuous-time stochastic SIS model, each infected individual recovers after an exponentially distributed duration of the illness. Once a person recovers from the disease, she becomes susceptible again. Because there is a finite population, the rate of new infections depends on both the number infected and the number susceptible; a new person falls ill at a rate proportional to $ I_t\cdot \frac{S_t}{N}$ . Because this continuous-time Markov chain would be absorbed into state $I_t = 0$ , it is common to include an exogenous infection rate proportional to just $\frac{S_t}{N}$ . We will refer to this model as the stochastic SIS with exogenous infections; Figure 3 shows the rate diagram for the transitions from infected to susceptible and from susceptible to infected. For the sake of comparison, we set the exogenous infection rate as $\eta^*$ , the epidemic infection rate as $\alpha$ , and the recovery rate as $\beta$ .

Figure 3. Stochastic SIS model with exogenous infections.

One can note that there are immediate similarities between this process and the ESEP. That is, new infections increase the infection rate while recoveries decrease it, and infections can be the result of either external or internal stimuli. However, the primary difference between these two models is that the SIS process has a finite population, whereas the ESEP does not. In Proposition 10 we find that as this population size grows large the difference between these models fades, yielding that the distribution of the number infected in the exogenously driven SIS model converges to the distribution of the queue length in the ESEP.

Proposition 10. Let $I_t$ be the number of infected individuals at time $t \geq 0$ in an exogenously driven stochastic SIS model with population size $N \in \mathbb{Z}^+$ , exogenous infection rate $\eta^* > 0$ , epidemic infection rate $\alpha > 0$ , and recovery rate $\beta > 0$ . Then, as $N \to \infty$ ,

\begin{align*}I_t\stackrel{D}{\Longrightarrow}Q_t,\end{align*}

where $Q_t$ is the active number in system at time t for an ESEP with baseline intensity $\eta^*$ , intensity jump $\alpha$ , and expiration rate $\beta$ .

Proof. Because the SIS model is a Markov process, one can use the infinitesimal generator approach to find a time derivative for the moment generating function of the number of infected individuals at time $t \geq 0$ . Thus, by noting that $S_t = N - I_t$ we have that

\begin{align*}\frac{\textrm{d}}{\textrm{d}t}{\mathbb{E}\!\left[e^{\theta I_t}\right]}&={\mathbb{E}\!\left[\frac{\alpha I_t S_t}{N} \left( e^{\theta}-1\right) e^{\theta I_t}+\beta I \left(e^{-\theta} - 1 \right) e^{\theta I_t}+\frac{\eta^* S_t}{N} \left(e^{\theta} - 1\right) e^{\theta I_t}\right]}\\&={\mathbb{E}\!\left[\frac{\alpha I_t (N - I_t)}{N} \left( e^{\theta}-1\right) e^{\theta I_t}\right]}+{\mathbb{E}\!\left[\beta I_t \left(e^{-\theta} - 1 \right) e^{\theta I_t}\right]}\\&\qquad +{\mathbb{E}\!\left[\frac{\eta^* (N- I_t)}{N} \left(e^{\theta} - 1\right) e^{\theta I_t}\right]},\end{align*}

which we can re-express in PDE form as

\begin{align*}\frac{\partial {\mathbb{E}\!\left[e^{\theta I_t}\right]}}{\partial t}&=\left(\alpha \left( e^{\theta}-1\right)+\beta\! \left(e^{-\theta} - 1 \right)-\frac{\eta^* }{N} \left(e^{\theta} - 1\right)\right)\frac{\partial {\mathbb{E}\!\left[e^{\theta I_t}\right]}}{\partial \theta}-\frac{\alpha }{N} \left( e^{\theta}-1\right)\frac{\partial^2 {\mathbb{E}\!\left[e^{\theta I_t}\right]}}{\partial \theta^2}\\&\quad+\eta^* \left(e^{\theta} - 1\right){\mathbb{E}\!\left[e^{\theta I_t}\right]}.\end{align*}

Now as the population size $N \to \infty$ , this converges to

\begin{align*}\frac{\partial {\mathbb{E}\!\left[e^{\theta I_t}\right]}}{\partial t}&=\left(\alpha \left( e^{\theta}-1\right)+\beta\! \left(e^{-\theta} - 1 \right)\right)\frac{\partial {\mathbb{E}\!\left[e^{\theta I_t}\right]}}{\partial \theta}+\eta^* \left(e^{\theta} - 1\right){\mathbb{E}\!\left[e^{\theta I_t}\right]},\end{align*}

which we can recognize as the PDE for the moment generating function of the ESEP through its own infinitesimal generator.

As a demonstration of this convergence, we plot the empirical steady-state distribution of the SIS process for increasing population size below in Figure 4. Note that in this example the distributions appear fairly close for populations of size $N = 1,000$ . On the scale of the populations of cities or even some larger high schools, this is quite small. At a medium university size of $N = 10,000$ , the distributions are essentially indistinguishable.

Figure 4. Steady-state distribution of the number infected in the exogenously driven SIS model for increasing population size N, where $\eta^* = 10$ , $\alpha = 2$ , and $\beta = 3$ .

We would be remiss if we did not note that connections from epidemic models to birth–death processes are not new. For example, [Reference Ball8] demonstrated that epidemic models converge to birth–death processes, and [Reference Singh and Myers56] even noted that the exogenously driven susceptible–infected–recovered (SIR) model—in which people cannot become reinfected—converges to a linear birth–death–immigration process; however, these works did not outright form connections to self-exciting processes. In [Reference Rizoiu54], the similarities between the Hawkes process and the SIR process are shown and formal connections are made, although this is through a generalization of the Hawkes process defined on a finite population rather than through increasing the epidemic model population size. Regardless, these prior works serve to expand the practical relevance of the ESEP, as they show that such epidemic models are also of use outside of public health. For example, these models have been used to study the phenomenon of contagion in areas such as product adoption, idea spread, and social influence. These all also naturally relate to the concept of self-excitement, and in Proposition 10 we observe that this connection can be formalized.

These connections take on added importance in the contemporary context of the COVID-19 public health crisis. It is worth noting that the convergence of distributions in Proposition 10 does not require our commonly assumed stability condition $\beta > \alpha$ . Thus, this convergence also covers potential pandemic scenarios in which $\alpha \geq \beta$ . In such settings, one can quickly observe that the mean number infected by time t is such that ${\mathbb{E}\!\left[N_t\right]} \in O(e^{(\alpha - \beta)t})$ , reproducing the exponential growth exhibited by the novel coronavirus in these early stages. The ESEP arrival process then captures the times of new infections, some portion of which may then be used as the arrival process to a queueing model for the cases that require hospitalization.