Learning Neural Jump Stochastic Differential Equations with Latent Graph for Multivariate Temporal Point Processes

1 Introduction

Event data typically consists of extensive asynchronous event sequences, playing a crucial role across various domains from scientific research to industrial applications, underscoring its undeniable importance. Each event within these sequences is characterized by a timestamp and a type mark, indicating both when and what the event occurred. For instance, in seismology, advancements in technology have enabled the collection of large-scale earthquake event data, aiding scientists in better understanding and predicting seismic activities [Reinhart, 2018]. In financial, financial transaction records serve as event data widely utilized for detecting anomalies and forecasting market trends [Du et al., 2016]. On social media platforms, capturing social interaction data provides valuable insights, assisting in understanding user behavior and conducting sentiment analysis [Tung et al., 2016]. Temporal Point Processes (TPPs) represent a crucial tool for studying event sequences, primarily focusing on the occurrence times of events and the time intervals between events. Furthermore, for event sequences with multiple types of events, multivariate TPPs (MTPPs) can be applied to reveal the correlations, patterns, and trends among different event types.

∗Corresponding author.

The traditional TPPs typically rely on certain assumptions when defining the intensity function. For instance, the Hawkes process defines a self-exciting process indicating that

rences [Hawkes, 1971]. The self-correcting process assumes that the occurrence of past events inhibits the occurrence of future events [Isham and Westcott, 1979]. In a Poisson point process, the probability of an event occurring at each time point is constant and independent of past events [Daley and Vere-Jones, 2003]. These traditional TPP models provide a foundational framework for studying event sequences, but they also share some common limitations. These limitations include simplifying assumptions about the correlations between events, oversimplification of the dynamic evolution of event sequences, and ignoring complex patterns and nonlinear relationships. Consequently, in recent years, methods based on the powerful representational capabilities of neural networks have gradually emerged.

Due to advancements in deep learning, Neural Point Processes (NPPs) have gained significant attention for employing neural networks to model the intricate dynamics of asynchronous event sequences. For instance, Omi et al. introduce a fully neural network-based method to model the integral of the intensity function and obtain the intensity function as its derivative [Omi et al., 2019]. Shchur et al. propose an intensity-free method that learns without explicitly modeling the intensity function, reducing reliance on specific forms [Shchur et al., 2019]. Then, transformerbased approaches are proposed to use attention mechanisms to model MTPPs [Zuo et al., 2020]. However, these approaches still exhibit limitations when modeling continuoustime processes, particularly in terms of ensuring continuity and smoothness. Consequently, techniques rooted in Ordinary Differential Equations (ODEs) and Stochastic Differential Equations (SDEs) have showcased competitive performance [Song et al., 2024; Zhang et al., 2024]. Nevertheless, these methods do not explore the impact of underlying relationships between events, whereas graphs play a crucial role in predicting the spatiotemporal evolution of complex systems [Pei et al., 2022]. Therefore, the main motivation of this paper is to enhance the modeling capability of MTPPs by incorporating a latent graph over different types of events.

This paper presents a novel jump-diffusion SDE method with a latent graph for modeling MTPPs. Firstly, inspired by the equivalent SDE formulation of the multivariate Hawkes process, multi-dimensional SDEs are proposed to represent the hidden state of the intensity function of each event type. Then, we integrate an encoder-free graph representation method and design interaction functions between events in the jump process to capture potential relationships between events, which can further enhance the representation and predictive abilities of event sequences. Moreover, the proof of existence and uniqueness of solutions enhances the reliability of our approach.

Our main contributions are as follows:

• Multi-dimensional SDEs are employed to model the intensity function for each event type, enhancing the representation capability of the underlying MTPPs.
• To further capture the underlying relationships between events, a graph generator without encoders is incorporated. Subsequently, interaction functions between event types in the jump process are also proposed.
• We prove the existence and uniqueness of the solutions of our approach. Meanwhile, experimental results show that the performance of our approach is better than other SOTA methods.

2 RelatedWorks

2.1 Neural Point Processes

