Neural Temporal Point Processes: A Review

1 Introduction

Many applications in science and industry are concerned with collections of events with timestamps. Earthquake occurrences in seismology, neural spike trains in neuroscience, trades and orders in a financial market, and user activity logs on the web, can all be represented as sequences of discrete (instantaneous) events observed in continuous time.

Temporal point processes (TPP) are probabilistic models for such event data [Daley and Vere-Jones, 2007]. More specifically, TPPs are generative models of variable-length point sequences observed on the real half-line-here interpreted as arrival times of events. TPPs are built on rich theoretical foundations, with early work dating back to the beginning of the 20th century, where they were used to model the arrival of insurance claims and telephone traffic [Brockmeyer et al., 1948; Cram´er, 1969]. The field underwent rapid development in the second half of the century, and TPPs were applied to a wide array of domains including seismology, neuroscience, and finance.

Nevertheless, TPPs entered the mainstream of machine learning research only very recently. One of the exciting ideas developed at the intersection of the fields of point processes and machine learning were neural TPPs [Du et al., 2016; Mei and Eisner, 2017]. Classical (non-neural) TPPs can only capture relatively simple patterns in event occurrences, such as self-excitation [Hawkes, 1971]. In contrast, neural TPPs are able to learn complex dependencies, and are often even computationally more efficient than their classical counterparts. As such, the literature on neural TPPs has witnessed rapid growth since their introduction.

Scope and structure of the paper. The goal of this survey is to provide an overview of neural TPPs, with focus on models (Sections 3-5) and their applications (Section 6). Due to limited space, we do not attempt to describe every existing approach in full detail, but rather focus on general principles and building blocks for constructing neural TPP models.

We also discuss the main challenges that the field currently faces and outline some future research directions (Section 7). For other reviews of TPPs for machine learning, we refer the reader to the tutorial by [Gomez-Rodriguez and Valera, 2018]; and two recent surveys by [Yan, 2019], who also covers non-neural approaches, and [Enguehard et al., 2020] who experimentally compare neural TPP architectures in applications to healthcare data. Our work provides a more detailed overview of neural TPP architectures and their applications compared to the above papers.

2 Background and Notation

A TPP [Daley and Vere-Jones, 2007] is a probability distribution over variable-length sequences in some time interval [0, T]. A realization of a marked TPP can be represented as an event sequence X = {(t1,m1), . . . , (tN,mN)}, where N, the number of events, is itself a random variable. Here, 0 < t1 < · · · < tN ≤ T are the arrival times of events and mi ∈ M are the marks. Categorical marks (i.e., M = {1, . . . ,K}) are most commonly considered in practice, but other choices, such as M = RD , are also possible. Sometimes, it is convenient to instead work with the interevent times τi = ti - ti-1, where t0 = 0 and tN+1 = T . For a given X, we denote the history of past events at time t as Ht = {(tj,mj) : tj < t}.

A distribution of a TPP with K categorical marks can be characterized by K conditional intensity functions λ∗k (t) (one for each mark k) that are defined as

where the ∗ symbol is used as a shorthand for conditioning on the history Ht. Note that the above definition also applies to unmarked TPPs if we set K = 1. While the conditional intensity is often mentioned in the literature, it is not the only way to characterize a TPP-we will consider some alternatives in the next section.

3 Autoregressive Neural TPPs

Neural TPPs can be defined as autoregressive models, as done in the seminal work by [Du et al., 2016]. Such autoregressive TPPs operate by sequentially predicting the time and mark of the next event (ti,mi). Usually, we can decompose this procedure into 3 steps (see Figure 1):

1. Represent each event (tj,mj) as a feature vector yj.
2. Encode the history Hti (represented by a sequence of feature vectors (y1, . . . , yi-1)) into a fixed-dimensional history embedding hi.
3. Use the history embedding hi to parametrize the conditional distribution over the next event Pi(ti,mi|Hti ).

