A parameter estimation method for multivariate binned Hawkes processes

It is often assumed that events cannot occur simultaneously when modelling data with point processes. This raises a problem as real-world data often contains synchronous observations due to aggregation or rounding, resulting from limitations on recording capabilities and the expense of storing high volumes of precise data. In order to gain a better understanding of the relationships between processes, we consider modelling the aggregated event data using multivariate Hawkes processes, which offer a description of mutually-exciting behaviour and have found wide applications in areas including seismology and finance. Here we generalise existing methodology on parameter estimation of univariate aggregated Hawkes processes to the multivariate case using a Monte Carlo expectation–maximization (MC-EM) algorithm and through a simulation study illustrate that alternative approaches to this problem can be severely biased, with the multivariate MC-EM method outperforming them in terms of MSE in all considered cases.


Introduction
Modern data acquisition systems in many applications collect vast amounts of data that can be characterised as time series. An area of particular interest is cyber-security, where network data are recorded and questions arise about the correlation structure and 'excitational' effects in and between the resulting time series. We aim to uncover trends within and between such time stamped event data, which we model using point processes. For this, we are particularly interested in the multivariate Hawkes process, which provides a model for 'mutually exciting' events. Introduced in Hawkes (1971a, b) primarily to model the occurrence of seismic activity, multivariate Hawkes processes have found wide application to many disciplines due to their ability to model crossexcitation. In the case of financial data, multivariate Hawkes processes have been used to model the joint dynamics of trades and mid-price changes of the NYSE (Bowsher 2007) and an overview of literature surrounding the application of these processes to finance is given in Bacry et al. (2015). Additionally Hawkes processes have been considered within cyber-security for modelling of computer network traffic (Mark et al. 2019; Price-Williams and Heard 2019) and social media activity, for example in Kobayashi and Lambiotte (2016) where propagation of 'Twitter cascades' is considered. As discussed in Shlomovich et al. (2022) practical limitations on storing or recording high resolution data results in an abundance of binned data, that is streams of counts of occurrences per time bin. We can use the counting process representation of a P-variate point process to model this binned data. Let N ( p) (t) be the pth continuous time count process ( p = 1, . . . , P) denoting the number of events until time t ∈ R where N ( p) (0) = 0, N ( p) (t) = N ( p) ((0, t]) for t > 0 and −N ( p) ((t, 0]) for t < 0. We further denote to be the binned (aggregated) process. Then the discrete-time vector process N t is given by where K = T /Δ and Δ is the bin width. Here we develop and extend the methodology detailed in Shlomovich et al. (2022) to handle parameter estimation of multivariate binned Hawkes processes. We do this by considering the superposition of the multivariate count process. This is an elegant method which manages to replicate both the marginal and the cross-correlation structure of the true process. Results from a simulation study, comparing the extended MC-EM algorithm to both the INAR( p) and binned log-likelihood approximations are presented in Sect. 3. This method is shown to produce estimates with lower MSE than currently available alternatives. We additionally derive the Hessian of a multivariate Hawkes process with exponential kernel, required for optimization, and present this in "Appendix B" along with the gradient derived in Ozaki (1979).

Multivariate Hawkes processes
Formally, the P-dimensional Hawkes process N(t) = N (1) (t), . . . , N (P) (t) is a class of stochastic process such that independently for each p ∈ {1, . . . P}, (Hawkes 1971a). It is characterized via its conditional intensity function (CIF) λ * ( p) (t), defined as where ν > 0 is a P-dimensional vector called the background intensity and g(u) is the non-negative excitation kernel such that g(u) = 0 for u < 0 and given by a P × P matrix of functions. In this way, the intensity at an arbitrary time-point is dependent on the history of the multivariate process allowing for both self and mutually exciting behavior. In simple cases, if g i j (u) = 0 for all u and i = j, then this is non-cross exciting behavior. That is the cross-covariances are equal to zero and whilst the processes may be self exciting, they are independent and thus not mutually exciting (Hawkes 1971a). Further, in the bivariate case if g 21 (u) = 0 for all u but g 12 (u) = 0 for all u then N (1) (t) does not affect the likelihood of events in N (2) (t), but the converse does not hold.
In this case we have one way interaction between the two processes involved.
The CIF from (1) can be written for an exponential kernel as where ν, α = α pm , β = β pm , p, m = 1, . . . P are known as the baseline, excitation and decay parameters, respectively. Assuming stationarity of the multivariate Hawkes process we have that the vector of stationary densities is where G(ω) is the Fourier transform of the excitation kernel g(·), given by Note that G(0) is commonly referred to as the branching ratio, typically denoted by γ , with the condition for stationarity being that the spectral radius ρ(γ ) < 1 (Hawkes 1971b).