Due to certain limitations in traditional TPPs, such as the definition of intensity functions or kernel functions, NPPs have gradually become a focus of research. Early methods are mainly based on Recurrent Neural Networks (RNNs) to model point processes by assuming a specific functional form for the time evolution of the intensity function of a point process [Du et al., 2016; Mei and Eisner, 2017]. To overcome this limitation, Omi et al. propose a method that directly models the cumulative hazard by neural networks and obtains the intensity function with the automatic differentiation techniques [Omi et al., 2019]. Subsequently, several methods have been proposed that leverage Transformer architecture networks to capture the dependencies between events [Zuo et al., 2020; Zhang et al., 2020; Yang et al., 2022; Meng et al., 2024]. With the development of neural ODEs and SDEs, methods based on modeling continuous systems rather than other architectures like RNNs can explicitly model the smooth transitions between states of intensity functions [Chen et al., 2018]. Jia et al. propose NJSDE, a neural ODEs-based hybrid system that encompasses both flow and jump processes, to model TPPs [Jia and Benson, 2019]. Song et al. introduce Dec-ODE to decouple the hidden state dynamics by neural ODEs [Song et al., 2024]. Zhang et al. present NJDTPP, where the intensity process is designed by a neural jump-diffusion SDE [Zhang et al., 2024]. It is worth noting that both NJSDE and NJDTPP also design the impact of event occurrences in their jump items. However, neither of them continues to deeply consider the underlying relationships between event types. Consequently, this oversight may lead to the magnification of the impact between unrelated events. Furthermore, NJSDE uses ODEs to model the overall system rather than each event type, while NJDTPP employs one-dimensional SDEs for each event type. In contrast, our innovative approach utilizes multidimensional SDEs to capture more complex temporal dynamics for each event type. These improvements enhance the predictive and modeling capabilities of our approach.

2.2 Temporal Point Processes with Graph

In traditional TPPs, besides establishing a probabilistic graphical model based on existing knowledge, some methods aim to infer relationships among events. For instance, Zhou et al. propose a method to infer latent connectives between users in social networks by a multivariate Hawkes process [Zhou et al., 2013]. Lemonnier et al. propose a Low-Rank Hawkes Process framework for fitting multivariate Hawkes processes in large-scale problems [Lemonnier et al., 2017]. With the emergence of NPPs, some NPPs models combining graphs have appeared. Zhang et al. adopt a variational inference framework to model the posterior relation of MTPP data for probabilistic estimation [Zhang and Yan, 2021]. Zhang et al. propose a probabilistic graph generator. The sampled graph can be easily integrated as a plug-in to modify an existing NPP model [Zhang et al., 2021]. Yoon et al. extract potential graph structures through data processing and encode multivariate event sequences into a sequence of graphs, enhancing the predictive capability of the model [Yoon et al., 2023]. Yang et al. partition event sequences into different sub-intervals based on prior knowledge and establish dynamic graph relationships using a Variational Autoencoder (VAE) to capture changes in relationships between different event types [Yang and Zha, 2024]. In a nutshell, NPPs with latent graphs have effectively enhanced the interpretability and predictive capability of the model. However, SDE-based methods have not yet explored this aspect, which is also the motivation of our method.

3 Preliminaries

3.1 Multivariate Point Processes

A MTPP is a stochastic process that describes the temporal evolution of multiple event types V = {1, 2, ..., V }, where V is the number of event types. Let S = {(vi, ti)}Li=1 be an event sequence, where the tuple (vi, ti) is the i-th event, vi ∈ V and ti ∈ [0, T] are its event type and timestamp, T is the size of the time window. A general MTPP can be characterized by the conditional intensity function of each event type k ∈ V:

where dNk(t) = Nk(t +∆t) -Nk(t) ∈ {0, 1} is the jump size, Nk(t) counts the number of the event type k prior to history H(t) until time t, H(t) = {(v′ , t′ )|t′ < t, v′ ∈ V} records the history. The log-likelihood of observing a sequence S is the sum of the log-likelihood of events and nonevents and can be expressed as Eq. (2). The optimization of MTPPs is to maximize the log-likelihood.