We will now discuss each of these steps in more detail.

3.1 Representing Events as Feature Vectors

First, we need to represent each event (tj,mj) as a feature vector yj that can then be fed into an encoder neural network (Section 3.2). We consider the features ytime j based on the arrival time tj (or inter-event time τj) and ymark j based on the mark mj. The vector yj is obtained by combining ytime j and ymark j , e.g., via concatenation.

Time features ytime j . Earlier works used the inter-event time τj or its logarithm log τj as the time-related feature [Du et al., 2016; Omi et al., 2019]. Recently, [Zuo et al., 2020] and [Zhang et al., 2020a] proposed to instead obtain features from tj using trigonometric functions, which is based on positional encodings used in transformer language models [Vaswani et al., 2017].

Mark features ymark j . Categorical marks are usually encoded with an embedding layer [Du et al., 2016]. Real-valued marks can be directly used as ymark j . This ability to naturally handle different mark types is one of the attractive properties of neural TPPs (compared to classical TPP models).

3.2 Encoding the History into a Vector

The core idea of autoregressive neural TPP models is that event history Hti (a variable-sized set) can be represented as a fixed-dimensional vector hi [Du et al., 2016]. We review the two main families of approaches for encoding the past events (y1, . . . , yi-1) into a history embedding hi next.

Recurrent encoders start with an initial hidden state h1. Then, after each event (ti,mi) they update the hidden state as hi+1 = Update(hi, yi). The hidden states hi are then used as the history embedding. The Update function is usually implemented based on the RNN, GRU or LSTM update equations [Du et al., 2016; Xiao et al., 2017b].

The main advantage of recurrent models is that they allow us to compute the history embedding hi for all N events in the sequence in O(N) time. This compares favorably even to classical non-neural TPPs, such as the Hawkes process, where the likelihood computation in general scales as O(N 2 ). One downside of recurrent models is their inherently sequential nature. Because of this, such models are usually trained via truncated backpropagation through time, which only provides an approximation to the true gradients [Sutskever, 2013].

Set aggregation encoders directly encode the feature vectors (y1, . . . , yi-1) into a history embedding hi. Unlike recurrent models, here the encoding is done independently for each i. The encoding operation can be defined, e.g., using self-attention [Zuo et al., 2020; Zhang et al., 2020a]. It is postulated that such encoders are better at capturing longrange dependencies between events compared to recurrent encoders. However, more thorough evaluation is needed to validate this claim (see Section 7). On the one hand, set aggregation encoders can compute hi for each event in parallel, unlike recurrent models. On the other hand, usually the time of this computation scales as O(N 2 ) with the sequence length N, since each hi depends on all the past events (and the model does not have a Markov property). This problem can be mitigated by restricting the encoder to only the last L events (yi-L, . . . , yi-1), thus reducing the time complexity to O(NL).

3.3 Predicting the Time of the Next Event

For simplicity, we start by considering the unmarked case and postpone the discussion of marked TPPs until the next section. An autoregressive TPP models the distribution of the next arrival time ti given the history Hti . This is equivalent to considering the distribution of the next inter-event time τi given Hti , which we denote as P ∗ i (τi). 1 The distribution P ∗ i (τi) can be represented by any of the following functions:

1. probability density function f ∗ i (τi)
2. cumulative distribution function F ∗ i (τi) = ∫ τi 0 f ∗ i (u)du
3. survival function S∗ i (τi) = 1 -F ∗ i (τi)
4. hazard function φ∗i (τi) = f ∗ i (τi)/S∗ i (τi)
5. cumulative hazard function Φ∗i (τi) = ∫ τi 0 φ∗i (u)du .

In an autoregressive neural TPP, we pick a parametric form for one of the above functions and compute its parameters using the history embedding hi. For example, the conditional PDF f ∗ i might be obtained as

Here f(·|θ) is some parametric density function over [0,∞) (e.g., exponential density) and W, b are learnable parameters. A nonlinear function σ(·) can be used to enforce necessary constraints on the parameters, such as non-negativity.