Multivariate MCEM for binned Hawkes processes
In a continuous time framework, maximum likelihood estimation (MLE) can be used to estimate the model parameters from a set of exact multivariate events on the interval [0, T ] (Ozaki 1979) denoted where T is the maximum observation or simulation time, t p l is the lth event in process p, and N ( p) (T ) is the total number of events in process p. When we observe an aggregation of these latent continuous times to a count process of events per time bin, we lose information and in particular the likelihood of our observed event times given a parameter set can no longer be computed exactly due to reliance on the underlying time-stamps. In particular, the Hawkes process which is defined by its CIF, depends on the history of the process, which if 'blurred' by binning or rounding requires correct handling in order to obtain meaningful parameter estimates. In the univariate setting, a method for parameter estimation based on the Whittle likelihood is presented in Cheysson and Lang (2022) and a modified Monte Carlo Expectation Maximisation (MC-EM) approach is presented in Shlomovich et al. (2022), termed the Binned Hawkes Expectation Maximisation (BH-EM) algorithm; which is extended here for multivariate data.

The Monte Carlo EM algorithm
The EM algorithm (Dempster et al. 1977) augments observed data by a latent quantity (Wei and Tanner 1990) to iteratively compute the maximizer of a likelihood. Here, the observed data are the multivariate event counts per unit The latent data, T are the unobserved, true event times. These are marked time-stamps which are not observed due to practical restrictions. We denote the parameter set to be Θ = {ν, α, β}, where ν is P × 1, and α and β are P × P. The following two steps are as detailed in Shlomovich et al. (2022): 1. In the E (Expectation) step, we compute where T denotes the sample space for T . 2. In the M (Maximization) step, maximize the conditional expectation in (2) to obtain the updated parameter estimate, Θ i+1 .
Monte Carlo methods can be used to numerically compute (2) if it is intractable, forming an algorithm known as MCEM (Wei and Tanner 1990). As we cannot sample the time-stamps directly we use importance sampling to simulate proposals T * for T from a feasible alternative distribution, denoted q(T | N, Θ i ).
Unique to the multivariate formulation is the need to retain the latent covariance structure between the P processes. To this end, we sample the times of the univariate superposition of the P-variate process. We then split the superposed simulated time-stamps to P processes to create proposals matching the multivariate counts. We refer to such viable proposals as consistent. This approach is detailed further in Sect. 2.2. Once sampled, each proposal is weighted according to 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 (2) with It is shown in Shlomovich et al. (2022) that if only proposing consistent event times where log p(T * (k) | Θ) is given in Daley and Vere-Jones (2003) by

