1 Introduction↩︎

Many application scenarios generate time-stamped event sequences, which can be effectively modeled using temporal point processes (TPPs) [1]. Examples include neural spike train data in neuroscience [2], ask and bid orders in high-frequency financial trading [3], as well as tweets and retweets on social media [4]. These event sequences, composed of asynchronous events, often influence one another and exhibit complex dynamics, making them more challenging to analyze compared to traditional synchronous time series problems. Investigating the underlying dynamic processes of such event sequences not only facilitates the prediction of future events but also helps uncover causal relationships.

In the statistics community, TPPs have a long-standing history of research, with numerous statistical TPP models proposed over the years. Examples include the classic Poisson process [5], Hawkes process [6], and self-correcting process [7], among others. Each of these models is particularly well-suited to specific applications. For instance, the Poisson process was used to model telephone call arrivals, while the Hawkes process, due to its ability to capture self-exciting characteristics, has been widely applied to model earthquakes and aftershocks.

Early TPP models are primarily parametric, requiring explicit specification of the parametric form of the model. However, this imposes limitations on their expressive power. To address these limitations, various nonparametric TPPs have been proposed within the statistics community, including approaches from both the frequentist and Bayesian frameworks, enabling more flexible modeling without the constraints of fixed parametric forms. In recent years, driven by rapid advancements in deep learning, the machine learning community has introduced approaches that combine neural network architectures with TPPs, referred to as neural TPPs. These models further enhance expressive power and are often more intuitive, simpler, and easier to train compared to statistical nonparametric TPPs. Over the past several years, the emergence of large language models (LLMs) has brought transformative changes to the field of artificial intelligence. With their rich contextual understanding and ability to process multimodal data, LLMs offer new possibilities for modeling and analyzing event sequences.

This paper reviews recent advances in TPPs based on Bayesian methods, deep learning, and LLMs, with a focus on model design and parameter estimation. Due to space limitations, we do not aim to cover every method in detail but instead emphasize fundamental principles and core ideas. We also revisit classic applications of TPPs and discuss key challenges and future research directions in the field. This taxonomy of recent advances in TPPs is illustrated in 1. This survey is expected to give comprehensive background knowledge, research trends and technical insights for TPPs.

Several surveys on TPPs in machine learning already exist, such as [8] and [9]. The former summarizes advances in statistical TPPs and neural TPPs, while the latter provides a more detailed overview of neural TPPs. Additionally, surveys in other domains, such as finance, include [10] and [11]. Compared to these works, this paper provides a comprehensive update. In the statistical TPPs domain, we focus on recent progress in Bayesian nonparametric TPPs, which has been largely overlooked—most prior reviews, such as [8], primarily cover frequentist approaches. For neural TPPs, [9] reviews works prior to 2020, but this paper discusses advances from 2020 to 2025. Furthermore, we review the emerging area of LLM-based TPPs, which has gained significant attention in recent years but has yet to be comprehensively reviewed.

=[ rectangle, draw=hidden-draw, rounded corners, text opacity=1, minimum height=6em, minimum width=5em, inner sep=2pt, align=center, fill opacity=.5, line width=1.0pt, ] =[my-box, minimum height=2.8em, fill=hidden-pink!80, text=black, align=left,font=, inner xsep=2pt, inner ysep=4pt, line width=1.2pt, ]

2 Background of TPPs↩︎

We first briefly review some of the core concepts of TPPs. For readers unfamiliar with TPPs, we recommend [12] for a comprehensive introduction.

2.1 Unmarked TPP↩︎

A TPP is a stochastic process that models the occurrence of events over a time window \([0, T]\). A trajectory from a TPP can be represented as an ordered sequence \(\mathcal{T} = (t_1, \ldots, t_{N})\), \(N(t) = \max\{n : t_n \leq t, t \in [0, T]\}\) represents the associated counting process. Here, \(N\) denotes the random number of events within the interval \([0, T]\). TPPs can be defined using different parameterizations. One approach is to specify the distribution of the time intervals between consecutive events. We denote the \(f(t_{n+1} \mid \mathcal{H}_{t_n})\) to be the conditional density function of \(t_{n+1}\) given the history of previous events \(t_1, \ldots, t_n\). In this work, \(\mathcal{H}_{t^-}\) denotes the history of events up to but excluding time \(t\), while \(\mathcal{H}_{t}\) includes whether an event occurs at time \(t\). The conditional density function sequentially specifies the distribution of all timestamps. Consequently, the joint distribution of all events can be factorized as: \[f(t_1, \ldots, t_{N})=\prod_{n=1}^{N} f(t_n\mid \mathcal{H}_{t_{n-1}}). \label{unmark95likelihood}\tag{1}\] A TPP can be defined by specifying the distribution of time intervals. For example, a renewal process assumes that the time intervals are independent and identically distributed (i.i.d.), i.e., \(f(t_n | \mathcal{H}_{t_{n-1}}) = g(t_n - t_{n-1})\), where \(g\) is a probability density function defined on \((0, \infty)\). If we further specify \(g(t_n - t_{n-1})\) to follow an exponential distribution, we obtain a homogeneous Poisson process, where each event occurs independently of the past.

The above approach can directly define some classic point process models. However, in general cases, event occurrences may depend on the entire history, making it less convenient to specify the model using the probability density function of time intervals. Instead, the conditional intensity function provides a more convenient way to describe how the occurrence of an event depends on its history. The conditional intensity function is defined as: \[\begin{align} \lambda^*(t)dt&=\frac{f(t \mid \mathcal{H}_{t_n})dt}{1-F(t \mid \mathcal{H}_{t_n})}\\ &=\frac{P(t_{n+1}\in [t,t+dt]\mid \mathcal{H}_{t_n})}{P(t_{n+1}\notin (t_n,t)\mid \mathcal{H}_{t_n})}\\ &=P(t_{n+1}\in [t,t+dt]\mid t_{n+1}\notin (t_n,t),\mathcal{H}_{t_n})\\ &=P(t_{n+1}\in [t,t+dt]\mid \mathcal{H}_{t^-})\\ &=\mathbb{E}[N([t,t+dt])\mid \mathcal{H}_{t^-}], \end{align} \label{intensity}\tag{2}\] where \(F(t \mid \mathcal{H}_{t_n}) = \int_{t_n}^t f(\tau \mid \mathcal{H}_{t_n}) d\tau\) denotes the cumulative distribution function. Following tradition, we use * to indicate that the conditional intensity function is based on the history. The conditional intensity function has an intuitive interpretation: it specifies the average number of events in a time interval, conditional on the history up to but not including \(t\). It is worth noting that the history \(\mathcal{H}_{t_n}\) in the conditional density function differs from the history \(\mathcal{H}_{t^-}\) in the conditional intensity function. This subtle distinction is often overlooked in many TPP works.

The conditional intensity function and the conditional density function are one-to-one⁵. This can be easily proven by inverting 2 to express the conditional density function in terms of the conditional intensity function: \[\begin{align} F(t \mid \mathcal{H}_{t_n})&=1-\exp\left(-\int_{t_n}^t \lambda^*(\tau)d\tau\right),\\ f(t \mid \mathcal{H}_{t_n})&=\lambda^*(t)\exp\left(-\int_{t_n}^t \lambda^*(\tau)d\tau\right). \end{align} \label{con95density}\tag{3}\]

This means we can define new TPP models by directly specifying a particular form of the conditional intensity function. For example, specifying a constant intensity defines a homogeneous Poisson process, while specifying a time-varying intensity function \(\lambda^*(t)=\lambda(t)\) defines an inhomogeneous Poisson process. We can also define a Hawkes process by specifying a conditional intensity function: \[\lambda^*(t) = \mu + \sum_{t_n < t} \phi(t - t_n), \label{hawkes}\tag{4}\] where \(\mu > 0\) is the baseline intensity, and \(\phi(\cdot): \mathbb{R}^+ \to \mathbb{R}^+\) is the triggering function⁶. The summation of influences from past events increases the likelihood of future events, making it suitable for modeling self-exciting effects. While many other forms of TPPs exist, we primarily focus on the Poisson process (history-independent) and the Hawkes process (history-dependent) in the following due to their widespread use. An illustration of the unmarked temporal point process, along with its conditional density function and conditional intensity function, is shown in 2.