It is important to ensure that the chosen parametrization defines a valid probability distribution. For instance, the PDF f ∗ i must be non-negative and satisfy ∫ ∞ 0 f ∗ i (u)du = 1. This corresponds to the cumulative hazard function Φ∗i being strictly increasing, differentiable and satisfying Φ∗i (0) = 0 and limτ→∞ Φ∗i (τ ) = ∞.

Some parametrizations of P ∗ i (τ) proposed in the literature fail to satisfy the above conditions. For example, the hazard function φ∗i (τ) = exp(wτ + b) [Du et al., 2016] fails to satisfy limτ→∞Φ∗i (τ) = ∞ if the parameter w is negative. One more example is the cumulative hazard function defined by a single-hidden-layer neural network with positive weights w, v ∈ RD+ as Φ∗i (τ) = softplus (∑Dd=1 vd tanh(wdτ + bd)) [Omi et al., 2019] that fails to satisfy both Φ∗i (0) = 0 and limτ→∞Φ∗i (τ) = ∞. Since above parametric functions do not define a valid distribution over the inter-event times, sampling methods have a non-zero probability of failing or producing invalid event sequences. Therefore, it's crucial to pick a valid parametrization when designing a neural TPP model.

Specifying one of the functions (1) - (5) listed above uniquely identifies the conditional distribution P ∗ i (τi), and thus the other four functions in the list. This, however, does not mean that choosing which function to parametrize is unimportant. In particular, some choices of the hazard function φ∗i cannot be integrated analytically, which becomes a problem when computing the log-likelihood (as we will see in Section 5). In contrast, it is usually trivial to obtain φ∗i from any parametrization of Φ∗i , since differentiation is easier than integration [Omi et al., 2019]. More generally, there are three important aspects that one has to keep in mind when specifying the distribution P ∗ i (τi):

1 We again use ∗ to denote conditioning on the history Hti .

• Flexibility: Does the given parametrization of P ∗ i (τi) allow us to approximate any distribution, e.g., a multimodal one?
• Closed-form likelihood: Can we compute either the CDF F ∗ i , SF S∗ i or CHF Φ∗i analytically? These functions are involved in the log-likelihood computation (Section 5), and therefore should be computed in closed form for efficient model training. Approximating these functions with Monte Carlo or numerical quadrature is slower and less accurate.
• Closed-form sampling: Can we draw samples from P ∗ i (τi) analytically? In the best case, this should be done with inversion sampling [Rasmussen, 2011], which requires analytically inverting either F ∗ i , S∗ i or Φ∗i . Inverting these functions via numerical root-finding is again slower and less accurate. Approaches based on thinning [Ogata, 1981] are also not ideal, since they do not benefit from parallel hardware like GPUs. Moreover, closedform inversion sampling enables the reparametrization trick [Mohamed et al., 2020], which allows us to train TPPs using sampling-based losses (Section 5.2).

Existing approaches offer different trade-offs between the above criteria. For example, a simple unimodal distribution offers closed-form sampling and likelihood computation, but lacks flexibility [Du et al., 2016]. One can construct more expressive distributions by parametrizing the cumulative hazard Φ∗i either with a mixture of kernels [Okawa et al., 2019; Zhang et al., 2020b] or a neural network [Omi et al., 2019], but this will prevent closed-form sampling. Specifying the PDF f ∗ i with a mixture distribution [Shchur et al., 2020a] or Φ∗i with invertible splines [Shchur et al., 2020b] allows to define a flexible model where both sampling and likelihood computation can be done analytically. Finally, parametrizations that require approximating Φ∗i via Monte Carlo integration are less efficient and accurate than all of the abovementioned approaches [Omi et al., 2019].

As a side note, a more flexible parametrization might also be more difficult to train or more prone to overfitting. Therefore, the choice of the parametrization is an important modeling decision that depends on the application.