Multivariate sampling via the superposition process
Generating consistent time-stamps for the latent multivariate process such that the true covariance structure is sufficiently captured is an important problem. We propose the following method for this, which allows us to extend the work in Shlomovich et al. (2022) for multivariate count data. Consider N t = [N 1,t , . . . , N P,t ] being a binned P-variate Hawkes process with exponential kernel. In order to generate possible multivariate, continuous-time proposal, T * , of the underlying event times, T , we first consider the superposed latent time-stamps, where T p represents the set of latent times for the pth process, for p = 1, . . . , P, and thus the superposed binned count process We map the multivariate parameter set Θ = {ν, α, β} to a corresponding setΘ = {ν,α,β} which approximates the univariate superposed process. It is given that the intensity of the superposed count pro-cessÑ can be decomposed as whereλ * (t) is the CIF for the superposed process. Using (5) and λ * p (t) defined in Eq.
(1) we have that We approximate the CIF of the superposed process with that of a univariate Hawkes process with exponential kernel and parametersΘ = {ν,α,β}. Doing so enables the cross and self excitation effects to be sufficiently captured by the internal correlation structure of the approximated superposed process, which is significantly simpler to generate proposals for. Specifically, whereλ * a (t) is the approximated CIF. By equating terms in Eqs. (6) and (7), it is clear that the baseline for the superposed process is simple to handle, and can be defined as ν = P p=1 ν p . The choice ofα andβ requires more care as, from the preceding two Equations, we approximate the contribution to the CIF from excitational effects as To aid the choice of appropriateα andβ we utilise the stationarity intensity of a univariate Hawkes process. We have that whereλ a is the stationary intensity forλ * a (t) andγ is the branching ratio forλ * a (t), equal toα/β due to the exponential form of the kernel. As the intensity of the superposed process can be decomposed, as in (5), the stationary intensities can similarly be decomposed. Therefore we have that The relationship in Eq. (9) aids the choice ofα andβ such that we retain a consistent stationary intensity under the approximation of the CIF with the exponential kernel given in Eq. (7). Further, the form in Eq. (9) is simple to compute as E{N p (T )} is well approximated by the mean number of events in each process. We find that a suitable choice forβ is where ω m is the proportion of total events which are observed in process m. This form is reached by considering Eq. (8).
The decay rate across events can approximately be given by the mean of the multivariate decays, weighted according to the proportion of events in each process. The benefit of weighting the mean can be seen by reasoning about the case when the processes have significantly different numbers of events. Combining (10) with the relationship forγ in Eq. (9) we havẽ It is important to note that as this reparameterisation does not rely on any latent continuous time-points, but rather just the observed counts and current multivariate parameter estimates, it is efficient to implement in practice. Empiricial studies show this reparameterization accurately recovers the true CIF of the superposed process. In this way we only simulate a univariate process, albeit a superposed version, and denote the simulated times asT * . In order to generate a realisation of the multivariate Hawkes process, T * , with cross-covariances, we can then uniformly sample the observed number of events in each bin for each process from theT * . In other words, we simulate a consistent set of continuous times for the superposed process matching the observed counts forÑ (t), and then uniformly assign points within each bin to each of the P processes. To reduce variance of the estimates, the allocation of events to the P processes can be conductedm times to generate multiple possible multivariate versions of the proposed realisation ofÑ . Ifm > 1, the sample which maximises the log-likelihood is selected, otherwise the single multivariate proposal is taken and the MLE is used to estimate the parameters of the multivariate process from the continuous-time proposed realisation. We present results form = 10 in Sect. 3.
This method provides us with an efficient way of artificially injecting cross-correlation into the consistent proposals whilst also retaining the marginal properties. In order to speed up the maximisation of the likelihood, we require the gradient and Hessian, given in "Appendix B". The full algorithm for this multivariate approach is derived and presented in "Appendix C".

Sampling method
The question remains of how to best sample the latent times. The sequential simulation method detailed in Shlomovich et al. (2022) is applicable in the multivariate extension due to the reparameterization step meaning we need only sample times for the univariate superposition process. In this case, the normalised weight of the kth Monte Carlo sample (k = 1, . . . , M) is given by where M is the number of Monte Carlo samples. Further we note that where log q T * (k) |Ñ ,Θ i is the log-likelihood of the sequentially sampled superposed times given the superposed counts and reparameterized univariate estimates, and Pr T * (k) |T * (k) , N denotes the probability of the random division of superposed time-points into a P-variate count process matching the observed counts.
Large differences between p T * (k) | N and q(T * (k) | N) can result in the weights being close to zero. Therefore we rescale the exponent term to avoid arithmetic underflow when computing the weights. We have Thus we can use Eq. (13) in place of Eq. (12).

