{'fraction_non_alphanumeric': 0.04867406796916662, 'fraction_numerical': 0.007902264022882629, 'mean_word_length': 4.222278481012658, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Monte Carlo EM Algorithm for the Parameter Estimation of Aggregated Hawkes Processes

A key difficulty that arises from real event data is imprecision in the recording of event time-stamps. In many cases, retaining event times with a high precision is expensive due to the sheer volume of activity. Combined with practical limits on the accuracy of measurements, aggregated data is common. In order to use point processes to model such event data, tools for handling parameter estimation are essential. Here we consider parameter estimation of the Hawkes process, a type of self-exciting point process that has found application in the modeling of financial stock markets, earthquakes and social media cascades. We develop a novel optimization approach to parameter estimation of aggregated Hawkes processes using a Monte Carlo Expectation-Maximization (MC-EM) algorithm. Through a detailed simulation study, we demonstrate that existing methods are capable of producing severely biased and highly variable parameter estimates and that our novel MC-EM method significantly outperforms them in all studied circumstances. These results highlight the importance of correct handling of aggregated data.


Introduction
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 ∈ 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, N t = N (∆(t + 1)) − N (t∆) , for some ∆ > 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,[7][8][9]. Propagation of social media events has also been modeled using Hawkes processes, in particular 'twitter cascades' are considered in [15,25]. 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 dN (t) = N (t + dt) − N (t). It is characterized via its conditional intensity function (CIF) λ * (t), defined as where ν 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 ν, 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 Θ = {ν, α, β}. Other kernels can be used, including a power-law function of form g(u) = αβ(1 + βu) −(1+c) 1 R+ (u), in which case Θ = {ν, α, β, c}. In the continuous time setting, parameter estimation for any of these kernels is straightforward.

Continuous Time Framework
Typically, maximum likelihood estimation (MLE) is used to estimate parameters of a point process from a set of exact event times T = {t 1 , . . . , t NT } ⊂ (0, T ]. The CIF, or hazard function, of a point process is formally defined as where f * (t) and F * (t) are the conditional PDF and CDF, respectively, of the next arrival time, given the history of the process. From (3) we have that where t N denotes the last observed time prior to t. From (4), the joint likelihood of the univariate observations over the window [0, T ] is Thus, from Proposition 7.2.III of [5], the log-likelihood is given by If specifically considering a Hawkes process with exponential excitation kernel of form α exp(−βt), this loglikelihood can be simplified and expressed recursively as shown in [16]. In this paper we consider methods that are required when we instead observe a binned sequence of event counts. An alternative but equivalent representation of this is a discretization or rounding of the latent time-stamps, but here we will consider the observed data as the aggregation of T to bins.

Aggregated Data
In the literature, the issue of aggregated data is handled in many ways, from uniformly redistributing events across the bin [3], to only retaining unique time-stamps and discarding the rest [19]. Here, we propose a novel Monte Carlo-Expectation Maximization (MC-EM) method and compare it to two existing approaches, evaluating the performance of parameter estimation for each. The methods compared are: 1. Approximating a binned Hawkes process as an integer-valued auto-regressive (INAR) process, a method developed in [12,13] and described in Section 1.3.
2. Formulating a binned log-likelihood which assumes a piecewise constant CIF within each interval as described in Section 1.4.

A novel MC-EM approach detailed in Section 2.
There are other methods which have been covered in the literature, but are not considered here due to lack of applicability to this problem. As an example, a significant amount of the literature which aims to work with aggregated data considers binning the time-points such that the process contains at most one event per bin as in [21]. This is inappropriate here as we do not have access to the latent event times and so cannot select an appropriate discretization level, ∆. Likewise, some literature considers modeling binned behavior as a Bernoulli process [4]. Again, this is invalid here as it fails to account for the number of events in a bin and thus will heavily bias results. There exist methods that handle missing data when we observe continuous time-points with gaps in the recording windows [17]. That is, when the data considered contains precise but intermittent recordings. This is a closely related issue, however differs in the fact that when handling aggregated data, we do not have any precise times to work with.
We now outline the two methods against which we will compare our novel MC-EM approach.

Hawkes INAR(p) Approximation
It is shown in [12,13] that the distribution of the bin-count sequence of Hawkes processes can be represented by an integer-valued autoregressive model, known as the INAR(p) model, further details of which can be found in [12]. By representing the binned Hawkes counting process as an INAR(p) process, a non-parametric estimator for kernel g(u) is then formulated in terms of conditional least squares (CLS).
Let ∆ > 0 be the bin width, the univariate Hawkes process bin-count sequence is denoted N = [N 1 , . . . , N K ], for K = ⌊T /∆⌋ and N i denotes the counts in the i th bin. Then, defining some support ∆ < s < T , the CLS-operator is used on the bin counts N, with maximal lag p = ⌈s/∆⌉. Thus, where the design matrix Z is given as and Y is the lagged bin-count sequence, being [N p+1 , . . . , N k ]. Then the entries ofĝ (∆,s) , are estimates for the excitation kernel at the corresponding time-points. Kernel parameters α and β are then estimated by fitting an exponential function to these points.
Simulation studies examining the effect of bin width ∆ and parameter s are presented in [12,14], where they determine ∆ to have the greatest bearing on the quality of the estimates. There are however two points to note with this method. Firstly, CLS requires the inversion of ZZ ⊤ . In this case, this matrix contains the event counts per bin, and so it is possible to have cases where this matrix is singular, in particular when the counting process is very sparse. Secondly, as it is currently presented, this method does not constrain the parameter estimates to be those of a stationary Hawkes process. Therefore it is possible to yield infeasible estimates.

Binned Likelihood
An alternative method, briefly mentioned in [20] and developed here, considers sampling λ * at each discrete timepoint j∆ (j = 1, ..., K with K = ⌊T /∆⌋), thus representing λ * as a piecewise constant function within each bin. Letting N j be the number of events occurring in the sampling interval ((j − 1)∆, j∆] and using (5) we have that the log-likelihood of the underlying Hawkes process is where λ * (j∆) ≡ λ * (j∆|H j∆ ) and H j∆ denotes the history of the process until time j∆.
The assumption of a piecewise constant CIF is equivalent to assuming N j ∼ Poisson(∆λ * ((j − 1)∆)). However, it is important to note that this assumption is not correct as it ignores the excitation within each bin and therefore will be biased, especially in cases where the intensity is high relative to the bin width, ∆. Nevertheless it provides us with a simple approximation. To estimate the process parameters we maximize (6) with constraints ensuring stationarity, as expressed by (1). In the case of an exponential excitation kernel, explicitly this implies ν, α > 0 and α/β < 1.
We will now propose an alternative method of parameter estimation which iteratively uses 'legal' sets of continuous candidate time-points and therefore does not assume a piecewise constant CIF.

Monte Carlo EM Algorithm for Aggregated Data
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 N = [N 1 , . . . , N K ], where N j denotes the counts in the j th bin (j = 1, ..., K). The latent data, denoted T are the unobserved, true event times which are rounded on recording and the set of parameters to be estimated is denoted Θ = {v, α, β}. The algorithm proceeds as follows: 1. In the E (Expectation) step, we evaluate where T * denotes the sample space for T . That is, we compute the expectation Q i+1 (Θ, Θ i ) of the logposterior log(p(Θ | N, T )) with respect to the conditional predictive distribution p(T | N, Θ i ), where Θ i is the current, i th approximation.
2. In the M (Maximization) step, we update the value of the conditional expectation with its maximizer Θ 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 T directly from p(T | N, Θ i ), then we can approximate the integral in (7) with where T * (k) is the k th Monte Carlo sample of T . However, no such sampling regime is possible in the Hawkes process setting. We therefore use importance sampling to simulate a legal proposal for T (that is, a set of event times that match the binned counts) from an alternative distribution q(T | N, Θ 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 T * (1) , . . . , T * (m) , we assign weights and approximate (7) with We note that the numerator of (8) can be expressed as is the number of event times in T * (k) lying in bin j, and 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 p(T * (k) | Θ) 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(T | N, Θ i ). Therefore for this method to be most efficient and to ensure meaningful weights, the alternative distribution q(T | N, Θ i ) should be as close as possible to the true distribution of the time-stamps.

Sampling Method
It is possible to uniformly redistribute the events across a bin in order to generate a legal set of event times. However, especially for Hawkes processes with high activity, this is not optimal as it leads to weights that are too small to compute (9). Therefore, we propose an alternative method which samples from a distribution q(T | N, Θ) that more closely matches p(T | N, Θ). Here, this is developed for the exponential kernel but it is easily extendible to other kernels (see Appendices A and B for details regarding power-law and rectangular kernels, respectively).
Without loss of generality, consider having already simulated until time-point t n , n ∈ {1, . . . , N T }. Then, suppose we know that there are m time-points in the next non-empty bin, being t n+1 , . . . , t n+m . The joint probability of these events can be expressed using factorization. That is For brevity, we refer to f * Tn+1,...,Tn+m (t n+1 , . . . , t n+m ) as f * Tn+1:n+m (t n+1:n+m ). Note that for the simplest case of an exponential kernel, we can express this as As we wish to simulate possible realizations of the events given observed counts, we should account for the the fact that each time-point is known to have occurred within a given interval. That is, we account for the observed interval range [b − , b + ] for events in a given bin by considering the truncated joint density. We require that b − < t n+1 < t n+2 < . . . < t n+m < b + . Therefore, we can express the conditional CDF over this region as Even in the simplest case of an exponential decay kernel, this appears intractable due to the form of the conditional intensity function for a Hawkes process. Therefore we truncate the PDF by considering the joint CDF. As with the joint PDF, we can use factorization to express the joint CDF as , where we can use the form given in (4) for the CDF of each successive time-point given the history of the process. That is, the joint CDF of time-points t n+1 , . . . , t n+m is given by In the case of an exponential decay kernel, this is Thus the joint truncated PDF of m time-points given the history, t 1 , . . . t n can be expressed as f * Tn+1,...,Tn+m (t n+1 , . . . , t n+m ) κ , where κ = F * Tn+1,...,Tn+m (t n+1 , t n+2 , . . . , t n+m−1 , b + ) − F * Tn+1,...,Tn+m (b − , t n+2 . . . , t n+m−1 , t n+m ). The set of m proposed time-stamps for the given bin {t n+1 , . . . ,t n+m } are those that maximize (10).
In this way we can sequentially simulate a continuous version of the observed aggregate Hawkes process by progressively handling each bin such that we jointly maximize this likelihood.
Given parameters ν, α, β and some maximum simulation time T , we can simulate realizations of a Hawkes process. The generated events are those which form the latent space T , and aggregating these to different ∆ allows us to simulate the count data 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 ν, α, and β 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, α, and the decay rate β 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.

Conclusion
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 loglikelihood 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, ∆ 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.

A Power-Law Kernel
The proposed MC-EM method can still be applied if intending to consider a regularized power-law kernel of the form for α, β > 0. In this case, to ensure a stationary Hawkes process, we have that Therefore the stationarity condition is met for α < c [2].
We also need to consider the proposed univariate sampling method and thus the complete log-likelihood of sampled time-points. Both of these points fundamentally rely on expressing the conditional PDF and CDF. Firstly, for the sampling method introduced in Section 2.1, we now have that Similarly, the joint conditional CDF is given by αβ (1 + β(t n+i − t j )) 1+c du .
Then, (10) gives the form for the truncated PDF, as previously. All that remains is to adjust the log-likelihood for the CIF with regularized power-law function when implementing the MC-EM algorithm. Using (5) that is, αβ (1 + β(u − t j )) 1+c du.

B Rectangular Kernel
We can also consider a rectangular kernel of the form for α, β > 0. In this case, stationarity holds if β α n β − α du = n < 1.