Lastly, we would like to point out that the view of a TPP as an autoregressive model naturally connects to the traditional conditional intensity characterization (Equation 1). The conditional intensity λ∗ (t) can be defined by stitching together the hazard functions φ∗i

In the TPP literature, the hazard function φ∗i is often called "intensity," even though the two are, technically, different mathematical objects.

3.4 Modeling the Marks

In a marked autoregressive TPP, one has to parametrize the conditional distribution P ∗ i (τi,mi) using the history embedding hi. We first consider categorical marks, as they are most often used in practice.

Conditionally independent models factorize the distribution P ∗ i (τi,mi) into a product of two independent distributions P ∗ i (τi) and P ∗ i (mi) that are both parametrized using hi [Du et al., 2016; Xiao et al., 2017b]. The time distribution P ∗ i (τi) can be parametrized using any of the choices described in Section 3.3, such as the hazard function φ∗i (τ). The mark distribution P ∗ i (mi) is a categorical distribution with probability mass function p∗i (mi = k). In this case, the conditional intensity for mark k is computed as

where (i-1) is the index of the most recent event before time t. Note that if we set K = 1, we recover the definition of the intensity function from Equation 3.

While conditionally independent models require fewer parameters to specify P ∗ i (τi,mi), recent works suggest that this simplifying assumption may hurt predictive performance [Enguehard et al., 2020]. There are two ways to model dependencies between τi and mi that we consider below.

Time conditioned on marks [Zuo et al., 2020; Enguehard et al., 2020]. In this case, we must specify a separate distribution P ∗ i (τi|mi = k) for each mark k ∈ {1, . . . ,K}. Suppose that for each k we represent P ∗ i (τi|mi = k) with a hazard function φ∗ik (τ). Then the conditional intensity λ∗k (t) for mark k is computed simply as

It is possible to model the dependencies across marks on a coarser grid in time, which significantly improves the scalability in the number of marks [T¨urkmen et al., 2019b].

Marks conditioned on time. Here, the inter-event time is distributed according to P ∗ i (τi), and for each τ we need to specify a distribution P ∗ (mi|τi = τ). We again assume that the time distribution P ∗ i (τi) is described by a hazard function φ∗i (τ), and P ∗ i (mi|τi = τ) can be parametrized, e.g., using a Gaussian process [Biloˇs et al., 2019]. In this case the conditional intensity λ∗k (t) is computed as

where we used notation analogous to Equation 4. The term φ∗i (t - ti-1) is often referred to as "ground intensity."

Other mark types. A conditionally independent model can easily handle any type of marks by specifying an appropriate distribution P ∗ i (mi). Dependencies between continuous marks and the inter-event time can be incorporated by modeling the joint density f ∗ i (τi,mi) [Zhu et al., 2020].

4 Continuous-time State Evolution

Another line of research has studied neural TPPs that operate completely in continuous time. Such models define a leftcontinuous state h(t) at all times t ∈ [0, T]. The state is initialized to some value h(0). Then, for each event i the

The Evolve(h(ti-1), ti-1, ti) procedure evolves the the state continuously over the time interval (ti-1, ti] between the events. Such state evolution can be implemented either via exponential decay [Mei and Eisner, 2017] or, more generally, be governed by an ordinary differential equation [Rubanova et al., 2019; Jia and Benson, 2019]. The Update(h(ti), yi) operation performs an instantaneous update to the state, similarly to the recurrent encoder from last section.

While the above procedure might seem similar to a recurrent model from Section 3.2, the continuous-time model uses the state h(t) differently from the autoregressive model. In a continuous-time model, the state h(t) is used to directly define the intensity λ∗k for each mark k as

where gk : RH → R>0 is a non-linear function that maps the hidden state h(t) to the value of the conditional intensity for mark k at time t. Such function, for example, can be implemented as gk(h(t)) = softplus(wT k h(t)) or a multilayer perceptron [Chen et al., 2021b].