2.2 Marked TPPs↩︎

The above discussion focuses on the unmarked TPP, also known as the univariate TPP. However, TPPs can be extended to the marked case, represented as a time-ordered marked sequence \(\mathcal{T} = ((t_1, k_1), \ldots, (t_{N}, k_{N}))\) over a time window \([0, T]\), where \(k_n\) is the mark of the \(n\)-th event. The mark space can be continuous or discrete. In practice, discrete marks are more common, often referred to as multivariate TPPs. Similar to the unmarked case, marked TPPs can also be described using the conditional density function: \[f((t_1, k_1), \ldots, (t_{N}, k_{N})) = \prod_{n=1}^{N} f(t_n, k_n \mid \mathcal{H}_{t_{n-1}}),\] where \(f(t, k \mid \mathcal{H}_{t_n})\) is the joint density of time and mark, conditional on history. The history \(\mathcal{H}_{t_n}\) now includes information about both the times and marks of past events.

Similarly, when event occurrences depend on the entire history, specifying the model using the probability density function becomes less convenient. In such cases, the conditional intensity function offers a more practical and expressive way to capture the dependence of events on historical information. The conditional intensity function is defined as: \[\lambda^*(t, k) dtdk = \frac{f(t, k \mid \mathcal{H}_{t_n}) dtdk}{1 - F(t \mid \mathcal{H}_{t_n})}=\mathbb{E}[N(dt\times dk)\mid \mathcal{H}_{t^-}],\] where \(F(t\mid\mathcal{H}_{t_n}) = \int_{t_n}^t \int_k f(\tau,k\mid\mathcal{H}_{t_n}) dk d\tau\) is the conditional cumulative distribution function and the mark is marginalized out. The conditional intensity function specifies the average number of events with a mark \(k\) in a time interval, conditional on the history up to but not including \(t\).

We can also define a marked TPP by specifying a particular conditional intensity function. A classic example is the multivariate Hawkes processes, where the mark \(k \in \{1, \ldots, K\}\) represents the event type. The conditional intensity function of multivariate Hawkes processes is given by: \[\lambda^*(t, k) = \mu_k + \sum_{t_n < t} \phi_{k, k_n}(t - t_n), \label{mhp}\tag{5}\] where \(\mu_k\) represents the baseline intensity for event type \(k\), and \(\phi_{k, k^\prime}\) captures the triggering effects from events of type \(k^\prime\) on type \(k\).

2.3 Inference↩︎

There are various methods for estimating the parameters of TPPs. The most common method is maximum likelihood estimation (MLE). In this section, we introduce MLE and Bayesian inference for parametric TPPs. Assume we observe a trajectory of a marked TPP \(\mathcal{T} = ((t_1, k_1), \ldots, (t_{N}, k_{N}))\) over the time window \([0, T]\). The unmarked case corresponds to the situation where there is only a single mark. Assume the marked TPP is specified by a parametric conditional intensity function \(\lambda_\theta^*(t, k)\). Then, the likelihood function is given by: \[\begin{align} f(\mathcal{T};\theta)=\prod_{n=1}^{N}\lambda_\theta^*(t_n,k_n)\exp\left(-\int_0^T \lambda_\theta^*(t)dt\right), \end{align} \label{mark95likelihood}\tag{6}\] where \(\lambda_\theta^*(t)=\int\lambda_\theta^*(t, k)dk\) is the ground intensity. The proof of 6 is straightforward. Simply substitute 3 into 1 to verify it for the unmarked case. The proof for the marked case follows a similar procedure. We can use numerical methods to maximize the log-likelihood to obtain parameter estimates.

MLE is a widely used parameter estimation method in the frequentist framework. It offers several advantages, such as consistency, asymptotic normality, and asymptotic efficiency. However, as a point estimation method, MLE cannot capture model uncertainty, which limits its applicability in high-stakes domains where understanding uncertainty is crucial. To address this issue, the Bayesian framework has been incorporated into TPPs [13]. In Bayesian TPPs, we impose suitable priors on the model parameters and then compute their corresponding posterior distribution, equipping the model with the ability to quantify uncertainty. Specifically, a Bayesian TPP is formally expressed as: \[f(\theta\mid \mathcal{T})=\frac{f(\mathcal{T}\mid \theta)f(\theta)}{\int f(\mathcal{T}\mid \theta)f(\theta)d\theta},\] where \(f(\mathcal{T}\mid \theta)\) is the likelihood in 6 , \(f(\theta)\) is the prior on model parameters, the denominator is the marginal likelihood, and \(f(\theta\mid \mathcal{T})\) is the posterior distribution of model parameters. In general, the inference for Bayesian TPPs is more challenging than for frequentist TPPs because the TPP likelihood is not conjugate to any prior. This means that the posterior does not have an analytical expression and can only be obtained through approximation methods, such as Markov chain Monte Carlo (MCMC) [14], [15], variational inference [16], [17], and Laplace approximation [18], [19], among others.

3 Bayesian Nonparametric TPPs↩︎

Early work on TPPs was limited to parametric models, whether in the frequentist or Bayesian framework. These methods rely heavily on model assumptions, making them inflexible and often performing poorly on complex datasets. To address this limitation, many studies have proposed nonparametric methods. [8] provides a comprehensive description of frequentist nonparametric methods, while this paper focuses on Bayesian nonparametric TPPs, which not only enhance model flexibility but also incorporate the ability to quantify model uncertainty.

3.1 Bayesian Nonparametric Poisson Process↩︎

Bayesian nonparametric TPPs do not parameterize the intensity function into a fixed form. Instead, they treat the intensity function itself as a model parameter with infinite dimensions and impose suitable prior on it. For instance, in the case of an inhomogeneous Poisson process, the intensity function \(\lambda(t)\) is treated as the parameter, and a prior \(f(\lambda(t))\) is placed on it. The goal is then to compute the posterior distribution. This prior, being a distribution over functions, is commonly modeled using a Gaussian process (GP). Consequently, the Bayesian nonparametric framework is formulated as follows: \[f(g(t) \mid \mathcal{T})=\frac{f(\mathcal{T}\mid \lambda(t)=l \circ g(t))\mathcal{GP}(g(t))}{\int f(\mathcal{T}\mid \lambda(t)=l \circ g(t))\mathcal{GP}(g(t))dg}, \label{inference}\tag{7}\] where \(l(\cdot): \mathbb{R} \to \mathbb{R}^+\) is a link function ensuring the intensity function is non-negative, e.g., exponential (log Gaussian Cox process [20]), scaled sigmoid (sigmoidal Gaussian Cox process [21]), square (permanental process [22]), etc. It is worth noting that computing 7 is highly challenging, as the posterior is doubly intractable due to an intractable integral over \(t\) in the numerator and another over \(g\) in the denominator. This is a well-known problem in the field of Bayesian nonparametric TPPs [23].

Many methods have been proposed to solve 7 . Some studies focus on utilizing MCMC methods. [21] proposed an MCMC inference method for a Poisson process with a sigmoidal GP prior. The core idea is to incorporate latent thinned points to make the posterior tractable. However, this method scales cubically with the number of data and thinned points. [24] extended this approach to multiple dependent Cox processes using multi-output Gaussian processes, and also derived an MCMC sampler for performing inference. Later, [25] leveraged inducing points, a common technique in GP for reducing time complexity [26], to derive an MCMC sampler that reduces the computational cost to linear w.r.t. the number of data points.

Some studies focus on methods based on the Laplace approximation. [27] proposed a Laplace approximation method that utilizes GP on a fixed grid to estimate the intensity function. However, strictly speaking, grid-based approaches are no longer fully nonparametric. [28] combined Kronecker methods with the Laplace approximation to enable scalable inference. [29] proposed a fast Laplace approximation method relying on the Mercer decomposition of the GP kernel. However, its tractability is limited to standard kernels such as the squared exponential kernel. To address this limitation, [30] introduced an alternative fast Laplace approximation leveraging the spectral representation of kernels. This approach retains tractability while accommodating a broader range of stationary kernels. Furthermore, [31] extended this method to non-stationary kernels by leveraging sparse spectral representations to overcome the limitations of stationary kernels. This approach provides a low-rank approximation of the kernel, effectively reducing computational complexity from cubic to linear.

