MHP-DDP: Multivariate Hawkes Process with Dependent Dirichlet Process

Venue: NeurIPS_cc_2024_Workshop_BDU
Published: 2024-09-06 05:48:08

Abstract: Multivariate Hawkes Processes (MHPs) model complex temporal dynamics among event sequences on multiple dimensions. Typically, strong parametric assumptions are made about the excitation functions of MHP, motivating the need for modeling flexible excitation patterns. Further, different excitation functions across dimensions often have strong similarities. Motivated by reasons above, we propose MHP based on dependent Dirichlet process (MHP-DDP), a hierarchical nonparametric Bayesian modeling approach for MHP. MHP-DDP flexibly estimates the excitation function via a mixture of scaled Beta distributions, and borrows strengths across dimensions by modeling such mixing distribution as a mixture of a shared Dirichlet process (DP) and a group-specific idiosyncratic DP. We develop two algorithms using Markov chain Monte Carlo (MCMC) and the stochastic variational inference (SVI) algorithm. We also conduct simulations to compare MHP-DDP to benchmark methods where total or no information is borrowed. We show that MHP-DDP outperforms the benchmark methods in terms of lower estimation error for both algorithms, with SVI being  computationally efficient than MCMC.

Main paper: 
 introduction: 
 Point processes (e.g., see Cox and Isham, 1980) are a widely used to model temporal event sequence data, i.e., data corresponding to the times of occurrence of a series of countable events.In particular, multivariate Hawkes Processes (MHPs, e.g., see Hawkes, 1971 and Liniger, 2009) have been widely used to model sequences of events that demonstrate selfand mutually-exciting behaviors, i.e., patterns in which the likelihood of events increase after the occurrence of others.Hawkes Processes have been widely applied in a wide range of fields, including seismology (Ogata, 1988) , finance (Bacry et al., 2015) , electronic health records (Choi et al., 2015; Sun et al., 2024) and social media analysis (Rizoiu et al., 2017) . MHPs can be characterized through their conditional intensity functions, which describe the instantaneous rate of arrivals of new events. The excitation function is an important module of the conditional intensity function, as it controls how past events cause the conditional intensity function to change and decay. Excitation functions are often modeled parametrically using exponential functions that assume monotonic excitation decay (Hawkes, 1971) . However, the inter-event times in many real-world examples tend to have complex patterns, motivating the need for flexible excitation functions (e.g., see Markwick, 2020 and Rodríguez et al., 2017) .In this paper, we introduce a Bayesian nonparametric model for MHPs that builds on the ideas of Donnet et al. (2020) and Markwick (2020) but addresses various practical questions that have so far remained open. One challenge associated with the estimation of MHPs is that the number of parameters that need to be estimated grows quadratically with the number of dimensions. Hence, with a Submitted to Workshop on Bayesian Decision-making and Uncertainty, 38th Conference on Neural Information Processing Systems (BDU at NeurIPS 2024). Do not distribute. dataset of moderate size, modeling each dimension of the process independently can lead to inefficiencies. Motivated by the observation that often excitation functions looks similar across different dimensions, we propose a hierarchical modeling approach based on mixtures of non-parametric mixtures. More specifically, we adapt the approach introduced in Müller et al. (2004) , which models each of the excitation functions as a mixture of an idiosyncratic component and a common component shared by all excitation functions, and study some of the properties of such formulation.A second challenge in implementing MHP models is computation. As is the case more generally, Markov chain Monte Carlo (MCMC) algorithms are the most common approach to computation for Bayesian models for Hawkes processes (e.g., see Rasmussen, 2013) . However, MCMC algorithms for Hawkes process models are often to be too slow even for moderate sample sizes because, except for special cases such as the exponential excitation function, their complexity is quadratic in both the number of observations and the number of dimensions. This challenge is amplified in the case of nonparametric models. To address it, we expand on previous work on the use of stochastic gradient methods for MHPs, and develop a scalable stochastic variational inference (SVI) algorithm that can be used to fit our model. An important part of this development involves a carefully comparison of the performace of SVI and MCMC methods, both in terms of accuracy and speed.In summary we make two key contributions in this paper: (1) we propose a novel and flexible model for linear MHPs in which the various excitation functions are assigned a joint nomparametric prior that allows us to efficiently borrow information, (2) we develop MCMC and SVI algorithms for estimation and prediction in the context of this nonparametric model, and thoroughly evaluate their relative performance. 