Simulation study
We conduct a simulation study to compare the performance of the multivariate MC-EM algorithm to the INAR( p) method introduced in Kirchner (2016). We also compare the results to an approximation which ignores inter-bin excitation, referred to as the binned log-likelihood method. This approach represents the CIF as a piecewise constant function within each bin, equivalently assuming that N j are the counts in the jth bin of the pth process (Mark et al. 2019).
Given parameters ν, a 2 × 1 matrix, and α and β both being 2 × 2 matrices, along with some maximum simulation time T , we can simulate realizations of a Hawkes process. The generated events T represent the underlying process and aggregating these to a chosen binning Δ allows us to simulate the count data {N t , t = 1, . . . , K }. We then apply each of the multivariate MC-EM methods, INAR( p) and binned loglikelihood approximation.
Boxplots for each of the ten estimated parameters used for characterising the bivariate Hawkes process are given in Fig. 1. The parameters used for simulation are ν = 0.3 0.3 , α = 0.7 0.9 0.6 1.0 , β = 1.5 2.0 2.0 3.5 , with Δ = 1 and T = 2000. The parameters have been chosen as a stationary case with ample cross-excitation and nonsymmetric self-excitation. The mean parameter estimates from repeated simulations is presented on the vertical axis. The INAR( p) approximation method can yield highly variable results, which is to be expected as the method is primarily intended for selecting a parametric kernel from continuous time data. The usual implementation of the INAR( p) involves selecting Δ such that there is approximately one event per bin, however for our application this is not possible and so the choice of Δ here is chosen to better reflect real data. A log-scale has been used in all four graphs relating to each of α and β. Both the binned log-likelihood, and particularly the INAR( p) method produced large outliers, resulting in very large MSE relative to the MC-EM approach. Tables 1 and 2 present summary statistics for the outlier-trimmed data. Specifically we remove the top and bottom 5% of parame- ter estimates for each of the three methods and report the relative bias, mean and standard deviation for each of the parameters and methods. We see that even after removing outliers and negative values from the INAR( p) parameter estimates, the variability is much higher than that of the MC-EM. The absolute value of the bias of both the INAR( p) and binned log-likelihood approaches is also larger than that of the MC-EM, in some cases significantly so. Overall, the proposed MCEM approach has much improved estimation performance than the other methods. Table 3 additionally presents the ratio between the RMSE using T = 2000 and T = 1000 for each of the estimation approaches, using 300 Hawkes process realizations for each method. By fixing Δ, this allows the effect of doubling the number of bins K to be explored. As with the MLE, for the proposed MCEM method, it is expected that the RMSE reduces by approximately a factor of 1/ √ 2 when doubling K . We find that the mean RMSE ratio across the 10 parameters is 0.762 for the MCEM method, and 0.725 for the MLE acting on the continuous time points. However both the INAR( p) and binned log-likelihood methods have a significantly higher mean ratio, of 0.875 and 0.990 respectively.
Evaluation of the methods is additionally explored via goodness of fit. This is an important aspect which allows us to check the validity of the estimates given the data. Often, the random change theorem, given in Daley and Vere-Jones (2003), is used for considering goodness of fit by transforming time-points using the compensator function The theoretical expected relationship is ≈ 1/ √ 2 ≈ 0.707. The mean ratio across all 10 parameters is also given for each approach where t p k is the kth event in process p. In practice, λ * ( p) (u) is estimated using parameter estimatesΘ, and goodness of fit is conducted by considering the distribution of the transformed times, defined as T † = {t † 1 , t † 2 , . . .} = {Λ(t 1 ), Λ(t 2 ), . . .}.
By the random time change theorem, T † is a realization of a unit rate Poisson process if and only if T is a realization from the point process defined by Λ(·). Figure 2 shows QQ-plots relating to count processes from the simulation study. The continuous times from a single realization of a Hawkes process generated in the simulation study are transformed using each of the parameter estimates from the three approaches considered. The QQ-plots then compare the distribution of the interarrival times between each of the transformed times and a unit rate exponential distribution. We see that the estimates produced by the MC-EM algorithm are indeed a viable parameter set for this realization under an exponential Hawkes model, with the transformed time-points being distributed very close to the theoretical distribution. We also note that as the binned log-likelihood method ignores the effect of inter-bin excitation, it is expected that as the average number of counts in a bin increases, the fit from using this method will worsen. Whilst the QQ-plots given in Fig. 2 demonstrate the viability of the estimated parameters for a specific realization, in order to illustrate goodness of fit across all realizations in the simulation study, Fig. 3 also presents the distribution of Kolmogorov-Smirnov (KS) test statistics. Here, the KS test is used to test the equality of the distri- bution of the transformed interarrival times to an E x p(1) distribution for each of the three methods considered. This is done for all realizations used in the simulation study, and the distribution of the resulting KS test statistics is presented in Fig. 3. The MLE is also used to produce parameter estimates from the true, continuous time stamps, with the resulting KS test statistics also given in Fig. 3 for easier interpretation. It is clear that the MC-EM approach is far more comparable to the MLE results than the binned log-likelihood or INAR( p) methods. Note, variance estimation has not been discussed here, but parametric bootstrapping would allow for approximations of this to be obtained (Efron 1987).