Another common approach is variational inference. [16] introduced the first fully variational inference scheme for permanental processes. However, similar to [29], its tractability is limited to certain standard types of kernels. [32] further extended the method in [16] to a multitask point process model, leveraging information from all tasks via a hierarchical GP. [33] expanded the approach in [16] to utilize the Fourier representation of the GP, enabling the use of more general stationary kernels. [34] presented a novel tractable representation of the likelihood through augmentation with a superposition of Poisson processes. This perspective enabled a structured variational approximation that captured dependencies across variables in the model. The method avoided discretization of the domain, did not require numerical integration over the input space, and was not limited to GPs with squared exponential kernels.

It is worth noting that, [35] introduced the data augmentation technique based on Pólya-Gamma variables to the field of Bayesian nonparametric TPPs. This technique is an improvement and extension of the method proposed by [21]. The method augments not only thinned points but also Pólya-Gamma latent variables for all data and thinned points in the likelihood. This enables the augmented likelihood to be conditionally conjugate to the GP prior. By leveraging the conditionally conjugacy, we can derive fully analytical Gibbs sampler, EM algorithm, and mean-field variational inference method. This method was later extended to jointly model multiple heterogeneous and correlated tasks—such as classification, regression, and point processes—using multi-output Gaussian processes. This extension facilitated information sharing across heterogeneous tasks while enabling nonparametric parameter estimation [36].

The works discussed above primarily use GP as prior. However, other forms of Bayesian nonparametric priors are also possible. For example, [37], [38] used a Dirichlet process mixture of Beta densities as a prior for the normalized intensity function of a Poisson process. Bayesian nonparametric TPPs based on Dirichlet process mixtures and those based on Gaussian processes represent two orthogonal modeling paradigms, both capable of achieving Bayesian nonparametric inference. Inference for Dirichlet process mixture-based TPPs typically relies on specialized MCMC or variational inference techniques developed within the Dirichlet process framework. A detailed discussion of these methods is beyond the scope of this paper.

3.2 Bayesian Nonparametric Hawkes Process↩︎

For the Hawkes process, as shown in 4 , the conditional intensity function consists of the baseline intensity \(\mu(\cdot)\)⁷ and the triggering function \(\phi(\cdot)\). Therefore, the Bayesian nonparametric Hawkes process typically places GP priors on both \(\mu(\cdot)\) and \(\phi(\cdot)\) (we take the unmarked case as an example): \[\begin{align} &f(g(t),h(\tau) \mid \mathcal{T}) \propto \\ &f(\mathcal{T}\mid \mu(t)=l \circ g(t),\phi(\tau)=l \circ h(\tau))\mathcal{GP}(g(t))\mathcal{GP}(h(\tau)), \end{align} \label{hawkes95inference}\tag{8}\] where \(l(\cdot)\) is a link function ensuring \(\mu(t)\) and \(\phi(\tau)\) are non-negative, similar to 7 . Computing 8 is more challenging than 7 because in the likelihood of Hawkes process, \(\mu(t)\) and \(\phi(\tau)\) are coupled together, which significantly complicates the inference process.