To summarize, in a continuous-time state model, the state h(t) is defined for all t ∈ [0, T], and the intensity λ∗k (t) at time t depends only on the current state h(t). In contrast, in an autogressive model, the discrete-time state hi is updated only after an event occurs.. Hence, the state hi defines the entire conditional distribution P ∗ i (τi,mi), and therefore the intensity λ∗k (t) in the interval (ti, ti+1].

Discussion. Continuous-time models have several advantages compared to autoregressive ones. For example, they provide a natural framework for modeling irregularlysampled time series, which is valuable in medical applications [Enguehard et al., 2020]. By modeling the unobserved attributes as a function of the state h(t), it is easy to estimate each attribute at any time t [Rubanova et al., 2019]. These models are also well-suited for modeling spatio-temporal point processes (i.e., with marks in M = RD ) [Jia and Benson, 2019; Chen et al., 2021b].

However, this flexibility comes at a cost: evaluating both the state evolution (Equation 7) and the model likelihood (Equation 9) requires numerically approximating intractable integrals. This makes training in continuous-time models slower than for autoregressive ones. Sampling similarly requires numerical approximations.

5 Parameter Estimation

5.1 MaximumLikelihood Estimation

Negative log-likelihood (NLL) is the default training objective for both neural and classical TPP models. NLL for a single sequence X with categorical marks is computed as

The log-likelihood can be understood using the following two facts. First, the quantity λ∗k (ti)dt corresponds to the probability of observing an event of type k in the infinitesimal interval [ti, ti + dt) conditioned on the past events Hti . Second, we can compute the probability of not observing any events of type k in the rest of the interval [0, T] as exp (- ∫ T 0 λ∗k (u)du). By taking the logarithm, summing these expressions for all events (ti,mi) and event types k, and finally negating, we obtain Equation 9.

Computing the NLL for TPPs can sometimes be challenging due to the presence of the integral in the second line of Equation 9. One possible solution is to approximate this integral using Monte Carlo integration [Mei and Eisner, 2017] or numerical quadrature [Rubanova et al., 2019; Zuo et al., 2020]. However, some autoregressive neural TPP models allow us to compute the NLL analytically, which is more accurate and computationally efficient. We demonstrate this using the following example.

Suppose we model P ∗ i (τi,mi) using the "time conditioned on marks" approach from Section 3.4. That is, we specify the distribution P ∗ i (τi|mi = k) for each mark k with a cumulative hazard function Φ∗ik (τ). By combining Equation 5 and Equation 9, we can rewrite the expression for the NLL as

Assuming that our parametrization allows us to compute Φ∗ik (τ) analytically, we are now able to compute the NLL in closed form (without numerical integration). Remember that the hazard function φ∗ik can be easily obtained by differentiation as φ∗ik (τ) = ∂ ∂τ Φ∗ik (τ). Finally, note that the NLL can also be expressed in terms of, e.g., the conditional PDFs f ∗ ik or survival functions S∗ ik (Section 3.3).

Evaluating the NLL in Equation 10 can be still computationally expensive when K, the number of marks, is extremely large. Several works propose approximations based on noise-contrastive-estimation that can be used in this regime [Guo et al., 2018; Mei et al., 2020].

Training. Usually, we are given a training set Dtrain of sequences that are assumed to be sampled i.i.d. from some unknown data-generating process. The TPP parameters (e.g., weights of the encoder in Section 3) are learned by minimizing the NLL of the sequences in Dtrain. This is typically done with some variant of (stochastic) gradient descent. In practice, the NLL loss is often normalized per sequence, e.g., by the interval length T [Enguehard et al., 2020]. Importantly, this normalization constant cannot depend on X, so it would be incorrect to, for example, normalize the NLL by the number of events N in each sequence. Finally, the i.i.d. assumption is not appropriate for all TPP datasets; [Boyd et al., 2020] show how to overcome this challenge by learning sequence embeddings.