Case study
Network flow data, referred to as NetFlow, assembles records exported by routers and describes communications between devices connected to an enterprise network. Monitoring and analysing NetFlow data has been successful at detecting a range of malicious network behavior (Turcotte et al. 2018).
Here we detect both self-exciting effects in the communication between a pair of network devices, termed here as an 'edge', and mutually-exciting activity between such edges in the Los Alamos National Lab (LANL) enterprise network. By modelling the activity in this way, insight into the correlation structure of communications in the network can be Counts obtained and monitored. We consider the methods outlined in Sect. 3 for parameter estimation of aggregated Hawkes processes as the NetFlow data is recorded at a 1 s resolution resulting in multiple events occurring simultaneously. Figure 4 presents a selected pair of edges in the LANL network with possibe mutually-exciting Hawkes behavior in the communications over a duration of 45 min, with a total of 171 events. We refer to the counts in the top plot as process 1 and the counts in the bottom plot as process 2. The window selected is chosen as a period of more frequent events surrounded by no activity.
Parameter estimates are found using the MC-EM, INAR( p) and binned log-likelihood methods and performance is considered using goodness of fit. This allows us to assess whether the estimates represent viable parameters for modelling the observed data as a mutually-exciting Hawkes process. To do this, we use the time rescaling theorem discussed in Sect. 3. Figure 5 shows QQ-plots for the observed times, uniformly redistributed within their respective bins so as to obtain continuous times for the purpose of assessing goodness of fit (Gerhard and Gerstner 2010). From Fig. 5 where denotes element-wise division and γ i j is the average number of events in process j directly triggered by each event in process i.
These results indicate that there is mutually-exciting behavior between these processes in one direction such that process 2 does not affect process 1 but process 1 does affect process 2, where process 1 is presented in blue in Fig. 4 and process 2 in red. There parameter estimates suggest that for each event in process 1, an average of 0.27 events will be triggered in process 2. The baseline parameters given by ν indicate the rate of the events is approximately one event The excitation kernel g(x) as estimated by the MC-EM algorithm. Note that only 3 curves are visible asĝ 12 (x) is practically zero every 100 s. Both processes are modelled with self-excitation, with each event in processes 1 and 2 triggering 0.34 and 0.28 events in their respective processes. Figure 6 shows the estimated excitation kernel components. From this we can see nature of the exciting effects, where g i j (x) illustrates the excitational effect which process j has on process i.

Conclusion
We have presented a novel method for generalising an MCEM algorithm to handle multivariate aggregated data. By reparameterizing our multivariate model in terms of the superposed process we can inject the necessary crosscovariance structure required for generating valid proposals. We also present closed form expressions for the gradient and Hessian of the log likelihood for increased computational efficiency. We further conducted a simulation study to compare this approach to the INAR( p) approximation detailed in Kirchner (2016) and a multivariate extension of the binned log likelihood method from Mark et al. (2019) and Shlomovich et al. (2022). As in the univariate case, the MCEM method out-performed both alternatives in the presented parameter set and moreover the bin width, Δ can vary provided the interval bounds are known. The multivariate extension can be applied for other Hawkes kernels.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
includes the optimal selection of p and exponential fit. The binned log-likelihood 0.46 s per realization.

Appendix B: Gradient and Hessian of multivariate continuous time Hawkes processes
The likelihood function of the multivariate Hawkes process is given by