To address this issue, a common approach is to augment a branching latent variable into the Hawkes process likelihood to indicate whether each event is triggered by itself via the baseline intensity or by a previous event via the triggering function [39]. The branching variable \(\mathbf{X}\) is a lower triangular matrix with Bernoulli entries, where \(x_{nm}\) indicates whether the \(n\)-th event is triggered by itself or a previous event \(m\): \[\begin{align} &x_{nn}=\left\{ \begin{aligned} &1 \;\;\text{if event n is a background event},\\ &0 \;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{otherwise},\\ \end{aligned} \right.\\ &x_{nm}=\left\{ \begin{align} &1 \;\;\text{if event n is caused by event m},\\ &0 \;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{otherwise}.\\ \end{align} \right. \end{align}\] After augmenting the branching latent variable, the joint likelihood is expressed as: \[\begin{align} &f(\mathcal{T},\mathbf{X} \mid \mu(t),\phi(\tau))=\underbrace{\prod_{n=1}^{N}\mu(t_n)^{x_{nn}}\exp{\left(-\int_0^T\mu(t)dt\right)}}_{\text{baseline intensity part}}\cdot\\ &\underbrace{\prod_{n=2}^{N}\prod_{m=1}^{n-1}\phi(t_n-t_m)^{x_{nm}}\prod_{n=1}^{N}\exp{\left(-\int_0^{T_\phi}\phi(\tau)d\tau\right)}}_{\text{triggering function part}}, \end{align}\] where the support of triggering function is assumed to be \([0,T_\phi]\). If the branching variable \(\mathbf{X}\) is marginalized out, we obtain the original likelihood.

It is clear that, after introducing the branching variable, the joint likelihood factorizes into two independent components: one corresponding to the baseline intensity and the other to the triggering function. These two components are linked through the branching variable \(\mathbf{X}\). To the best of our knowledge, [39] was the first to identify this structure and subsequently proposed an EM algorithm that leverages it: the E-step computes the posterior distribution of the branching variable, while the M-step estimates the parameters of both the baseline intensity and the triggering function. Later, [40] extended this approach to the frequentist nonparametric setting by treating the baseline intensity and the triggering function as two flexible, unconstrained functions. They imposed Good’s roughness penalty [41] on both functions and derived the solutions using the Euler–Lagrange equation. [42] further extended the method in [40] to multivariate Hawkes processes by assuming that the triggering functions in the multivariate setting are linear combinations of a set of basis functions. Each basis function was then estimated using the Euler–Lagrange equation. There also exist frequentist nonparametric approaches that do not rely on the branching variable. For example, [43] proposed an estimation method based on solving a Wiener–Hopf equation that relates the triggering function to the second-order statistics. Additionally, [44], [45] attempted to estimate the triggering function by minimizing a quadratic contrast function using a grid-based representation of the triggering function. [46] extended the approach of [47], which formulated Poisson process intensity estimation within a reproducing kernel Hilbert space (RKHS) framework, to the more complex setting of Hawkes processes.

In the Bayesian nonparametric Hawkes process setting, a similar strategy can be adopted. By introducing the branching variable, the Hawkes likelihood factorizes into two independent components, each of which can be viewed as an independent Poisson process. Since we place independent GP priors on these two components, we can directly apply the methods discussed in 3.1 to compute the posteriors of \(\mu(t)\) and \(\phi(\tau)\). This naturally leads to an iterative algorithm where, at each iteration, the posterior of \(\mathbf{X}\) is used to update the posteriors of \(\mu(t)\) and \(\phi(\tau)\), which are then used to update the posterior of \(\mathbf{X}\) in turn.

In recent years, many GP-based Bayesian nonparametric Hawkes process studies have adopted this iterative framework. [48] derived a Gibbs sampler and a maximum a posteriori (MAP) EM algorithm to estimate a nonparametric triggering function. [49] extended the variational inference method from [16] to the Hawkes process for estimating a nonparametric triggering function. Further, [50] applied variational inference to simultaneously estimate the nonparametric baseline intensity and triggering function. [51], [52] extended the data augmentation method based on Pólya-Gamma variables from [35] to the Hawkes process. This work derived fully analytical Gibbs sampler, EM algorithm, and mean-field variational inference method. This approach was subsequently extended to nonlinear Hawkes processes to account for excitation and inhibition effects between interacting variables [53]. [54] analyzed the theoretical properties of this approach; specifically, it provided concentration rates for the posterior distribution of the parameters under mild assumptions on both the prior and the model, and established consistency guarantees for the inferred Granger causal graph.

The works discussed above primarily use GP as prior. However, other forms of Bayesian nonparametric priors are also possible. For instance, [55]–[57] extended the Hawkes process to the Bayesian nonparametric setting using Dirichlet process mixture priors. Building on this line of research, [58] proposed an online inference method using sequential Monte Carlo techniques. Additionally, [59] introduced priors based on piecewise constant functions with either regular or random partitions, as well as priors defined via mixtures of Beta distributions. A comprehensive discussion of alternative approaches lies beyond the scope of this paper.

4 Neural TPPs↩︎

Benefiting from the rapid development of deep learning, another way to enhance the flexibility of TPPs is by using deep models to model TPPs. Compared to frequentist/Bayesian nonparametric TPPs, neural TPPs offer more intuitive and straightforward modeling and parameter estimation. [9] provided a comprehensive review of neural TPPs, but it focuses on work prior to 2020. This paper focuses more on the latest advancements from 2020 to 2024. In the following, we categorize neural TPPs into three major types and introduce each in detail. A schematic illustration of these models is shown in 3.

4.1 Recurrent Neural TPPs↩︎

The earliest work on neural TPPs can be traced back to [60], which was the first to use a recurrent neural network (RNN) to model TPPs. In that work, each event in \(\{(t_n, k_n)\}_{n=1}^N\) is passed through an embedding layer to obtain a more compact representation \(\mathbf{y}_n\). At each event location, a history embedding \(\mathbf{h}_n\) is designed to capture historical information. Then, when a new event occurs, the history embedding is updated based on the feature of the current event: \[\mathbf{h}_n = \text{Update}(\mathbf{h}_{n-1}, \mathbf{y}_n).\] Finally, the history embedding \(\mathbf{h}_n\) is used to parameterize the conditional distribution of the next event. In [60], \(\mathbf{h}_n\) is used to represent the conditional intensity function after \(t_n\): \[\lambda^*(t, k) = \exp{\left(\mathbf{v}_k^\top\mathbf{h}_n + w_k(t - t_n) + b_k\right)}, \label{rnn95intensity}\tag{9}\] where \(\mathbf{v}_k\), \(w_k\), and \(b_k\) are learnable parameters for \(k\)-th mark. The first term captures the influence of past events through the history embedding, the second term is an extrapolation component that models the intensity at time \(t\), and the third term is a bias. The exponential function ensures that the conditional intensity function remains positive. As discussed in 2, although the conditional intensity function is a common parameterization, other formulations for characterizing the conditional distribution of the next event are also feasible.

Due to the gradient vanishing or exploding problems of traditional RNNs, [60] was unable to model long-range dependencies. To address this issue, [61] proposed a long short-term memory [62] (LSTM)-based TPP model, which mitigates the shortcomings of traditional RNNs and achieves improved performance. [63] introduced a dual-LSTM framework, where one LSTM models the TPP and the other models a time-series covariate. The history embeddings from both networks are then fused to predict the next event. By leveraging the additional information from covariates, TPPs can achieve improved prediction performance. [64] further extended the LSTM-based point process model to the spatio-temporal domain. [65] used RNNs to encode history, but modeled the cumulative conditional intensity instead of the conditional intensity itself to avoid costly numerical integration during MLE. This enables efficient training via differentiation rather than integration, as discussed in 4.4. [66] proposed UniPoint, an RNN-based model that serves as a universal approximator for point process intensity functions. They theoretically proved that RNNs can approximate any valid intensity by leveraging the Stone-Weierstrass theorem. [67] addressed missing events by using two RNNs to model the generative processes of observed and missing events, with the latter treated as latent variables. They proposed an unsupervised training method based on variational inference to jointly learn both models. [68] established excess risk bounds for RNN-based TPPs under various TPP settings. They showed that a four-layer RNN-TPP can achieve vanishing generalization error.

The main advantage of recurrent models lies in their computational efficiency. During the prediction phase, once the history embedding \(\mathbf{h}_n\) is obtained, the model can predict the next event with constant time complexity \(O(1)\), regardless of the sequence length. This makes recurrent models particularly suitable for online or streaming settings where fast real-time prediction is essential. Moreover, during training, the time complexity scales linearly with the number of events, i.e., \(O(N)\), since the history embedding is updated in an iterative, step-by-step manner along the event sequence. However, recurrent models also suffer from several notable limitations. First, due to their inherently sequential architecture, they cannot be efficiently parallelized during training, which significantly slows down the training process, especially for long event sequences. Second, traditional RNN-based models struggle to capture long-range dependencies due to issues such as vanishing or exploding gradients, which can degrade the model’s ability to learn complex temporal dynamics over extended horizons. Although architectures such as LSTMs and GRUs [69] partially mitigate these issues, they do not fundamentally resolve them. These limitations have motivated the development of alternative architectures, such as attention-based and Transformer-style models, which offer better support for parallelism and more effective modeling of long-range dependencies. In 1, we summarize the advantages and disadvantages of recurrent models and autoregressive models for comparison.

It is worth noting that in recent years, several powerful recurrent architectures have been proposed in the field of sequential models, such as the Receptance Weighted Key Value (RWKV) [70], the Structured State Space Sequence (S4) [71], and Mamba [72]. These works aim to design efficient recurrent architectures to replace Transformers, as Transformers have a high time complexity of \(O(N^2)\) in training and \(O(N)\) in prediction. All these models are recurrent architectures, so they share the same advantages as RNNs, such as \(O(1)\) time complexity for prediction and \(O(N)\) time complexity for training. However, their key distinction from RNN lies in their ability to support parallelized training and capture long-range dependencies. Integrating these novel and efficient recurrent architectures with TPPs is an important future research direction. This integration could significantly enhance the scalability of neural TPPs for training and prediction on large-scale datasets. Currently, works in this area are limited. [73] and [74] explored the idea of combining Mamba or deep state space models with TPPs, offering a promising path forward for scalable TPP modeling.

4.2 Autoregressive Neural TPPs↩︎

As stated in 4.1, due to the limitations of RNNs, such as the inability to support parallel training and capture long-range dependencies, a large number of studies since 2020 have explored using Transformer architectures to model TPPs. The earliest works include [75] and [76], which share similar ideas but differ slightly in certain details.

Each event in the sequence \(\{(t_n, k_n)\}_{n=1}^N\) is encoded as a feature vector \(\mathbf{y}_n \in \mathbb{R}^M\), combining both temporal and mark information. We adopt sinusoidal positional encodings for the temporal component, defined as: \[z_j(t_n) = \begin{cases} \cos\left(t_n / 10000^{\frac{j-1}{M}} \right), & \text{if } j \text{ is odd},\\ \sin\left(t_n / 10000^{\frac{j}{M}} \right), & \text{if } j \text{ is even}, \end{cases}\] where \(z_j(t_n)\) denotes the \(j\)-th entry of the time embedding vector \(\mathbf{z}(t_n) \in \mathbb{R}^M\), and \(j = 0, \ldots, M-1\). The complete time embedding matrix is denoted by: \[\mathbf{Z} = [\mathbf{z}(t_1), \ldots, \mathbf{z}(t_N)]^\top \in \mathbb{R}^{N \times M}.\] Let \(\mathbf{U} \in \mathbb{R}^{M \times K}\) be a learnable mark embedding matrix, where \(K\) is the total number of event types. The mark \(k_n\) is first converted into a one-hot vector \(\mathbf{k}_n \in \mathbb{R}^K\), and its embedding is given by: \[\mathbf{e}(k_n) = \mathbf{U} \mathbf{k}_n \in \mathbb{R}^M.\] Collecting all event mark embeddings yields: \[\mathbf{E} = [\mathbf{e}(k_1), \ldots, \mathbf{e}(k_N)]^\top \in \mathbb{R}^{N \times M}.\] The final embedding for each event is obtained by summing its temporal and mark embeddings: \[\mathbf{Y} = \mathbf{Z} + \mathbf{E} \in \mathbb{R}^{N \times M},\] where each row \(\mathbf{y}_n\) in \(\mathbf{Y}\) represents the complete embedding of the \(n\)-th event in the sequence. The matrix \(\mathbf{Y}\) is then multiplied by corresponding weight matrices to compute the query, key, and value matrices: \(\mathbf{Q} = \mathbf{Y}\mathbf{W}_Q\in\mathbb{R}^{N\times M_K}, \quad \mathbf{K} = \mathbf{Y}\mathbf{W}_K\in\mathbb{R}^{N\times M_K}, \quad \mathbf{V} = \mathbf{Y}\mathbf{W}_V\in\mathbb{R}^{N\times M_V}.\) Finally, the attention score is computed as: \[\mathbf{S}=\text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^\top}{\sqrt{M_K}}\right)\mathbf{V}. \label{attention}\tag{10}\] To ensure causality and prevent future events from affecting past events, it applies a mask to the upper triangular entries of \(\mathbf{Q}\mathbf{K}^\top\).