3.2 Multivariate Hawkes Processes

We then introduce a classic and important MTPP known as the Multivariate Hawkes Process (MHP). Let N(t) = (N1(t), . . . , NV (t))T be a MTPP, where its intensity function λk(t) for each event type k is defined as follows:

where µk ∈ R+ represents the base intensity, and g(t) ∈ R+ denotes a kernel function that signifies the extent of influence of an event. The entry αjk indicates how event type j effects event type k.

3.3 Stochastic Differential Equations

A basic SDE can be defined as:

where x0 is the initial value, f is the drift coefficient function, g is the diffusion coefficient function, {W(t)}t≥0 is a Brownian motion. Then, X={X(t)}t≥0 is called a solution of the SDE, and it satisfies the Itˆo integral equation.

In the context of neural SDEs, the selection of the diffusion function is especially crucial. To make training stable, this paper also considers a Linear Noise SDE (LNSDE) [Oh et al., 2024], presented as follows.

It has been proved that, with appropriate neural network designs, the LNSDEs in Eq. (6) have their unique strong solutions. Furthermore, neural LNSDEs also demonstrate robustness in maintaining excellent performance and effectively preventing overfitting.

3.4 Equivalent SDE Formulation for MHPs

To establish a connection between SDEs and MTPPs, the following derived equivalent SDE formulation is introduced.

Theorem 1. Suppose that the kernel function g(t) = exp(-β(t)), the MHPs N(t) and its intensity functions

λ(t) = (λ1(t), . . . , λV (t))T can be equivalently expressed as the solution to the jump SDEs. Each λk(t) follows

Proof. The proof is similar to the univariate Hawkes Process discussed in [Zhang et al., 2024], which can be found in Appendix A.1.

4 Methodologies

4.1 Overview

As shown in Fig. 1, the proposed method involves SDEs with jumps designed for MTPPs with a latent graph. To better capture the internal states within each event, every event type is associated with multidimensional hidden states that evolve according to an SDE system. Moreover, we employ a graph generator without any encoder and integrate the influence of the latent graph into the jump process.

4.2 Proposed SDEs for Intensity function

Due to the diversity of underlying dynamical systems, we introduce two different variants of SDEs for the hidden state in our model as Eqs. (8) and (9), which denote whether the diffusion term is time-independent or not. In our model design, we assume that each type of event k ∈ V has its corresponding hidden state Zk, and Zk(t) evolves as a SDE. To more accurately capture the temporal effects, the continuous SDE between events [ti-1, ti) is designed as follows:

where f : RdZ × R+ × R+ → RdZ considers not only the hidden state Zk(t) but also the time ti-1 and the time since the last event occurred t - ti-1, while the diffusion term g : RdZ → RdZ is dependent on the hidden state Zk(t), {Wk(t)}t≥0 is a dZ-dimensional Brownian motion. Also, the version of neural LNSDE can be expressed as follows.

where the function g : R+ → RdZ is dependent on the time t-ti-1. For the convenience of discussion, we use Eq. (8) as the default. Furthermore, the proposed approach can directly adopt Eq. (6).

Then, the intensity function λk(t) of each event type k is defined as follows.

where η : RdZ → R, softplus(x) = log(1 + exp(ζx))/ζ ensures the λk(t) > 0 and ζ is a hyper-parameter.

4.3 Latent Graph in Jump Process

From Theorem 1, it can be intuitively observed that the influence between different types of events takes effect when events occur. Therefore, compared to incorporating the effects of node interactions into the drift coefficient function, our proposed approach emphasizes the latent relationships between different types of events in the jump process. Moreover, this method can reduce computational costs by passing messages between nodes solely at the times of event occurrences, rather than continuously throughout all time.

Firstly, let us consider the situation without any graph structure. In this case, the occurrence of an event of one type will have a certain impact on all other types of events. Therefore, when an event type j occurs at time t, then the jump size ∆Zk(t) = Zk(t+)-Zk(t) of event type k can be formulated as Eq. (11).