multivariate hawkes processes: 
 Let N (1) (t), . . . , N (K) (t) be a collection of K point processes defined on the positive real line R + , where N (k) (t) represents the number of events on dimension k to occur on the interval [0, t]. We denote a generic set of observations from this process by X = {(t i , d i ) : i = 1, . . . , n}, where t i ∈ R + represents the timestamp at which the i-th event occurs and d i ∈ {1, . . . , K} represents the dimension in which the event occurs. Then, X follows a multivariate Hawkes process (Hawkes, 1971; Liniger, 2009) if the conditional intensity function on dimension k has the following form:λ k (t) ≡ lim h→0 E N (k) (t + h) -N (k) (t) | H t h = µ k + K k=1 ti<t,di=k α ℓ,k φℓ,k (t -t i ) ,where we µ k > 0 is the background intensity for dimension k, α ℓ,k > 0 be the parameter that controls the strength by which past events from dimension ℓ influence the occurrence of new events on dimension k, and φℓ,k (•) : R + → R + be the (normalized) excitation function that controls how such influence decays over time. Note that we require that ∞ 0 φℓ,k, (s)ds = 1, which ensures that φℓ,k are identifiable. An alternative construction of the MHP is as a multivariate branching process in which the first generation of events in dimension k (often called 'ìmmigrants" in the literature)arise from a homogeneous Poisson process with rate µ k , and the points in subsequent generations are generated from non-homogenous Poisson processes with rates given by the α ℓ,k s and the interarrival times are controlled by the ϕ ℓ,k s. In the sequel, we use the binary matrix B whereB j,j = 1 j-th event is an immigrant 0 otherwise , B i,j = 1 j-th event is generated from i-th event 0 otherwiseto encode the latent branching structure associated with a realization of an MHP. The augmented likelihood for the data is then given byL(X, B | α, ϕ) = K k=1 K ℓ=1     |O k,ℓ | (log α k,ℓ ) - i<j di=k,dj =ℓ B i,j φk,ℓ (t j -t i )     + K ℓ=1 |I ℓ | log µ ℓ - K ℓ=1 µ ℓ T - K k=1 K ℓ=1 α k,ℓ i:di=k Φk,ℓ (T -t i ), (1)where|I ℓ | = di=ℓ I(B ii = 1) and |O k,ℓ | = di=k,dj =ℓ,i<j I(B ij = 1) denote the number of immigrants for dimension k and the number of offpring on dimension k who arise from points on dimension ℓ.3 Nonparametric Bayesian modeling of excitation functions for multivariateFirst, we note that for all k, ℓ = 1, . . . , K, φk,ℓ can be interpreted as a probability distribution function on R + , with corresponding cumulative distribution function Φk,ℓ . Denote H as the space of all density functions on R, MHP-DDP defines a prior on φk,ℓ over H as a mixture of scaled Beta distributions, with respect to a random measure G k,ℓ (•):EQUATIONwhere the kernel density functionf Beta (• | a, b, T 0 ) := Γ(a + b) Γ(a)Γ(b) t T 0 a-1 1 - t T 0 b-1 T -1 0 , a > 0, b > 0, 0 < t < T 0is a scaled Beta density function with support on (0, T 0 ), indexed by shape parameters a and b. The Beta mixture allows the excitation function to have a flexible form, including multi-modal structures and heavy tails.We let G k,ℓ (•) be the mixing distribution for the two shape parameters that corresponds to φk,ℓ , which is a distribution on R 2 . Denote M(R 2 ) as the space of all probability distributions on R 2 , we define a prior on G k,ℓ over M(R 2 ) using the following hierarchical model, following Müller et al. (2004) :G k,ℓ (•) = εH 0 (•) + (1 -ε)H k,ℓ (•), k, ℓ = 1, . . . , K, H 0 , H k,ℓ ∼ DP(α DP , Gamma (c a , d a ) × Gamma (c b , d b )) k, ℓ = 1, . . . , K, ε ∼ Beta(1, 1).(3)Under the hierarchical prior, G k,ℓ can be expressed as a mixture of two random measures: the common component H 0 that is shared across all dimension pairs and the idiosyncratic component H k,ℓ that characterizes the dimension-pair specific behaviors of φk,ℓ . We let ε ∈ [0, 1] be the weighting factor that controls how much prior information on the characteristics of the distributions is borrowed across all dimension pairs, and assume it follows a uniform prior between 0 and 1. For implementation, we consider a finite number truncation approximation on the number of mixtures (Ishwaran and James, 2002) for the DP mixture considered in (2). Further, we assign Gamma priors to the Hawkes process parameters:µ ℓ | a ℓ , b ℓ i.i.d ∼ Gamma (a ℓ , b ℓ ) , α k,ℓ | e k,ℓ , f k,ℓ i.i.d ∼ Gamma (e k,ℓ , f k,ℓ ) , k, ℓ = 1, . . . , K.For the computational implementation of the model, we developed a Markov chain Monte Carlo (MCMC) method and a stochastic gradient variational inference (SVI) algorithm, which updates the variational parameters based on minibatches of the dataset.The data is generated from a multivariate Hawkes process with K = 2 dimensions and the following parameter settings: α = 0.6 0.15 0.3 0.6 , µ = 0.05 0.1 , T = 15000. For the triggering kernels, we consider a scenario where the Beta mixture component is generated from a mixture of two Beta kernels:φk,ℓ (t) = ε true exp(-t) + (1 -ε true ) exp(-λ j,k t),where λ 11 λ 12 λ 21 λ 22 = 2 0.8 0.8 2 . Additionally, ε true = {0, 0.5, 1} represents the true degree of information borrowing across dimensions. We fit the model using two methods: a Metropolis-within-Gibbs exact full-batch sampler based on Markov chain Monte Carlo (MCMC) and a stochastic gradient variational inference (SVI) algorithm based on minibatches. We also compare our method with random ε ('RANDOM'), to two benchmark versions of both MCMC and SVI where there is no information borrowing ('IDIO') or the triggering kernels are identical ('COMMON'). Additionally, we also consider a frequentist benchmark method based on piecewise basis kernels using the EM algorithm (EM-BK, see Zhou et al., 2013) . We evaluate both the point and uncertainty estimation accuracy for the triggering kernels, using a set of performance metrics. We use the root mean integrated sqaured error (RMISE) as a metric for point estimation accuracy:RMISE(ϕ) = 1 K 2 K k=1 K ℓ=1 +∞ 0 ϕ true k,ℓ (x) -φk,ℓ (x) 2 dxFigure 1 shows the results and their specific scenarios. It can be shown that our methods have the lowest estimation error for both methods under most scenarios (except for MCMC when ε = 1, where the difference is very small). For the computation costs, SVI took only less than hour to converge while the MCMC algorithm takes over a day, further suggesting that the usage of scalable 1 : RMISE as point estimation metric for all methods under true information-borrowing ratios. The values in the grid cells are averaged over 10 independently datasets, and the standard deviation is shown in the brackets. The numbers in bold refers to the best-performing method among 'RAN-DOM', 'IDIO' and 'COMMON' for both MCMC and SVI methods. 
conclusion: 
 In this work, we developed a novel Multivariate Hawkes processes model for complex temporal event data. Especially, we flexibly modeled the decaying patterns of the triggering kernels using a Dependent Dirichlet process mixture of Beta distributions. Our model shows favorable results in both point and uncertainty prediction methods compared to benchmark models, and could serve as a basis for forecasting and decision making. We developed MCMC and SVI methods for computation, with SVI being computationally efficient and scalable to large datasets. 
 
 Appendix: We first prove the Kullback-Leibler property for a generic Dirichlet process mixture of Beta kernels model in Theorem 1, and then extend to our model setting in Corollary A.1.Theorem 1. Let H 0 be the set of all continuous densities on [0, 1]. If f 0 (x) ∈ H 0 , Π is the prior induced by the Beta likelihood mixture kernel and DP(α DP , Gamma (c a , d a ) × Gamma (c b , d b )) on H 0 , Then f 0 ∈ KL(Π), i.e. there exists ε > 0 such that Π({g : KL(f 0 , g)}) < ε.Proof. We show the KL property of f 0 by proving Conditions A1-A3 from Theorem 1 in Wu and Ghosal (2008) holds. Since we don't have ϕ in our case, Condition A2 is automatically satisfied.The remaining proof is similar to Theorem 11 in (Wu and Ghosal, 2008) , except that the mixing distribution is drawn from a Dirichlet process, and that the two location parameters are also being mixed. For ∀f 0 ∈ H 0 and ∀ε > 0, there exists a finite Beta mixture function f Pε with H mixtures whereEQUATIONandEQUATIONsuch that1 0 f 0 (x) log f 0 (x) f Pε (x) dx < ε.Thus Condition A1 holds. We then show Condition A3 also holds. First, we define C ω ⊂ S H-1 , C h ab ⊂ R 2 , h = 1, . . . , H as sets such thatEQUATIONwherex M = max d, 1 -d, a 0 h -a h a 0 h -a h + b 0 h -b h . Finally, we let C ab = ⊕ H h=1 C h ab , ω := (ω 1 , . . . , ω H ) ∈ S H-1 , a := (a 1 , . . . , a H ) ∈ R H , b := (b 1 , . . . , b H ) ∈ R Hand let W ⊂ H be the set of finite mixture distributions induced by C ω and C ab :W := P ∈ W | P = H k=1 ω k δ a,b (a k , b k ) , ω ∈ C ω , (a h , b h ) ∈ C h ab .We note that for the chosen ω 0 , a 0 , b 0 , for any ω ∈ C ω , we haveEQUATIONWe then note thatEQUATIONThus for any (a h , b h ) ∈ C h ab , consider the Cauchy remainder form of the first-order expansion forl(x, a 0 h , b 0 h ) := log f Beta (x; a 0 h , b 0 h ) at (a h , b h ), we have l(x; a 0 h , b 0 h ) = l(x; a h , b h ) + ∇l(x; a h , b h ) T a 0 h -a h b 0 h -b h + θ 2 2 a 0 h -a h b 0 h -b h ∇ 2 l(x; a * h , b * h ) a 0 h -a h b 0 h -b h < l(x; a h , b h ) + (a 0 h -a h ) (log x + ψ(a h + b h ) -ψ(a h )) + (b 0 h -b h ) (log(1 -x) + ψ(b h + a h ) -ψ(b h )) < l(x; a h , b h ) + (a 0 h -a h ) log x M + (b 0 h -b h ) log(1 -x M ) < l(x; a h , b h ) + ε 4 , ∀x ∈ (d, 1 -d).(9) Note that the chain of inequalities are based on equations ( 6) and ( 8), and the fact that ∇ 2 l(x; a, b) ≺ 0, a > 0, b > 0. Finally, note that ( 9) is equivalent toEQUATIONand if for h = 1, . . . , H, (a h , b h ) ∈ C h ab , for any (ω 1 , . . . , ω H ) ∈ C ω , we haveH j=1 ω h f Beta x | a 0 h , b 0 h H j=1 ω h f Beta (x | a h , b h ) < e ε 4 , ∀x ∈ (d, 1 -d).Combining equations ( 7) and ( 10), we have, for any ωEQUATION)We then consider the scenario where x ∈ (0, d] [1 -d, 1). We want to show that the likelihood f Beta x | a 0 h , b 0 h f Beta (x | a h , b h ) = Be(a h , b h ) Be(a 0 h , b 0 h ) sup x∈(0,d] [1-d,1) x a 0 h -a h (1 -x) b 0 h -b h = Be(a h , b h ) Be(a 0 h , b 0 h ) d a 0 h +b 0 h -a h -b h ≤ Be(a h , b h ) Be(a 0 h , b 0 h ) , which follows from the fact that a h < a 0 h , b h < b 0 h . Thus we have sup a,b∈C ab ω∈Cω sup x∈(0,d] [1-d,1) H j=1 ω 0 h f Beta x | a 0 h , b 0 h H j=1 ω h f Beta (x | a h , b h ) ≤ e -ε 4 sup a,b∈C ab Be (a h , b h ) Be (a 0 h , b 0 h ) ≤ sup a,b∈C ab Be (a h , b h ) Be (a 0 h , b 0 h ) := M < +∞.Thus we haved 0 f 0 (x) log f Pε (x) f P (x) dx + 1 1-d f 0 (x) log f Pε (x) f P (x) dx < M (F 0 (d) + 1 -F 0 (1 -d)) .Thus, we can choose d small enough such that M (F 0 (d) + 1 -F 0 (1 -d)) < ε 2 , such that we have, Π * {f 1 , . . . , f K : KL(f k 0 , f k ) < ε, k = 1, . . . , K} > 0.Proof. Let f 0 0 (x) := 1 K K k=1 f k 0 (x). For ε > 0, we letG 0 = {g 0 ∈ H 0 : KL(f 0 0 , g 0 ) < ε 2 }, G k = {g k ∈ H 0 : KL(f k 0 , g k ) < ε 2 }.Note that from Theorem 1.1 we have Π g k g k ∈ G k > 0, k = 0, 1, . . . , K. Let M := max 1≤j,k≤M sup g 0 ∈G0 KL(f 0 ij , g 0 ) < +∞, we have, based on the convexity of KL divergence:KL(f k 0 , δg 0 + (1 -δ)g k ) ≤ δKL(f k 0 , g 0 ) + (1 -δ)KL(f k 0 , g k ) < δM + (1 -δ) ε 2 < ε,for all ∀δ < ε 2M -ε , g k ∈ G k , k = 0, 1, . . . , K. Thus Π * f 1 , . . . , f K : KL f k 0 , f k < ε, k = 1, . . . , K =Π * g 0 , g 1 , . . . , g K , δ : KL f k 0 , δg 0 + (1 -δ)g k < ε, k = 1, . . . , K >Π * g 0 , g 1 , . . . , g K , δ :g k ∈ G k , k = 0, 1, . . . , K, δ < ε 2M -ε = 0≤k≤K Π g k (g k ∈ G k )Π δ (δ < ε 2M -ε ) > 0. 



Reviews: review: title: title: Blind Review
main_review: main_review: This paper proposes two algorithms to improve the efficiency of the MPH algorithm. The authors experimentally verify their effectiveness compared to the state of the art.

I am inclined to accept this paper. It covers a relevant topic area, the design seems novel compared to recent work, and the presentation is mostly clear. Performance is tested on previously used synthetic datasets, and the model is well-described.

A potential weakness of the paper is that the evaluation is limited to a single dataset. It would have been better to include higher-dimensional datasets, even if they were synthetic or from meme tracker datasets.

Other notes:

Line 120: There is a typo in “RMISE,” which stands for "Root Mean Integrated Squared Error (RMISE)."
title: title: It's a short paper
main_review: main_review: The paper proposes a novel approach called MHP-DDP for modeling complex temporal dynamics among event sequences on multiple dimensions using Multivariate Hawkes Processes (MHPs). The proposed approach is based on a dependent Dirichlet process (DDP) and allows for flexible estimation of the excitation function by using a mixture of scaled Beta distributions. The paper presents two algorithms for implementing the MHP-DDP model and compares its performance to benchmark methods through simulations. The results show that MHP-DDP outperforms the benchmark methods in terms of lower estimation error, with the stochastic variational inference (SVI) algorithm being more computationally efficient than the Markov chain Monte Carlo (MCMC) algorithm.
Pros:
This paper presents a novel and adaptable model for Multivariate Hawkes Processes (MHPs) that enhances flexibility by allowing the sharing of information across dimensions and by estimating the excitation function using a mixture of scaled Beta distributions. This approach offers greater adaptability compared to traditional parametric models.
The study introduces two algorithms, Markov Chain Monte Carlo (MCMC) and Stochastic Variational Inference (SVI), for the implementation of the MHP-DDP model. Notably, the SVI algorithm proves to be computationally efficient and capable of scaling to accommodate large datasets.
The paper conducts a comprehensive evaluation of the proposed model, comparing it against established benchmark methods through detailed simulations. The findings highlight the MHP-DDP model's superior performance, particularly in achieving lower estimation errors.
Cons:
Please provide a deeper analysis of the limitations and potential drawbacks of the MHP-DDP model. What assumptions or conditions might impair the model’s performance?
Could you elaborate on the theoretical properties of the MHP-DDP model? Specifically, how does it compare to or differ from existing models cited in the literature?
To validate the MHP-DDP model's performance, please consider utilizing real-world applications, case studies, or datasets. Additionally, could you outline any plans for future research in this area?
I request more comprehensive details regarding the simulations used to compare the MHP-DDP model with benchmark methods. This should include the number of data points analyzed, the specific metrics applied for evaluation, and any underlying assumptions. Could you clarify which metrics were employed and explain their calculation methods?
Please provide additional information on the motivation behind the proposed approach and delineate how it diverges from existing methodologies.
Can you discuss the scalability of the proposed SVI algorithm with respect to data size and dimensional complexity? Have there been any tests of the algorithm on larger datasets?
Please offer further insights into the theoretical basis and advantages of employing a dependent Dirichlet process and a mixture of scaled Beta distributions within the MHP-DDP model.