A Monte Carlo EM Algorithm for the Parameter Estimation of Aggregated Hawkes Processes

Point processes are extensively used to model event data and have found wide applications in many fields including seismology [23] and cyber-security [24]. One representation of a point process is via the counting process $\{N(t),\,t\in\mathbb{R}\}$, where $N(t)$ denotes the number of events up to time $t$, and $N(0)=0$. Due to limited recording capabilities and storage capacities, retaining event times with a high precision is expensive and often infeasible. Therefore, in much real-world data, it is common to instead observe a times series of the aggregated process,

for some $\Delta>0$, which we refer to as the bin width. This aggregated process may arise from a predetermined binning of the data into counts, or equivalently from the rounding of event times. In the context of network traffic data, for example, the resolution of the recorded times can be anywhere from milliseconds to seconds (as is the case with the Los Alamos National Laboratory (LANL) NetFlow data [26]), or even coarser. In this setting the value of this aggregated process at each time point is the number of events with that rounded time-stamp.

Intuitively, when aggregating data we lose information and essentially ‘blur’ our view of the continuous time point process, making it potentially problematic to apply methods which assume a continuous time framework. Thus, the problem we consider here is to infer upon the underlying continuous process from the observed aggregated data.

The Hawkes process is a type of ‘self-exciting’ process which provides us with a model for contagious event data. Their flexibility and real-world relevancy has resulted in a host of applications. In the case of financial data for example, this allows propagation of stock crashes and surges to be modeled [1, 3, 8, 9, 7]. Propagation of social media events has also been modeled using Hawkes processes, in particular ‘twitter cascades’ are considered in [25, 15]. Further applications include the modeling of civilian deaths due to insurgent activity in Iraq [18], and predicting origin times and magnitudes of earthquakes [22].

Formally, the Hawkes process is a class of stochastic process with the defining property that

where ${\rm{d}}N(t)=N(t+{\rm{d}}t)-N(t)$. It is characterized via its conditional intensity function (CIF) $\lambda^{\ast}(t)$, defined as

where $\nu$ is called the background intensity and $g(u)$ is the excitation kernel. This means the intensity at an arbitrary time-point is dependent on the history of the process, producing self-exciting behavior. Depending on the kernel $g(u)$, the excitation may be quite local, or have longer term effects [10].

We can also consider a Hawkes process as a branching process of time-stamped events. From this viewpoint, formalised in [11], events can be seen to arrive either via immigration or birth. That is, an event can be triggered by the background intensity rate $\nu$, in which case the event is seen as an immigrant. Alternatively, an event which is cause by self-excitation can be considered a descendant, referred to as being generated ‘endogenously’. Unlike a homogeneous Poisson point process, where events happen independently and at a constant rate, self-exciting processes are such that there is a higher likelihood of events happening in the near future of an arbitrary event and this is due to endogenous triggering [25]. In this way we can consider the branching ratio of a Hawkes process, defined as

The inequality above ensures that the process does not ‘explode’, a case in which we have infinite events occurring in some finite time interval [16].

As introduced in Hawkes’ original paper, exponential decay is a common choice for the excitation kernel due to the simplifications it provides for the theoretical derivations [10, 16]. In this case we can write the excitation kernel as a sum of $L$ exponential decays,

Here, and as is most common, we let $L=1$ when considering the exponential kernel. In this case, the branching ratio defined in (1) becomes

Given this model, we wish to estimate the parameter set $\Theta=\{\nu,\alpha,\beta\}$. Other kernels can be used, including a power-law function of form $g(u)=\alpha\beta{(1+\beta u)^{-(1+c)}}\mathds{1}_{\mathbb{R}_{+}}(u)$, in which case $\Theta=\{\nu,\alpha,\beta,c\}$. In the continuous time setting, parameter estimation for any of these kernels is straightforward.

The EM algorithm [6] is an iterative method for the computation of the maximizer of a likelihood. The idea of this algorithm is to augment the observed data by a latent quantity [27]. In the case considered here, the observed data are the event counts per unit time. We denote this by ${\mathbf{N}}=[N_{1},\ldots,N_{K}]$, where $N_{j}$ denotes the counts in the $j^{\rm{th}}$ bin ($j=1,...,K$). The latent data, denoted $\mathcal{T}$ are the unobserved, true event times which are rounded on recording and the set of parameters to be estimated is denoted $\Theta=\{v,\alpha,\beta\}$. The algorithm proceeds as follows:

1. 1.

   In the E (Expectation) step, we evaluate

   $\displaystyle Q_{i+1}(\Theta,\Theta^{i})=\int_{\mathbb{T}^{\ast}}\log(p(\Theta\mid{\mathbf{N}},\mathcal{T}))p(\mathcal{T}\mid{\mathbf{N}},\Theta^{i}){\rm{d}}\mathcal{T},$

   (7)

   where $\mathbb{T}^{\ast}$ denotes the sample space for $\mathcal{T}$. That is, we compute the expectation $Q_{i+1}(\Theta,\Theta^{i})$ of the log-posterior $\log(p(\Theta\mid{\mathbf{N}},\mathcal{T}))$ with respect to the conditional predictive distribution $p(\mathcal{T}\mid{\mathbf{N}},\Theta^{i})$, where $\Theta^{i}$ is the current, $i^{\rm th}$ approximation.