The attention score is then used to generate the history embedding \(\mathbf{h}_n\), which serves as the history representation for parameterizing the conditional distribution of the next event. In [75], [76], \(\mathbf{h}_n\) is used to represent the conditional intensity function after \(t_n\): \[\lambda^*(t, k) = \text{softplus}\left(\mathbf{v}_k^\top \mathbf{h}_n + \frac{w_k(t - t_n)}{t_n} + b_k\right), \label{attention95intensity}\tag{11}\] where \(\mathbf{v}_k\), \(w_k\), and \(b_k\) are learnable parameters for \(k\)-th mark. However, other parameterizations for characterizing the conditional distribution of the next event are also feasible. For example, [77] used the history embedding \(\mathbf{h}_n\) obtained from the Transformer to model the conditional density function of the next event \(f(t \mid \mathbf{h}_n)\). This approach makes sampling the next event more efficient.

The advantages and disadvantages of autoregressive models stand in contrast to those of recurrent models. A key strength of autoregressive models lies in their ability to support parallel training and effectively capture long-range dependencies via self-attention mechanisms. However, these benefits come at the cost of computational efficiency. During training, the time and memory complexity scales quadratically with the sequence length, i.e., \(\mathcal{O}(N^2)\), due to the self-attention computation across all event pairs. During prediction, autoregressive models typically require all previous hidden states to compute the next event, resulting in a time complexity of \(\mathcal{O}(N)\). While key-value (KV) caching techniques can reduce the prediction time to constant time complexity \(\mathcal{O}(1)\) by storing intermediate representations, they incur significant memory overhead, which can become prohibitive for long sequences or memory-constrained environments. In contrast, recurrent models maintain a hidden state that is updated incrementally, leading to linear time complexity \(\mathcal{O}(N)\) for training and constant time \(\mathcal{O}(1)\) for prediction, but at the cost of sequential processing and limited ability to model long-range dependencies. In 1, we provide a comparison of the advantages and disadvantages of autoregressive and recurrent models.

Between 2020 and 2024, numerous studies have proposed improvements to Transformer TPP models. For instance, [78] modified the computation of the attention score in 10 , replacing the commonly used dot-product attention with a flexible non-linear attention score based on Fourier kernels, enabling the capture of more complex event similarities. [79] extended Transformer-based point process models to the spatio-temporal domain by first using a Transformer to encode the event history. The resulting history embedding is then mapped to a latent stochastic process, which is subsequently used to generate a set of temporal and spatial kernels. These kernels are linearly combined to construct the spatio-temporal conditional intensity function. [80] further improved the attention mechanism proposed in [75], [76] by introducing future-event-specific query vectors that incorporate continuous positional encodings of time \(t\). As \(t\) increases, the attention weights over past events vary smoothly, enabling a more adaptive and temporally aware modeling of historical influence. [81] proposed a sparse Transformer TPP model based on a sliding window mechanism to reduce the quadratic time and space complexity of Transformers. [82] introduced improvements in the combination of time and mark embeddings, the computation of \(\mathbf{Q}\) and \(\mathbf{K}\) matrices, and the modeling of the conditional intensity function. These enhancements allowed the Transformer Hawkes process to perfectly align with statistical nonlinear Hawkes processes, thereby improving its interpretability. [83] combined Transformer TPPs with federated learning, enabling collaborative learning from large amounts of distributed event sequence data.

4.3 Differential Equation-based Neural TPPs↩︎

Recurrent and autoregressive neural TPPs share a fundamental limitation: as discrete-time models, they can only compute the history embedding \(\mathbf{h}_n\) and the conditional intensity function \(\lambda^*(t_n)\) at discrete event times. However, they cannot directly characterize the conditional intensity function over the continuous intervals between events. To address this issue, many studies introduce extrapolation terms to approximate the intensity function over such intervals. For example, in 9 , the term \(w_k(t - t_n)\), and in 11 , the term \(w_k(t - t_n)/t_n\), are both designed to serve as extrapolation mechanisms. However, these extrapolation components adopt fixed parametric forms, which limit the expressiveness of the conditional intensity function over event intervals. Specifically, the extrapolation behavior in 11 is approximately linear with respect to \(t\), as the \(\text{softplus}\) function closely resembles a ReLU. This limitation is visually demonstrated in [82], where the authors plot the conditional intensity function over intervals to highlight this issue.

Differential equation-based TPPs represent another line of research. These models, being continuous-time, can model the conditional intensity function over continuous time, thus avoiding the above issue. Specifically, these models utilize differential equations (often stochastic differential equations (SDEs)) to describe a history-dependent left-continuous hidden state \(\mathbf{h}^*(t)\) over \([0, T]\) with an initial state \(\mathbf{h}^*(0)\). For instance, in [84], the SDE is defined as: \[d\mathbf{h}^*(t)=\text{NN}_{\theta_1}(\mathbf{h}^*(t),t)dt+\text{NN}_{\theta_2}(\mathbf{h}^*(t),t)dN(t), \label{sde}\tag{12}\] where we take the unmarked case as an example. The functions \(\text{NN}_{\theta_1}\) and \(\text{NN}_{\theta_2}\) are two neural networks that govern the flow and jump of \(\mathbf{h}^*(t)\), respectively. \(N(t)\) is the counting process that records the number of events up to time \(t\). When no event occurs, \(\mathbf{h}^*(t)\) evolves smoothly according to \(\text{NN}_{\theta_1}\). When an event occurs, \(\mathbf{h}^*(t)\) undergoes a jump at the event’s timestamp, governed by \(\text{NN}_{\theta_2}\). Then, the history-dependent left-continuous hidden state \(\mathbf{h}^*(t)\) is used to define the conditional intensity function: \[\lambda^*(t) = \text{NN}_{\theta_3}(\mathbf{h}^*(t)), \label{sde95intensity}\tag{13}\] where \(\text{NN}_{\theta_3}\) is another neural network ensuring a non-negative output. Finally, we use MLE to estimate the model parameters. Clearly, since \(\mathbf{h}^*(t)\) can vary flexibly over the event intervals—being parameterized by \(\text{NN}_{\theta_1}\)—the resulting conditional intensity function \(\lambda^*(t)\) also exhibits flexible dynamics between events. This effectively overcomes the limited expressiveness of recurrent and autoregressive neural TPPs in modeling the conditional intensity function over event intervals.

The earliest work, to the best of our knowledge, that integrates differential equations with point processes is [85], which considered a simple inhomogeneous Poisson process whose intensity function evolves according to an ordinary differential equation (ODE). However, this model does not account for the discontinuities in the conditional intensity function caused by historical events. To address this limitation, [84] proposed a method based on SDE, as shown in 12 13 , to model the conditional intensity of history-dependent TPPs. This approach allows the conditional intensity function to exhibit jumps at the occurrence of events. A contemporaneous work, [86], combined ODEs with RNNs. In this framework, the hidden state \(\mathbf{h}^*(t)\) evolves smoothly according to an ODE when no event occurs, and undergoes a jump at the event time, governed by the RNN update. This enables the model to capture jumps in the conditional intensity function at event times. In fact, [86] shares strong similarities with [61], but is more general. Specifically, [61] assumes that \(\mathbf{h}^*(t)\) evolves across event intervals following an exponential law, whereas [86] allows \(\mathbf{h}^*(t)\) to evolve flexibly over event intervals according to arbitrary ODEs. The ODE-based approach was later extended to spatio-temporal TPPs by [87]. [88] adopted a SDE-based approach to directly model the conditional intensity function of TPPs. Additionally, they used an SDE to model users’ opinions, enabling the incorporation of user feedback into the TPP framework. This results in a closed-loop model that jointly captures the dynamics of users’ opinions and event generation. More recently, [89] proposed a novel SDE-based method for modeling TPPs. Instead of modeling a latent hidden state via an SDE and mapping it to the conditional intensity function, their approach directly models the conditional intensity function with an SDE, and further provides a theoretical analysis on the existence and uniqueness of its solution.