where t+ is the right-limit of t, the function ρ : RdZ ×RdZ → RdZ is used to capture the interaction between different types of events. Furthermore, in general, the form of the jump size can be extended as

However, our proposed approach incorporates a latent graph in the jump process to emphasize capturing potential relationships between event types. Trying to capture events that are practically unrelated may result in a decrease in the model's representation and predictive capabilities. Therefore, we adopt a method to integrate a latent graph without an encoder inspired by [Pan et al., 2024]. Additionally, we design interaction functions to enhance the model's ability to capture relationships between different event types.

The latent graph 1 is constructed as follows.

where a ∈ {0, 1}, Ψ0 and Ψ1 ∈ R are trainable parameters, and γ is the inverse temperature, and A1 = 11T - A0 . Instead of using the initial states Z(0) to build the latent graph, this way can decouple the graph structure from the hidden states to enhance training stability, especially when dealing with long sequences.

Then, the jump size of event k is designed by incorporating the latent graph as shown in Eq. (14).

where ˆ Zjk(t) is the interaction term between event type vj and vk under the latent graph, ρv : RdZ → RdZ is the jumpwise function, ρn : RdZ × RdZ → RdZ is the aggregationwise function and ρe(a) : RdZ ×RdZ → RdZ is the edge-wise function for adjacency matrix Aa .

4.4 Theoretical Analysis

Based on the two sections mentioned above, our proposed approach can be expressed as follows.

where Z(t)=(Z1(t), ...,ZV (t))T is a collection of latent state vectors, W(t)=(W1(t), ...,WV (t))T is a collection of Brownian motions, its components Wv(t) are pairwise independent, ∆Z(t)=(∆Z1(t), ...,∆ZV (t))T are calculated by Eq. (14), t and ti-1 are also column vectors. We then proceed to investigate the existence and uniqueness of the solution {Z(t)}t≥0 and {λ(t)}t≥0 to our proposed approach.

Theorem 2. Assuming that f , g, ρv , ρn , ρe , η are measurable functions and satisfy the Lipschitz continuity condition, ρv , ρn , ρe are continuous, then a strong unique solution Z := {Z(t)}t≥0 exists, once initial values z := {zk}Vk=1 is fixed, and the resulting process λ := {λ(t)}t≥0 also has a strong unique solution.

1 In this context, the term "graph relationships" does not refer to the ground-truth causal relationships but rather to Granger causality, which aims to assess the utility of one variable in predicting other variables.

4.5 Model Training

Throughout the training process, we utilize the EulerMaruyama scheme [Kloeden et al., 1992] with a constant step size for its computational efficiency. This scheme subdivides the time interval (ti-1, ti] into D sub-intervals, where ti-1=τ i0<· · · <τ id<· · · <τ iD=ti with the fixed step size ∆id = (ti - ti-1)/D. Then, for each hidden state Zk, Eq. (8) on (ti-1, ti] can be discredited by the recursive equation as Eq. (16).

for d ∈ {0, 1, ..., D - 1}, and Zk(τ i0 ) = Zk(t+i-1 ), ∆W i d = W(τ id+1 ) -W(τ id ) is sampled fromN(0,∆id ) for numerical computation. Finally, the log-likelihood can be approximated as follows.

where τ 1 0 =0, τ L+1 D =T, λk(t) can be calculated by Eqs. (10) and (15).

4.6 Time and Event Prediction

To predict the subsequent event time and event type based on the historical events H(ti), the conditional density function of the next event time for ˆ ti+1 ≥ ti is estimated as Eq. (18).

where λ(t) = ∑Vk=1 λk(t).

For the next event type prediction ˆ vi+1, it can be estimated by Eq. (19).ˆ vi+1 = argmaxk∈V λk(ti+1)/λ(ti+1). (19)

5 Experiments

5.1 Datasets