2. 2.

   In the M (Maximization) step, we update the value of the conditional expectation with its maximizer $\Theta^{i+1}$.

When (7) is analytically intractable we require Monte Carlo methods for numerical computation. This is known as MC-EM [27]. If we are able to sample $\mathcal{T}$ directly from $p(\mathcal{T}\mid{\mathbf{N}},\Theta^{i})$, then we can approximate the integral in (7) with

where $\mathcal{T}^{\ast(k)}$ is the $k^{\rm{th}}$ Monte Carlo sample of $\mathcal{T}$. However, no such sampling regime is possible in the Hawkes process setting. We therefore use importance sampling to simulate a legal proposal for $\mathcal{T}$ (that is, a set of event times that match the binned counts) from an alternative distribution $q(\mathcal{T}\mid{\mathbf{N}},\Theta^{i})$ which is simple to sample from (see Section 2.1 for details). Each of these proposals is then weighted depending on the probability it came from the desired distribution. That is, given a set of $m$ samples $\mathcal{T}^{\ast(1)},\ldots,\mathcal{T}^{\ast(m)}$, we assign weights

and approximate (7) with

We note that the numerator of (8) can be expressed as

where $\sum_{i=1}^{n}\mathds{1}_{[b_{i}^{-},b_{i}^{+})}\left(\mathcal{T}_{i}^{\ast(k)}\right)$ is the number of event times in $\mathcal{T}^{\ast(k)}$ lying in bin $j$, and $\mathds{1}_{B}(y)$ is 1 if $y$ is in the set $B$ and 0 otherwise. Therefore, if only proposing legal event times, we have that

where, $\log\left(p(\mathcal{T}^{\ast(k)}\mid\Theta)\right)$ is given by (5).

However, the question remains of how to best sample the latent times. Here, the density we would ideally like to sample from is that of the missing event times given the bin counts and the model parameters $p(\mathcal{T}\mid{\mathbf{N}},\Theta^{i})$. Therefore for this method to be most efficient and to ensure meaningful weights, the alternative distribution $q(\mathcal{T}\mid{\mathbf{N}},\Theta^{i})$ should be as close as possible to the true distribution of the time-stamps.

Given parameters $\nu,\alpha,\beta$ and some maximum simulation time $T$, we can simulate realizations of a Hawkes process. The generated events are those which form the latent space $\mathcal{T}$, and aggregating these to different $\Delta$ allows us to simulate the count data $\mathbf{N}$. We can then apply each of the three methods detailed: the binned log-likelihood, INAR($p$) approximation, and the MC-EM method. Figures 1-6 show boxplots for the estimates of each of $\nu,\alpha,$ and $\beta$ for 20 realizations of a Hawkes process with the ground truth parameters specified. The mean value of each parameter estimate is presented on the vertical axis. We clearly see that the INAR($p$) approximation method can yield highly variable results. In Figure 1, the boxplots for both the excitation rate, $\alpha$, and the decay rate $\beta$ have been presented on a log-scale in order to show the results on one axis. Clearly, the INAR($p$) method has resulted in outliers that are factors of 10 away from the ground truth. The remaining figures all likewise show the MC-EM method to perform better than either of the two alternative approaches considered for a range of different parameter sets. In particular, again in Figure 5 we note that the excitation rate estimates have been presented on a log-scale due to the extreme outliers in the INAR($p$) method. The binned log-likelihood method, whilst suffering less from extreme outliers, performs less well than the MC-EM method in all cases.

In Figure 7 we also consider the bias across different levels of aggregation. That is, for each of the realizations of a self-exciting process for a given parameter set, we can aggregate the data to different levels and compare the bias in the parameter estimates. The right hand plot in Figure 7 presents the bias on a log-scale. It is evident that the INAR($p$) method is more biased for larger bin widths. The binned log-likelihood performs better, however, still not as well as the MC-EM method which most consistently exhibits a low bias.

Here we presented a new technique for handling aggregated data using an MC-EM algorithm. By sampling from a distribution close to that of the latent event times given the observed times and current parameter estimate, we have proposed a surrogate legal set of candidate time-points. This allows estimation of parameters using methods for continuous time-points. We further compared this to the INAR($p$) approximation proposed in [12] and a binned log-likelihood method which assumes a piecewise constant CIF within each interval. For the parameter sets considered, the MC-EM method has appeared to out-perform both alternatives. The MC-EM approach also provides us with additional flexibility in that the level of aggregation does not need to be consistent across the dataset. That is, provided the interval bounds are known, $\Delta$ can vary. The issues arising from aggregated data could also be handled via a MCMC (MCMH) algorithm and exploring this is the subject of future work.

Leigh Shlomovich is funded by an Industrial Strategy Engineering and Physical Sciences Research Council Scholarship and Lekha Patel is funded by an Imperial College President’s Scholarship.