Although differential equation-based neural TPPs offer greater flexibility than recurrent and autoregressive models in characterizing the conditional intensity function over event intervals, they suffer from significant computational inefficiencies during both training and sampling. During training, when using MLE, the model requires numerical integration to compute the integral of the conditional intensity function. This necessitates evaluating the intensity at multiple time points, which can only be obtained by solving the ODE or SDE defined in 12 using a numerical solver. As a result, training is considerably slower compared to recurrent or autoregressive neural TPPs. Similarly, during the sampling phase—such as when applying the thinning algorithm—intensity values at various time points must also be computed via ODE or SDE solvers, further increasing the computational cost and slowing down the sampling process.

4.4 Parametrization of Neural TPPs↩︎

Both recurrent and autoregressive neural TPPs extract a history embedding from past event information and use it to model the conditional distribution of the next event. As discussed in 2, there are multiple ways to represent the conditional distribution of the next event, such as the conditional density function \(f(t \mid \mathcal{H}_{t_n})\) in 1 , the cumulative distribution function \(F(t \mid \mathcal{H}_{t_n})\) in 2 , the conditional intensity function \(\lambda^*(t)\) in 2 , and the cumulative intensity function \(\Lambda^*(t) = \int_0^t \lambda^*(\tau)d\tau\). All of these parameterizations are equivalent and can be used to characterize the conditional distribution of the next event, since, as proved in 3 , the conditional intensity function and the conditional density function are in a one-to-one correspondence. A schematic illustration of the four common parameterizations of TPPs is shown in 4.

Most of the works discussed above adopt the conditional intensity function as the primary parameterization. The main advantage of this approach lies in its conceptual and implementation simplicity—the model only needs to ensure that the output is non-negative, which can be easily achieved using functions like ReLU, softplus, or exponential mappings. However, a notable drawback arises during MLE: the log-likelihood requires evaluating the integral of the conditional intensity function over the observation window. In most practical cases, this integral does not admit a closed-form solution and must be approximated using numerical integration techniques, such as Monte Carlo sampling or quadrature. This introduces additional computational overhead and can compromise both the estimation accuracy and the overall training efficiency, especially for complex or high-dimensional models.

In recent years, several studies have explored alternative parameterizations for TPPs. For example, [90] and [77] modeled the conditional density function \(f(t \mid \mathcal{H}_{t_n})\) directly as a mixture of log-normal distributions, using history embeddings as inputs. The key difference between the two lies in how the history embeddings are extracted: [90] used a RNN, while [77] employed a Transformer-based encoder. One major advantage of modeling the conditional density directly is that it eliminates the need to compute the integral of the conditional intensity function during MLE, earning it the name “intensity-free” modeling. Another benefit is efficient sampling: given the history, the next event can be sampled directly from a mixture of log-normal distributions, which admits a simple and tractable form. [91] modeled the inverse of the cumulative distribution function \(F^{-1}(t \mid \mathcal{H}_{t_n})\) using the history embedding with monotonic rational-quadratic splines. This method also avoids numerical integration and supports efficient sampling. [65], [92], [93] proposed modeling the cumulative intensity function \(\Lambda^*(t)\) conditioned on history embeddings using monotonic neural networks or splines. This enables a reformulation of the log-likelihood in 6 , which, for the unmarked case, becomes: \[\log f(\mathcal{T};\theta) = \sum_{n=1}^{N} \log\frac{d}{dt}\Lambda_{\theta}^*(t_n^-) - \Lambda_\theta^*(T), \label{CCIF-likelihood}\tag{14}\] where \(t_n^-\) denotes the left-hand limit of the derivative at \(t_n\). By modeling the cumulative intensity function directly, this parameterization transforms the integral in the log-likelihood into a derivative, allowing for efficient computation using automatic differentiation. Consequently, it eliminates the need for numerical integration during MLE training, improving both accuracy and efficiency.

5 LLM-based TPPs↩︎

LLM-based TPPs are essentially neural TPPs. However, given the significant impact of LLMs in recent years, we dedicate a separate section to this topic. Recent studies have explored integrating LLMs with TPPs along with benchmarks and models for LLM-inspired techniques, such as retrieval-augmented TPPs. A schematic illustration of LLM-based TPPs is shown in 5.

5.1 LLM-inspired TPPs↩︎

[97] proposed PromptTPP, which incorporates continual learning into neural TPPs, enabling efficient adaptation to streaming event data. This framework was inspired by recent advances in prompt learning methods for LLMs where learnable prompt parameters, dubbed temporal prompts, are prepended to the input tokens to improve the model prediction performance. The key innovation lies in integrating a continuous-time retrieval prompt pool with the base TPP model. The prompts are small learnable parameters stored in memory space and jointly optimized with the base TPP, ensuring sequential learning without buffering past examples or task-specific attributes. This approach addresses critical challenges in real-world applications including privacy constraints, memory limitations, and distribution shifts over time. The framework demonstrates consistent state-of-the-art performance across three real user behavior datasets, showcasing its effectiveness in handling streaming event sequences where traditional TPP models often suffer from catastrophic forgetting.

[95] introduced a framework, dubbed LAMP that leverages pretrained LLMs for abductive reasoning to improve future event prediction. The framework operates through a multi-step process: first, the event sequence model proposes predictions on future events given historical data; second, a fine-tuned LLM learns to suggest possible causes for each prediction proposal through few-shot learning; third, a search module identifies previous events that match the suggested causes; finally, a scoring function evaluates whether the retrieved events could actually cause the proposed future events. This approach effectively combines the temporal modeling capabilities of TPPs with the sophisticated reasoning abilities of LLMs. They show improvements over state-of-the-art event sequence models across several datasets, highlighting that the reasoning capability of LLMs can enhance traditional TPP prediction tasks.

5.2 Direct LLM-TPP Integration↩︎

[94] developed TPP-LLM, a novel framework integrating LLMs with TPPs to capture both semantic and temporal aspects of event sequences. TPP-LLM directly utilizes textual descriptions to capture rich semantic information instead of categorical event information. It incorporates temporal embeddings through conventional positional encoding and employs parameter-efficient fine-tuning (PEFT) methods to learn temporal dynamics without extensive retraining, significantly improving predictive accuracy and computational efficiency. Experimental results on real-world datasets demonstrate superior performance in sequence modeling and event prediction tasks, establishing the benefits of combining semantic understanding with temporal modeling for enhanced TPP capabilities.

Similarly, [96] proposed Language-TPP, which also integrates TPPs with LLMs for enhanced event sequence modeling but introduces a fundamentally different temporal encoding approach. The key innovation is a novel temporal encoding mechanism that converts continuous time intervals into specialized byte-tokens, enabling seamless integration with standard LLM architectures without requiring extensive architectural modifications. This byte-token representation allows the model to process temporal information as naturally as textual data, maintaining the LLM’s inherent processing capabilities while incorporating temporal dynamics. The framework achieves state-of-the-art performance across multiple TPP tasks including event time prediction, type prediction, and intensity estimation on five diverse datasets. Additionally, they also explored a novel task of event description generation enabled by LLMs. The authors demonstrate that incorporating temporal information significantly improves the quality of generated event descriptions, showcasing the bidirectional benefits of the integration.

5.3 Other Extensions↩︎

[99] introduced TPP-Embedding for temporal event sequence retrieval from textual descriptions, targeting applications in e-commerce behavior analysis, social media monitoring, and criminal incident tracking. The authors developed TESRBench, a comprehensive benchmark with diverse real-world datasets and synthesized textual descriptions. Their model leverages the TPP-LLM [94] framework to integrate LLMs with TPPs, encoding both event texts and temporal information through pooling representations and contrastive loss to align sequence-level embeddings with textual descriptions. This approach outperforms baseline models across TESRBench datasets and establishes foundations for retrieval-augmented generation in TPP domains.

Extending beyond textual data, [98] introduced DanmakuTPPBench, a comprehensive multi-modal TPP benchmark addressing the gap in datasets with temporal, textual, and visual information. The benchmark comprises two components: DanmakuTPP-Events, derived from Bilibili’s user-generated bullet comments (Danmaku) that form multi-modal events with timestamps, textual content, and video frames; and DanmakuTPP-QA, a question-answering dataset constructed via a multi-agent pipeline using state-of-the-art LLMs and multi-modal LLMs (MLLMs). Targeting temporal-textual-visual reasoning tasks requiring multi-modal event dynamics understanding, extensive evaluations of classical TPP models and recent MLLMs reveal significant performance gaps in modeling multi-modal event sequences. This work establishes baselines and opens research directions for integrating TPP modeling into the multi-modal language modeling landscape. An illustration of DanmakuTPPBench is shown in 6.