We evaluate our proposed approaches and other baselines on five real-world benchmark datasets used in [Yang and Zha, 2024; Zhang et al., 2024], including New York Motor Vehicle Collisions (NYMVC), two stack exchange datasets (MathOF and AskUbuntu), Taobao and Taxi. The statistical characteristics of these datasets are shown in Table 1. Additionally, to scale down the time, a logarithmic transformation is applied to all time values. We also ensure that at most one event occurs at any given time to facilitate the training and prediction processes for all models. More detailed information of datasets can be found in Appendix B.

5.2 Baselines

We compare our approaches with the following eight state-ofthe-art models. One RNN-based model: Recurrent Marked Temporal Point Process (RMTPP) [Du et al., 2016]. One TPP with fully neural networks (FullyNN) [Omi et al., 2019]. Two attention-based models: Transformer Hawkes Process (THP) [Zuo et al., 2020], Attentive Neural Hawkes Process (AttNHP) [Yang et al., 2022] and Interpretable Transformer Hawkes Processes (ITHP) [Meng et al., 2024]. One log-normal mixture distribution-based model: Variational Autoencoder for MTPPs with Dynamic Latent Graphs (VAETPP) [Yang and Zha, 2024]. One ODEs-based method: (NJSDE) [Jia and Benson, 2019]. One SDEs-based method: Neural Jump-Diffusion Temporal Point Processes (NJDTPP) [Zhang et al., 2024]. The detailed description and comparison of baselines can be found in Appendix C.

5.3 Experimental Settings

For our approaches, both f and g in our approach are implemented as three-layer Multilayer Perceptions (MLPs), which contain an input layer, a hidden layer, and an output layer. In the jump process, ρv , ρn , and ρe are two-layer MLPs, each consisting of an input layer and an output layer. The activation function default used is Tanh. The dimension dZ of Z is set to 32, and the dimension of hidden layers is uniformly set to 64. The learning rate of parameters in neural networks is set to 0.001 for the NYMVC dataset and 0.0001 for other datasets, and the learning rate of initial hidden states z and Ψ is set to 0.1. We employ the Adam optimizer [Kingma and Ba, 2015] to optimize the parameters. During the training process, the number of sub-intervals D in the Euler-Maruyama scheme is set to 10, while during the prediction process, D is set to 1000 for a more precise forecast. More details can be found in Appendix D. As for the other comparative baselines, their settings are adopted from the original papers.

5.4 EvaluationMetrics

We utilize the negative log-likelihood (NLL) as a metric to assess the modeling capability of event sequences in real-world datasets. Furthermore, we assess performance in the standard next-event prediction task within TPPs, predicting each subsequent event (vi, ti) based on the history H(ti). Event time prediction is evaluated using Root Mean Square Error (RMSE), while event type prediction is assessed through accuracy and the weighted F1 score. The accuracy in multiclass classification represents the proportion of correct classifications. The weighted F1 score provides a better assessment of class imbalance situations. These metrics are shown in Eqs. (20), (21) and (22).

where C denotes the number of classes, wi denotes the sample weight for class i, and F1i represents the F1 score for class i.

5.5 Experimental Results

In Table 2, the performance of the baselines and our proposed methods is presented in terms of negative log-likelihood and time prediction RMSE on various datasets. Concerning NLL, our proposed approaches consistently achieve some of the highest rankings across all datasets, showcasing the robust modeling capabilities of our approach for MTPPs. NJDTPP and VAETPP also demonstrate commendable performance. In the realm of time prediction, our proposed methods and VAETPP exhibit strong competitiveness. VAETPP utilizes a mixture of log-normal distributions to model the distribution of one-dimensional inter-event times, enhancing its capability for time prediction. Consequently, its time prediction accuracy leads on some datasets, while our approaches also perform well on the majority of datasets.