5.2 Alternatives toMLE

TPPs can be trained using objective functions other than the NLL. Often, these objectives can be expressed as

Such sampling-based losses have been used by several approaches for learning generative models from the training sequences. These approaches aim to maximize the similarity between the training sequences in Dtrain and sequences X generated by the TPP model p(X) using a scoring function f(X). Examples include procedures based on Wasserstein distance [Xiao et al., 2017a], adversarial losses [Yan et al., 2018; Wu et al., 2018] and inverse reinforcement learning [Li et al., 2018].

Sometimes, the objective function of the form (11) arises naturally based on the application. For instance, in reinforcement learning, a TPP p(X) defines a stochastic policy and f(X) is the reward function [Upadhyay et al., 2018]. When learning with missing data, the missing events X are sampled from the TPP p(X), and f(X) corresponds to the NLL of the observed events [Gupta et al., 2021]. Finally, in variational inference, the TPP p(X) defines an approximate posterior and f(X) is the evidence lower bound (ELBO) [Shchur et al., 2020b; Chen et al., 2021a].

In practice, the gradients of the loss (Equation 11) w.r.t. the model parameters usually cannot be computed analytically and therefore are estimated with Monte Carlo. Earlier works used the score function estimator [Upadhyay et al., 2018], but modern approaches rely on the more accurate pathwise gradient estimator (also known as the "reparametrization trick") [Mohamed et al., 2020]. The latter relies on our ability to sample with reparametrization from P ∗ i (τi,mi), which again highlights the importance of the chosen parametrization for the conditional distribution, as described in Section 3.3.

On a related note, sampling-based losses for TPPs (Equation 11) can be non-differentiable, since N, the number of events in a sequence, is a discrete random variable. This problem can be solved by deriving a differentiable relaxation to the loss [Shchur et al., 2020b].

6 Applications

The literature on neural TPPs mostly considers their applications in web-related domains, e.g., for modeling user activity on social media. Most existing applications of neural TPPs fall into one of two categories:

• Prediction tasks, where the goal is to predict the time and / or type of future events;
• Structure discovery, where the tasks is to learn dependencies between different event types.

We now discuss these in more detail.

6.1 Prediction

Prediction is among the key tasks associated with temporal models. In case of a TPP, the goal is usually to predict the times and marks of future events given the history of past events. Such queries can be answered using the conditional distribution P ∗ i (τi,mi) defined by the neural TPP

model. Nearly all papers mentioned in previous sections feature numerical experiments on such prediction tasks. Some works combine elements of neural TPP models (Section 3) with other neural network architectures to solve specific realworld prediction tasks related to event data. We give some examples below.

Recommender systems are a recent application area for TPPs. Here, the goal is to predict the next most likely purchase event, in terms of both the time of purchase and the type (i.e., item), given a sequence of customer interactions or purchases in the past. Neural TPPs are especially well-suited for this task, as they can learn embeddings for large item sets using neural networks, similar to other neural recommendation models. Moreover, representing the temporal dimension of purchase behavior enables time-sensitive recommendations, e.g., used to time promotions. For example, [Kumar et al., 2019] address the next item prediction problem by embedding the event history to a vector, but are concerned only with predicting the next item type (mark). This approach is extended to a full model of times and events, with a hierarchical RNN model for intra- and inter-session activity represented on different levels, in [Vassøy et al., 2019].

Another common application of event sequence prediction is within the human mobility domain. Here, events are spatiotemporal, featuring coordinates in both time and space, potentially along with other marks. Examples include predicting taxi trajectories in a city (e.g., the ECML/PKDD 2015 challenge 2 ), or user check-ins on location-based social media. [Yang et al., 2018] address this task with a full neural TPP, on four different mobility data sets. This extends the approach in DeepMove [Feng et al., 2018], where the authors similarly use an RNNs to compute embeddings of timestamped sequences, but limit the predictions to the next location alone.