Overall, this research direction is rapidly evolving, with approaches ranging from LLM-inspired techniques to direct integration and multi-modal extensions. The field shows promise for further exploration in areas such as retrieval-augmented TPPs, multi-modal event understanding, and more sophisticated temporal reasoning capabilities.

6 Model Training↩︎

In this section, we focus on frequentist methods for estimating model parameters in neural TPPs and LLM-based TPPs. The training objective is to align the model distribution \(f_\theta(\mathcal{T})\) with the empirical data distribution \(f(\mathcal{T})\), where \(\mathcal{T}\) denotes an observed event sequence. This is typically formulated as a divergence minimization problem: \[\hat{\theta} = \mathop{\mathrm{arg\,min}}_\theta \text{D}(f(\mathcal{T}) \,\|\, f_\theta(\mathcal{T})),\] where \(\text{D}(\cdot \,\|\, \cdot)\) is a statistical distance or divergence measure. Different choices of \(\text{D}\) correspond to different estimation principles. Below we summarize four representative approaches.

6.1 KL Divergence↩︎

KL divergence is the most commonly used training criterion. Minimizing the KL divergence \(\text{KL}(f(\mathcal{T}) \,\|\, f_\theta(\mathcal{T}))\) is equivalent to maximizing the log-likelihood: \[\hat{\theta} = \mathop{\mathrm{arg\,max}}_\theta \log f_\theta(\mathcal{T}).\] This approach necessitates explicit evaluation of the conditional intensity function and its integral over time, which is typically intractable and must be approximated via numerical integration. Despite its computational burden, MLE remains asymptotically efficient and statistically optimal. Due to its solid theoretical foundation and general applicability, MLE has been adopted in the vast majority of TPP studies [100], [101].

6.2 Wasserstein Distance↩︎

To circumvent the need for numerical integration inherent in MLE, [102] proposed leveraging the Wasserstein distance within a Wasserstein GAN (WGAN) framework for TPPs. In this setup, the TPP model acts as a generator, and a critic network learns to distinguish between real and generated event sequences by minimizing the Wasserstein distance. Compared to KL divergence, Wasserstein distance provides smoother gradients and better handles mode collapse. However, the adversarial nature of WGANs may lead to training instability and increased complexity. To further improve upon this approach, [103] introduced a likelihood-free training method based directly on the Wasserstein distance between point processes. Unlike the earlier WGAN-based framework, which primarily learns the aggregate intensity over a dataset, their method enables individual-level, in-sample forward prediction of event sequences conditioned on historical context. Moreover, Wasserstein distance captures the underlying geometric structure of event sequences more effectively than KL divergence, leading to more robust alignment between the generated and real sequences.

6.3 Noise Contrastive Estimation↩︎

Noise contrastive estimation (NCE) reframes parameter learning as a binary classification problem that distinguishes between observed data and artificially generated noise samples. This approach bypasses direct likelihood computation, making it particularly suitable for models with intractable likelihoods, such as TPPs. Both [104] and [105] applied NCE techniques to estimate the parameters of TPPs. The key difference between their approaches lies in the specific NCE variant they adopt. [104] employed the original binary classification formulation proposed by [106], which treats the learning task as discriminating between real and noise sequences. In contrast, [105] adopted a ranking-based NCE variant introduced by [107], which is better suited for modeling conditional distributions—a natural fit for TPPs where future events are conditioned on historical ones.

6.4 Fisher Divergence↩︎

Fisher divergence, also known as score matching [108], quantifies distributional differences based on score functions (i.e., gradients of log-densities). This method avoids the need to compute normalizing constants, making it appealing for models with intractable likelihoods. Several prior works have introduced score matching into the context of point processes. For example, [109] derived score matching estimators for classical Poisson processes. [110] extended this framework to deep covariate-based spatio-TPPs, while [111] further generalized it to multivariate Hawkes processes. These efforts have significantly advanced the applicability of score matching in point process modeling. However, recent work by [112] highlights a critical limitation: the estimators proposed in these earlier studies are incomplete and only valid for specific classes of point processes. In more general cases—including some simple parametric models—these methods fail to produce accurate parameter estimates. To address this issue, [112] proposed a weighted (autoregressive) score matching estimator that generalizes to a broader class of point process models, offering improved theoretical soundness and practical applicability.

6.5 Comparison between Different Estimators↩︎

Among these methods, only MLE requires computing the integral of the conditional intensity function. The others avoid this step, potentially improving training efficiency. All methods are consistent under mild conditions, but MLE remains asymptotically optimal in terms of variance. Alternative methods may exhibit slower convergence and higher asymptotic variance, but offer better scalability and computational simplicity.

7 Applications↩︎

TPPs have a wide range of applications, including in seismology, finance, neuroscience, social networks, and epidemiology. Broadly speaking, these applications can be categorized into two main types: event prediction and causal discovery. An overview of applications of TPPs is shown in 7.

7.1 Application in Event Prediction↩︎

Event prediction leverages historical data to forecast the timing, frequency, and types of future events, with applications spanning social networks, epidemiology, earthquake forecasting, finance, and recommendation systems.

In social networks, Hawkes processes and neural TPPs are widely employed to model temporal interactions and information diffusion. For instance, [113] introduced a neural TPP for detecting coordinated behavior, while [114], [115] applied Hawkes processes to retweet dynamics prediction, and [116] further analyze the higher-order social interactions. Recent work enhances interpretability [82] and predicting information popularity [117]. Additionally, these methods aid in anomaly detection, such as identifying fake accounts [118], and [119] integrated sentiment analysis to detect fake news trends during the COVID-19 pandemic.

In epidemiology, Hawkes processes effectively model disease propagation, as demonstrated by [120]. Studies such as [121] incorporated mobility data to predict COVID-19 transmission patterns, while [122] leveraged similar data for early outbreak forecasting. Further contributions include estimating transmission times [123] and assessing pandemic impacts using self-exciting processes [124].

In earthquake forecasting, spatio-temporal Hawkes processes are instrumental in capturing aftershock sequences. For instance, [125] introduced spatio-temporal TPP models for earthquake occurrences, with recent advancements enhancing flexibility [126] and refining decay rate modeling [127].

In financial markets, Hawkes processes analyze market microstructure and limit order book dynamics [128]. Neural Hawkes processes [129] and [130] further improve predictive accuracy in high frequency financial data.

Recommendation systems leverage Hawkes processes to model sequential user behavior. [131] combined them with attention mechanisms for sequential recommendations. These techniques can also be used to predict a user’s future shopping times and item types based on their past purchase history, enabling targeted promotional strategies [61], [82].

7.2 Application in Causal Discovery↩︎

In Hawkes process, causal discovery aim to recover the causal structure among difference event types from observed event sequence data. Applications in this category are prevalent in areas like neuroscience, finance, AI for operations, social network, healthcare, and cybersecurity. Here, the focus is not on predicting future events but on uncovering dependencies between event types, often referred to as Granger causality [132]. These causal relationships enable better decision-making and provide mechanistic understanding of complex event dynamics. For example, in neuroscience, each neuron can be considered as an event type, and its spike train forms a univariate point process. The spike trains of multiple neurons naturally constitute a multivariate point process. The goal is to determine whether there exists functional connectivity between neurons [2], [53]. Similarly, in high-frequency financial trading, a large number of asks (sell orders) and bids (buy orders) occur within short periods. Here, all sell orders are treated as one event type, and all buy orders as another. The primary interest lies in understanding the mutual influence between buy and sell orders in the order book [3]. In AI operations, TPPs identify system failure root causes by distinguishing primary triggers from secondary effects through their causal structure, guiding prioritized fixes [133]. Social network analysis similarly leverages Hawkes processes to quantify mutual influence patterns, where user actions (posts, likes, shares) as distinct event types reveal how influential users trigger reaction cascades [134]. Healthcare applications employ TPPs to analyze drug reactions and symptom interactions for improved treatment strategies [135]. Cybersecurity implementations further demonstrate their value in attack pattern analysis for enhanced defense mechanisms [136], [137].