Table 3 illustrates the accuracy of predicting event types using various methods on real-world datasets. The results indicate that our methods, SDE and LNSDE variants, outperform other baselines in terms of prediction accuracy across all datasets, underscoring the effectiveness of our methods tailored for MTPPs. One reason for the lower accuracy of event prediction in NJSDE is that it models the entire system based on ODEs, making it difficult to accurately represent the hidden states of each event. Furthermore, SDE and LNSDE variants demonstrate varying performances across different datasets, highlighting the importance of careful consideration when selecting diffusion terms. Suitable diffusion term selections will further enhance prediction accuracy. Additionally, the competitive performance of VAETPP and NJDTPP is also noteworthy. The discernible gap between accuracy and weighted F1 arises from the imbalance in event categories, particularly evident in the AskUbuntu dataset.

Overall, our proposed methods have demonstrated strong capabilities in modeling, time prediction, and event-type prediction on various real-world datasets, further illustrating the effectiveness of our proposed methods.

5.6 Ablation Study

Our ablation experiments are designed as follows: ODE-NA signifies that the version of our approach without both the Brownian motion and the latent graph in the jump process, as depicted in Eq. (11). The ρ function is also implemented through three-layer MLPs. ODE denotes the version of our approach without Brownian motion.

The results shown in Table 4 show that the performance is significantly improved by incorporating a latent graph in the jump process. This also demonstrates that considering the dependency levels among different event types can enhance the modeling and predictive capabilities of the model. Furthermore, utilizing SDE will further enhance the accuracy of modeling and prediction.

5.7 Hyperparameter Analysis

Number of Sub-intervals

For the solver of SDEs (e.g., Euler-Maruyama method in our approach), a large number of sub-intervals results in a more accurate approximation of the integral but also linearly increases the time complexity. Following our preliminary experiments, we generally found that there is a good balance between training time and effectiveness when the number of sub-intervals D=10. For instance, the impacts on our LNSDE variant on NYMVC are in Table 5.

Dimensions of Hidden State Zk(t)

We also conduct experiments regarding the dimensions of the hidden state Zk(t). For instance, the performances of our LNSDE variant on NYMVC and Mathof are in Table 6.

The results show that increasing the dimensionality provides significant improvement in predicting events. However, overly large dimensions can harm time prediction, since the model may memorize noise or over-fit instead of capturing true temporal patterns.

5.8 Model Interpretability

The learned latent dependencies among events can offer some interpretability to the model. As a result, we conduct an experiment on a synthetic network to present the interpretability of our model.

We initially create a directed graph with 20 nodes. Then, we simulate a multivariate Hawkes process 100 times. Here, the parameters µ and α are randomly sampled from uniform distributions U[0, 0.02] and U[0.1, 0.3], respectively, while the time delays follow an exponential distribution with a parameter of β = 2.5. After training, we acquire the edge probabilities in the same way in [Pan et al., 2024]. Subsequently, we compare these edge probabilities to the groundtruth graph structure by the Receiver Operating Characteristic (ROC) curve, and the Area under the ROC Curve (AUC) value. Furthermore, we conduct similar testing on the static version of VAETPP, a SOTA NPPs model with a latent graph. As shown in Fig. 2, the AUC obtained by our approach is better than VAETPP, indicating a more precise capture of the latent relationships between different event types. More detailed information can be found in Appendix E.

6 Conclusion

In this study, we propose a novel approach that leverages multi-dimensional SDEs to enhance the representation capability of the underlying MTPPs. Additionally, we incorporate a graph generator to capture the dependency structure between event types and improve the interaction function between event types in the jump process. Theoretical analysis proves the existence and uniqueness of the solutions of our proposed approaches. Furthermore, experimental results have demonstrated that our method outperforms state-of-theart methods in terms of modeling ability and prediction performance. In the future, we aim to expand the scope of our research by extending our model to encompass the dynamics of non-stationary networks and addressing the scalability of NPPs in larger-scale data, especially concerning the largescale latent graph.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Nos. 62261136549, 62471403, U22A2098, 62271411), the Technological Innovation Team of Shaanxi Province (No. 2025RS-CXTD-009), the International Cooperation Project of Shaanxi Province (No. 2025GH-YBXM-017), the Fundamental Research Funds for the Central Universities (Nos. G2024WD0151, D5000240309).