Other applications include clinical event prediction [Enguehard et al., 2020], predicting timestamped sequences of interactions of patients with the health system; human activity prediction for assisted living [Shen et al., 2018], and demand forecasting in sparse time series [T¨urkmen et al., 2019a].

6.2 Structure Discovery &Modeling Networks

In prediction tasks we are interested in the conditional distributions learned by the model. In contrast, in structure discovery tasks the parameters learned by the model are of interest.

For example, in latent network discovery applications we observe event activity generated by K users, each represented by a categorical mark. The goal is to infer an influence matrix A ∈ RK×K that encodes the dependencies between different marks [Linderman and Adams, 2014]. Here, the entries of A can be interpreted as edge weights in the network. Historically, this task has been addressed using non-neural models such as the Hawkes process [Hawkes, 1971]. The main advantage of network discovery approaches based on neural TPPs [Zhang et al., 2021] is their ability to handle more general interaction types.

Learning Granger causality is another task that is closely related to the network discovery problem. [Eichler et al.,

2 http://www.geolink.pt/ecmlpkdd2015-challenge/

2017] and [Achab et al., 2017] have shown that the influence matrix A of a Hawkes process completely captures the notion of Granger causality in multivariate event sequences. Recently, [Zhang et al., 2020b] generalized this understanding to neural TPP models, where they used the method of integrated gradients to estimate dependencies between event types, with applications to information diffusion on the web.

Neural TPPs have also been used to model information diffusion and network evolution in social networks [Trivedi et al., 2019]. The DyRep approach by Trivedi et al. generalizes an earlier non-neural framework, COEVOLVE, by [Farajtabar et al., 2017]. Similarly, Know-Evolve [Trivedi et al., 2017] models dynamically evolving knowledge graphs with a neural TPP. A related method, DeepHawkes, was developed specifically for modeling item popularity in information cascades [Cao et al., 2017].

Other applications. Neural TPPs have been featured in works in other research fields and application domains. For instance, [Huang et al., 2019] proposed an RNN-based Poisson process model for speech recognition. [Sharma et al., 2018] developed a latent-variable neural TPP method for modeling the behavior of larval zebrafish. Performing inference in the model allowed the authors to detect distinct behavioral patterns in zebrafish activity. Lastly, [Upadhyay et al., 2018] showed how to automatically choose timing for interventions in an interactive environment by combining a TPP model with the framework of reinforcement learning.

7 Open Challenges

We conclude with a discussion of what, in our opinion, are the main challenges that the field of neural TPPs currently faces.

7.1 Experimental Setup

Lack of standardized experimental setups and high-quality benchmark datasets makes a fair comparison of different neural TPP architectures problematic.

Each neural TPP model consists of multiple components, such as the history encoder and the parametrization of the conditional distributions (Sections 3 and 4). New architectures often change all these components at once, which makes it hard to pinpoint the source of empirical gains. Carefullydesigned ablation studies are necessary to identify the important design choices and guide the search for better models.

On a related note, the choice of baselines varies greatly across papers. For example, papers proposing autoregressive models (Section 3) rarely compare to continuous-time state models (Section 4), and vice versa. Considering a wider range of baselines is necessary to fairly assess the strengths and weaknesses of different families of approaches.

Finally, it is not clear whether the datasets commonly used in TPP literature actually allow us to find models that will perform better on real-world prediction tasks. In particular, Enguehard et al. point out that two popular datasets (MIMICII and StackOverflow) can be "solved" by a simple historyindependent baseline. Also, common implicit assumptions, such as treating the training sequences as i.i.d. (Section 5), might not be appropriate for existing datasets.

To conclude, open-sourcing libraries with reference implementations of various baseline methods and collecting large high-quality benchmark datasets are both critical next steps for neural TPP researchers. A recently released library with implementations of some autoregressive models by [Enguehard et al., 2020] takes a step in this direction, but, for instance, doesn't include continuous-time state models.