The modeling framework for these applications in causal discovery builds upon multivariate Hawkes processes (MHP) where a \(K\)-varaite MHP can be formulated as a collection of \(K\) univariate TPPs with the conditional intensity function taking the form of 5 . Crucially, we say process \(k^\prime\) does not Granger-cause process \(k\) if and only if \(\phi_{k,k^\prime}(\cdot) = 0\). Therefore, estimating these triggering functions thus directly translates to learning the causal structure.

While MLE offers a straightforward solution for learning causal structure, it often produces spurious connections due to finite samples and the lack of sparsity constraints. To tackle this issue, various types of sparsity approaches have been developed, typically categorized into constraint-based and score-based methodologies.

Constraint-based approaches address the problem through statistical testing to prune spurious edges. For instance, [138] proposed a general constraint-based framework for learning causal structure in time series data using conditional independence tests, but it is only applicable to discrete-time processes. Later, [139] proposed a screening algorithm for the Hawkes process, extending it to the continuous-time domain.

In contrast, score-based approaches learn causal structure by optimizing a well-defined criterion with various sparsity regularizations to enforce structural sparsity. For instance, [140] proposed a nonparametric Hawkes process model with group sparsity regularization, while [134] employed both nuclear norm and \(\ell_1\) norm as sparsity regularizations. [141] used \(\ell_0\)-regularization via an \(\epsilon\)-sparsity approach [142]. Alternatively, based on data compression techniques, [143] proposed a minimum description length (MDL) criterion for Hawkes processes. Using a similar approach, [144] introduced a minimum message length (MML) criterion, which extends the MDL framework. This extension incorporates prior distributions—such as expert knowledge—over model parameters, enhancing flexibility in structure-related penalization.

By further assuming process stability, [145] demonstrated that the inference procedure can be accelerated using a cumulant matching strategy based on the analytical form of cumulants [146], thereby eliminating the need to estimate the triggering function.

Recently, several deep point process-based methods have been proposed. For example, [147] introduced an attribution method to uncover Granger causality, and [148] explored instance-wise causal structures using Transformer self-attention, aligning the mechanism with Granger causality principles.

Another line of research addresses more realistic scenarios where processes may exhibit nonstationarity, topological dependence, and insufficient temporal resolution. First, nonstationary dynamics frequently emerge when event interactions and background intensities vary over time [149]. [150] further demonstrated that the causal mechanisms can change over time. Second, network effects often exist, where events are influenced not only by their own history but also by their topological neighbors. Failure to account for these dependencies can lead to biased causal estimates. To handle topological dependencies, [133] proposed the Topological Hawkes Process (THP), which extends temporal convolution to graph-time convolution while employing an EM-based inference approach. Subsequent improvements by [151] enhanced THP’s scalability through gradient-based optimization, and [152] further generalized the framework using causal-attention Transformers to capture complex network relationships.

For low-resolution scenarios, [153] demonstrated the robustness of cumulant-based methods to observational noise. Later, [154] incorporated delayed effects using a delay-aware MDL criterion to handle observation delays. Furthermore, [155] showed that low resolution may lead to instantaneous effects and thus developed a discrete-time structural Hawkes process to handle such instantaneous causal relationships.

8 Challenges↩︎

Although TPPs have achieved significant advancements in recent years, several critical challenges remain. While [9] provided an overview of experimental issues, we focus here on open problems related to model design and practical applications.

8.0.0.1 Data and Model

A key challenge is the lack of a unified experimental setup across existing neural TPP studies, which hampers reproducibility and limits fair performance comparisons. To address this, [156] conducted a large-scale empirical study evaluating various architectures under a consistent setup across diverse synthetic and real-world datasets. Their results provide valuable insights into the effects of event encoding, history encoding, and decoder parameterization. Nevertheless, a standard benchmark protocol for model development and comparison is still lacking.

8.0.0.2 Model Interpretability

Traditional statistical TPPs (e.g., Poisson, Hawkes) offer interpretable components with clear physical meanings, such as base intensity and excitation kernels. In contrast, neural TPPs, while more expressive, typically operate as black-box models with limited transparency. Although recent studies have proposed structured designs or post-hoc explanations to improve interpretability, this remains an underexplored area. Designing inherently interpretable neural architectures is a promising direction for future research.

8.0.0.3 Model Scalability

Neural TPPs often adopt Transformer-based architectures to capture long-range dependencies. However, standard attention mechanisms incur quadratic time and space complexity \(\mathcal{O}(N^2)\), which severely limits their scalability on long sequences. Although recent architectures such as state space models and Mamba offer linear-complexity alternatives, their adaptation to TPP modeling is still in its early stages and requires further theoretical and empirical investigation.

8.0.0.4 Sampling Efficiency

Most existing works emphasize expressiveness and training efficiency, while sampling efficiency has been largely overlooked. Current sampling procedures, such as thinning for intensity-based models and inverse transform sampling for density-based models, are inherently sequential, limiting real-time and parallel applications. Recent advances in speculative decoding for autoregressive models could inspire new directions for parallel sampling in TPPs.

8.0.0.5 Multimodal Modeling

Conventional TPPs typically model simple event features such as timestamp and categorical mark. However, many real-world applications involve multimodal information (e.g., text, image, audio) associated with events. Inspired by the success of multimodal models in vision and language, recent research has begun to explore integrating rich modalities into TPPs. In particular, LLM-based TPPs offer a promising paradigm for modeling and generating multimodal event sequences, though this area remains in its infancy and requires further development in representation alignment and efficient fusion mechanisms.

9 Conclusions↩︎

TPPs provide a powerful mathematical framework for modeling asynchronous event sequences across diverse domains such as neuroscience, finance, and social media. Over the years, the field has evolved from traditional parametric models to increasingly flexible nonparametric and neural approaches. In this survey, we reviewed recent progress in three major directions of TPP research: Bayesian TPPs, neural TPPs, and LLM-based TPPs. Each paradigm offers distinct modeling philosophies and strengths—Bayesian methods emphasize uncertainty quantification and principled inference, neural methods prioritize expressive power and scalability, while LLM-based approaches open new avenues for handling complex multimodal sequences. We also provided a taxonomy of representative works, highlighted core modeling principles and estimation strategies, and discussed ongoing challenges including model interpretability, scalability, and sampling efficiency. In particular, we emphasized recent developments since 2020 that have not been fully covered in earlier surveys, including Bayesian nonparametric models and the emerging trend of applying large language models to TPPs. Looking ahead, we believe that continued progress will hinge on deeper integration across statistical rigor, neural flexibility, and language-model capabilities. Addressing key challenges and enabling broader application of TPPs in real-world, multimodal, and high-resolution settings remain central to advancing the field.

Acknowledgments↩︎

This work was supported by NSFC Project (No. 62106121), the MOE Project of Key Research Institute of Humanities and Social Sciences (22JJD110001), the fundamental research funds for the central universities, and the research funds of Renmin University of China (24XNKJ13).

References↩︎

---

1. Feng Zhou and Cheng Wan are with Center for Applied Statistics and School of Statistics, Renmin University of China, Beijing, China. Email: feng.zhou@ruc.edu.cn, wancheng0256@ruc.edu.cn;↩︎

2. Quyu Kong is with Alibaba Cloud, Hangzhou, China. Email: kongquyu@gmail.com;↩︎

3. Jie Qiao and Ruichu Cai are with School of Computer Science, Guangdong University of Technology, Guangzhou, China. Email: qiaojie.chn@gmail.com, cairuichu@gmail.com;↩︎

4. Yixuan Zhang is with School of Statistics and Data Science, Southeast University, Nanjing, China. Email: zh1xuan@hotmail.com.↩︎

5. The conditional intensity function must satisfy certain conditions.↩︎

6. It requires \(\int_0^\infty \phi(t) \, dt < 1\) to ensure the Hawkes process does not explode.↩︎

7. Some studies treat the baseline intensity as a constant, but here we consider it as a more general function.↩︎