7.2 EvaluationMetrics

Many of the metrics commonly used to evaluate TPP models are not well suited for quantifying their predictive performance and have subtle failure modes.

For instance, consider the NLL score (Section 5), one of the most popular metrics for quantifying predictive performance in TPPs. As mentioned in Section 3.4, the NLL consists of a continuous component for the time ti and a discrete component for the mark mi. Therefore, reporting a single NLL score obscures information regarding the model's performance on predicting marks and times separately. For example, the NLL score is affected disproportionately by errors in marks as the number of marks increase. Moreover, the NLL can be "fooled" on datasets where the arrival times ti are measured in discrete units, such as seconds-a flexible model can produce an arbitrarily low NLL by placing narrow density "spikes" on these discrete values [Uria et al., 2013].

More importantly, the NLL is mostly irrelevant as a measure of error in real-world applications-it yields very little insight into model performance from a domain expert's viewpoint. However, other metrics, such as accuracy for the mark prediction, and mean absolute error (MAE) or mean squared error (MSE) for inter-event time prediction are even less suited for evaluating neural TPPs. These metrics measure the quality of a single future event prediction, and MAE/MSE only take a point estimate of τi into account. One doesn't need to model the entire distribution over τi (as done by a TPP) to perform well w.r.t. such metrics. If only single-event predictions are of interest, one could instead use a simple baseline that only models P ∗ i (mi) or produces a point estimate τ pred i . This baseline can be trained by directly minimizing absolute or squared error for inter-event times or crossentropy for marks. TPPs are, in contrast, probabilistic models trained with the NLL loss; so comparing them to pointestimate baselines using above metrics is unfair.

The main advantage of neural TPPs compared to simple "point-estimate" methods is their ability to sample entire trajectories of future events. Such probabilistic forecasts capture uncertainty in predictions and are able to answer more complex prediction queries (e.g., "How many events of type k1 will happen immediately after an event of type k2?"). Probabilistic forecasts are universally preferred to point estimates in the neighboring field of time series modeling [Gneiting and Katzfuss, 2014; Alexandrov et al., 2020]. A variety of metrics for evaluating the quality of such probabilistic forecasts for (regularly-spaced) time series have been proposed [Gneiting et al., 2008], but they haven't been generalized to marked continuous-time event sequences. Developing metrics that are based on entire sampled event sequences can help us unlock the full potential of neural TPPs and allow us to better compare different models. More generally, we should take advantage of the fact that TPP models learn an entire distribution over trajectories and rethink how these models are applied to prediction tasks.

7.3 Applications

While most recent works have focused on developing new architectures for better prediction and structure discovery, other applications of neural TPPs remain largely understudied.

Applications of classical TPPs go beyond the above two tasks. For example, latent-variable TPP models have been used for event clustering [Mavroforakis et al., 2017; Xu and Zha, 2017] and change point detection [Altieri et al., 2015]. Other applications include anomaly detection [Li et al., 2017], optimal control [Zarezade et al., 2017; Tabibian et al., 2019] and fighting the spread of misinformation online [Kim et al., 2018].

Moreover, most papers on neural TPPs consider datasets originating from the web and related domains, such as recommender systems and knowledge graphs. Meanwhile, traditional application areas for TPPs, like neuroscience, seismology and finance, have not received as much attention. Adapting neural TPPs to these traditional domains requires answering exciting research questions. To name a few, spike trains in neuroscience [Aljadeff et al., 2016] are characterized by both high numbers of marks (i.e., neurons that are being modeled) and high rates of event occurrence (i.e., firing rates), which requires efficient and scalable models. Applications in seismology usually require interpretable models [Bray and Schoenberg, 2013], and financial datasets often contain dependencies between various types of assets that are far more complex than self-excitation, commonly encountered in social networks [Bacry et al., 2015].

To summarize, considering new tasks and new application domains for neural TPP models is an important and fruitful direction for future work.