Fully Neural Network based Model for General Temporal Point Processes

Introduction

The activity of many diverse systems is characterized as a sequence of temporally discrete events. The examples include financial transactions, communication in a social network, and user activity at a web site. In many cases, the occurrences of the event are correlated to each other in a certain manner, and information on future events may be extracted from the information of past events. Therefore, the appropriate modeling of the dependence of the event occurrence on the history of past events is important for understanding the system and predicting future events.

A temporal point process is a useful mathematical tool for modeling the time series of discrete events. In this framework, the dependence on the event history is characterized using a conditional intensity function that maps the history of the past events to the intensity function of the point process. The most common models, such as the Poisson process or the Hawkes process , assume a specific parametric form for the conditional intensity function. Recently, Du et al. (2016) proposed a model based on a recurrent neural network (RNN) for point processes , and the variant models were further developed . In this approach, an RNN is used to obtain a compact representation of the event history. The conditional intensity function is then modeled as a function of the hidden state of the RNN. Consequently, the RNN based models outperform the parametric models in prediction performance.

Although such RNN based models aim to capture the dependence of the event occurrence on the event history in a general manner, a specific functional form is usually assumed for the time course of the conditional intensity function (see for exception). For example, the model in assumed that the conditional intensity function exponentially decreases or increases with the elapsed time from the most recent event until the next event. However, using such an assumption can limit the expressive ability of the model and potentially deteriorate the predictive skill if the employed assumption is incorrect. We herein generalize RNN based models such that the time evolution of the conditional intensity function is represented in a general manner. For this purpose, we formulate the conditional intensity function based on a neural network rather than assuming a specific functional form.

However, exactly evaluating the log-likelihood function for such a general model is generally intractable because the log-likelihood function of a temporal point process contains the integral of the conditional intensity function. Although some studies, which considered a general model of the intensity function, used numerical approximations to evaluate the integral , numerical approximations can deteriorate the fitting accuracy and can also be computationally expensive. To overcome this limitation, we first model the integral of the conditional intensity function using a feedforward neural network rather than directly modeling the conditional intensity function itself. Then, the conditional intensity function is obtained by differentiating it. This approach enables us to exactly evaluate the log-likelihood function of our general model without numerical approximations. Finally, we show the effectiveness of our proposed model by analyzing synthetic and real datasets.

Method

A temporal point process is a stochastic process that generates a sequence of discrete events at times $\{t_{i}\}_{i=1}^{n}$ in a given observation interval $[0,T]$. The process is characterized via a conditional intensity function $\lambda(t|H_{t})$, which is the intensity function of the event at the time $t$ conditioned on the event history $H_{t}=\{t_{i}|t_{i}<t\}$ up to the time $t$, given as follows:

If the conditional intensity function is specified, the probability density function of the time $t_{i+1}$ of the next event, given the times $\{t_{1},t_{2},\ldots,t_{i}\}$ of the past events, is obtained as follows:

where the exponential term in the right-hand side represents the probability that no events occur in $[t_{i},t_{i+1})$. The probability density function to observe an event sequence $\{t_{i}\}_{i=1}^{n}$ is then obtained as follows:

The most basic example of a temporal point process is a stationary Poisson process, which assumes that the events are independent of each other. The conditional intensity function of the stationary Poisson process is given as $\lambda(t|H_{t})=\lambda$. Another popular example is the Hawkes process , which is a simple model of a self-exciting point process. The conditional intensity function of the Hawkes process is given as $\lambda(t|H_{t})=\mu+\sum_{t_{i}<t}g(t-t_{i})$, where $g(s)$ is a kernel function ($g(s)=0$ if $s<0$) that represents the triggering effect from the past event.

2 Recurrent neural network approach to a temporal point process

The conditional intensity function, which maps the event history to the intensity function, plays a major role in the modeling of point processes. Du et al. (2016) proposed to use an RNN to model the conditional intensity function . In this approach, an input vector $\boldsymbol{x}_{i}$, which extracts the information of the event time $t_{i}$, is first fed into the RNN. A simple form of the input is the inter-event interval as $\boldsymbol{x}_{i}=(t_{i}-t_{i-1})$ or its logarithm as $\boldsymbol{x}_{i}=(\log(t_{i}-t_{i-1}))$. A hidden state $\boldsymbol{h}_{i}$ of the RNN is updated as follows:

where $W^{h}$, $W^{x}$, and $\boldsymbol{b}^{h}$ denote the recurrent weight matrix, input weight matrix, and bias term, respectively, and $f$ is an activation function. We here treat the hidden state of the RNN as a compact vector representation of the event history. The conditional intensity function is then formulated as a function of the elapsed time from the most recent event and the hidden state of the RNN, given as follows:

where $\phi$ is a non-negative function referred to as a hazard function.

Du et al. (2016) assumed the following form for the hazard function :

The exponential function in the above equation is used to ensure the non-negativity of the intensity. In this model, the conditional intensity function exponentially decreases or increases with the elapsed time $\tau$ from the most recent event until the next event.

A simplified model, where the conditional intensity function is constant over the period between the successive events, was also considered in . We here formulate such a model as a special case of the model of eq. (6), given as follows:

In this model, the inter-event interval $\tau_{i}=t_{i+1}-t_{i}$ follows the exponential distribution with mean $1/\exp(\boldsymbol{v}^{\phi}\cdot\boldsymbol{h}_{i}+b^{\phi})$.

The log-likelihood function of the RNN based model can be obtained from eq. (3) as follows:

The parameter values of the model are estimated by maximizing the log-likelihood function. For this purpose, the backpropagation through time (BPTT) is employed to obtain the gradient of the log-likelihood function. In the BPTT, the RNN is unfolded to a feedforward network, whose weights are shared across layers, and the backpropagation is applied to the unfolded network. Although the hidden state $\boldsymbol{h}_{i}$ of the RNN originally depends on all the preceding events, fully considering the history dependence for a long sequence can be problematic due to gradient vanishing or explosion. Only the dependence on a fixed number $d$ of the most recent events is considered herein as was done in , where $d$ is a hyperparameter called the truncation depth. Namely, for each time index $i$, the hidden state $\boldsymbol{h}_{i}$ is obtained by feeding the inputs $\{\boldsymbol{x}_{j}\}_{j=i-d+1}^{i}$ from the $d$ most recent events to the RNN.

3 Problem

The main problem of the previous studies is that a specific functional form is usually assumed for the time course of the hazard function $\phi(\tau|\boldsymbol{h}_{i})$ as in eqs. (6) or (7), which can miss the general dependence of the event occurrence on the past events. One may want to exploit a more complex model for the hazard function $\phi(\tau|\boldsymbol{h}_{i})$ to generalize the model. However, such a complex model is generally intractable because the log-likelihood function in eq. (8) includes the integral of the hazard function. Although the integral may be approximately evaluated using numerical methods , the numerical approximations can deteriorate the fitting accuracy and be computationally expensive. This is the main limitation in the flexible modeling of the hazard function.

4 The proposed model A source code is available online. https://github.com/omitakahiro/NeuralNetworkPointProcess

Rather than directly modeling the hazard function, we herein propose to model the cumulative hazard function $\Phi(\tau|\boldsymbol{h}_{i})$, defined as follows:

The hazard function itself can be then obtained by differentiating the cumulative hazard function with respect to $\tau$ as follows:

The log-likelihood function is reformulated as follows using the cumulative hazard function:

Now, the log-likelihood function does not include the integral term in contrast to eq. (8) and can be exactly evaluated even for a complex model of the cumulative hazard function.

In the present study, we model the cumulative hazard function using a feedforward neural network (a cumulative hazard function network; Fig. 1) for flexible modeling. The cumulative hazard function is a monotonically increasing function of $\tau$ and is positive-valued. The cumulative hazard function network is designed to reproduce these properties. The positivity of the network output can be ensured using an output unit, in which the activation function is positive-valued. Considering monotonicity, we employ the idea used in . To summarize, the weights of the particular network connections are constrained to be positive (Fig. 1).

The detail of the cumulative hazard function network is described below. In the network, each unit receives the weighted sum of the inputs and applies an activation function to produce the output. The first hidden layer in the network receives the elapsed time $\tau$ and the hidden state $\boldsymbol{h}_{i}$ of the RNN as the inputsWe may input $\log\tau$ to the first layer rather than $\tau$ if the variation of the inter-event interval is large.. The weights of the connections from the elapsed time $\tau$ to the first hidden layer and all the connections from the hidden layers are constrained to be positiveIn order to enforce the weights to be positive, if a weight is updated to be a negative value during training, we replace it with its absolute value.. The connections, in which the weights are constrained to be positive, are indicated by the red line in Fig. 1. The activation functions of the hidden units and the output unit are set to be the $tanh$ function and the $softplus$ function, $\log(1+\exp(\cdot))$, respectively. In this setting, the network output is monotonically increasing with respect to the elapsed time $\tau$ and takes only a positive value, which mimics the cumulative hazard function.

The cumulative hazard function $\Phi(\tau|\boldsymbol{h}_{i})$ and the hazard function $\phi(\tau|\boldsymbol{h}_{i})$ are now formulated as follows based on the output $Z_{i}(\tau)$ of the cumulative hazard function network:

The differentiation term $\partial Z_{i}(\tau)/\partial\tau$, which is the derivative of the network output with respect to the network input, is computed using automatic differentiation (see Supplementary Material for more details). Automatic differentiation is a method to calculate the derivative of an arbitrary function, and it can be easily carried out using neural network libraries such as TensorFlow and PyTorch. The log-likelihood function in eq. (11) of our model is then given based on the network output via the eqs. (12) and (13). The gradient of the log-likelihood function with respect to the parameters is obtained using backpropagation (see Supplementary Material for more details).

We note that our model can also efficiently generate a prediction of the timing of the coming event in the following way. The predictive probability density function $p^{*}(t|t_{1},t_{2},\ldots,t_{i})$ of the time $t_{i+1}$ of the coming event given the past events $\{t_{1},t_{2},\ldots,t_{i}\}$ is calculated from eq. (2). We herein use the median $t_{i+1}^{*}$ of the predictive distribution $p^{*}$ to predict $t_{i+1}$. To obtain the median $t_{i+1}^{*}$, we use the relation $\Phi(t_{i+1}^{*}-t_{i}|\boldsymbol{h}_{i})=\log(2)$. This relation is derived from the property that the integral of the intensity function over $[t_{i},t_{i+1}]$ follows the exponential distribution with mean 1, or is derived by directly integrating eq. (2). Then, the median $t_{i+1}^{*}$ can be efficiently obtained by solving the above relation using a root finding method (e.g., the bisection method); it takes only a second for our model to generate predictions for 20000 events. Therefore, the cumulative hazard function also plays a crucial role in generating a median predictor.

Related works

The RNN based point process models were proposed in . Most previous studies assumed a specific functional form for the time-course of the hazard function. The exponential hazard function in eq. (6) is commonly assumed . Some studies assumed the constant hazard function as in eq. (7), which is equivalent to assume that the inter-event intervals follow the exponential distribution. In contrast to these studies, our model does not assume any specific functional form for the hazard function, and the time course of the hazard function is formulated in a general manner based on a neural network.

A few studies addressed the general modeling of the hazard function. Jing and Somla (2017) proposed to discretize the continuous hazard function to a piecewise constant function , given as follows:

for $j=1,2,\ldots,\tau_{max}/l$ for some choice of $l$ and $\tau_{max}$. Mei and Eisner (2017) proposed a continuous-time long short-term memory (LSTM) where the output continuously evolves in time, and the conditional intensity function is given as a function of the output . This model used the Monte Carlo method to approximate the integral. In all cases, numerical approximations are used to evaluate the integral in the log-likelihood function in eq. (8). However, numerical approximations can be computationally expensive and can also affect the fitting accuracy. In contrast to these studies, the log-likelihood function of our general model can be exactly evaluated without any numerical approximations because the integral of the hazard function is modeled by a feedforward neural network in our approach. Therefore, a more accurate estimate can be efficiently obtained by our approach.

Experiments

In this section, we conduct experiments using synthetic and real data. We herein evaluate the predictive performances of the four RNN based point process models. The number of units in the RNN is fixed to 64 for all the models. The first two models are equipped with the constant hazard function in eq. (7) (the constant model) and the exponential hazard function in eq. (6) (the exponential model), respectively. The third model employs the piecewise constant hazard function in eq. (14) (the piecewise constant model). We set $\tau_{max}$ to the maximum value of the inter-event interval in each dataset and use the condition $l=\tau_{max}/128$. The fourth model employs the neural network based hazard function proposed in this study (the neural network based model). For this model, we use two hidden layers for the cumulative hazard function network, and the number of units in each layer is 64. In this setting, the numbers of the parameters are almost the same between the third and fourth models.

Each dataset is divided into the training and test data. In the training phase, the model parameters are estimated using the training data. For this optimization, the Adam optimizer with the learning rate $0.001$, $\beta_{1}=0.9$, and $\beta_{2}=0.999$ is used , and the batch size is 256. We also choose the truncation depth $d$ from $$using$20\%$of the training data (see Sec. 2.2 for more details on truncation depth). In the test phase, we evaluate the predictive performances of the trained models using the test data. In each time step, the probability density function$p^{*}(t|t_{1},t_{2},\ldots,t_{i})$of the time of the coming event given the past events$\{t_{1},t_{2},\ldots,t_{i}\}$is calculated from eq. (2), and scored by the negative log-likelihood$-\log p^{*}(t_{i+1}|t_{1},t_{2},\ldots,t_{i})$for the actually observed$t_{i+1}$ (a smaller score means a better predictive performance). The score is finally averaged over the events in the test data. We performed these computations under a GPU environment provided by Google Colaboratory.

We use the synthetic data generated from the following stochastic processes. In this experiment, 100,000 events are generated from each process, and 80,000/20,000 events are used for training/testing.

Stationary Poisson Process (S-Poisson): The conditional intensity function is given as $\lambda(t|H_{t})=1$.

Non-stationary Poisson Process (N-Poisson): The conditional intensity function is given as $\lambda(t|H_{t})=0.99\sin(2\pi t/20000)+1$.

Stationary Renewal Process (S-Renewal): In this process, the inter-event intervals $\{\tau_{i}=t_{i+1}-t_{i}\}$ are independent and identically distributed according to a given probability density function $p(\tau)$. We herein use the log-normal distribution with a mean of 1.0 and a standard deviation of 6.0 for $p(\tau)$. In this setting, a generated sequence exhibits a bursty behavior: multiple events tend to occur in a short period and are followed by a long silence period like the burst firing of biological neurons.

Non-stationary Renewal Process (N-Renewal): A sequence $\{t_{i}\}$ following a non-stationary renewal process is obtained as follows : we first generate a sequence $\{t^{\prime}_{i}\}$ from a stationary renewal process, and then we rescale the time according to $t^{\prime}_{i}=\int_{0}^{t_{i}}r(t)dt$ for a non-negative trend function $r(t)$. We use the gamma distribution with a mean of 1.0 and a standard deviation of 0.5 to generate the stationary renewal process and set the trend function to $r(t)=0.99\sin(2\pi t/20000)+1$. In this process, an inter-event interval tends to be followed by the one with similar length, but the expected length gradually varies in time.

Self-correcting Process (SC): The conditional intensity function is given as $\lambda(t|H_{t})=\exp(t-\sum_{t_{i}<t}1)$.

Hawkes Processes (Hawkes1 and Hawkes2): We use the Hawkes process, in which the kernel function is given by the sum of multiple exponential functions: the conditional intensity function is given by $\lambda(t|H_{t})=\mu+\sum_{t_{i}<t}\sum_{j=1}^{M}\alpha_{j}\beta_{j}\exp\{-\beta_{j}(t-t_{i})\}$. For the Hawkes1 model, we set $M=1,\mu=0.2,\alpha_{1}=0.8,\mbox{and }\beta_{1}=1.0$. For the Hawkes2 model, we set $M=2,\mu=0.2,\alpha_{1}=0.4,\beta_{1}=1.0,\alpha_{2}=0.4,\mbox{and }\beta_{2}=20.0$. Compared to the Hawkes1 model, the kernel function of the Hawkes2 model rapidly varies in time for small $\tau$.

In addition to the four RNN based models, we evaluate the predictive performance of the true model, i.e., the model that generated the data, and use it as a reference. The scores of the RNN based models are standardized by subtracting the score of the true model. The value of in the standardized score corresponds to the score of the true model.

Figure 2 summarizes the performances of the four RNN based models for the synthetic datasets (a smaller score means a better performance). We first find that the proposed neural network based model achieves a competitive or better performance against performances of the other models. The neural network based model also performs robustly for all the datasets: the performance of the neural network based model is always close to that of the true model. These results demonstrate that (i) the performance is improved by employing the neural network based hazard function and that (ii) our model can be applicable to a diverse class of data generating processes.

The performances of the constant model and the exponential model critically depend on whether the hazard function is correctly specified. The constant hazard function is correct for the S-Poisson and N-Poisson processes. The exponential hazard function is correct for the S-Poisson, N-Poisson, and SC processes, and is approximately correct for the Hawkes1 process (a constant term is included in the hazard function of the Hawkes1 process but not in the exponential hazard function). In fact, these models perform similarly to the true model for the cases where the hazard function is correctly specified but perform poorly for the other cases. Figure 3 shows the estimated conditional intensity function and clearly demonstrates that the exponential model captures well the true conditional intensity function for the self-correcting process where the exponential hazard function is valid; however, it fails for the Hawkes2 process where the exponential hazard function is not valid. In contrast, our neural network based model can reproduce well the true model for both cases. In this manner, the constant and exponential models are sensitive to model misspecification.

The performance of the piecewise constant model is much worse than the neural network based model, particularly for the S-Renewal, N-Renewal, and Hawkes2 processes. For these processes, the variability of the inter-event intervals is large, and the conditional intensity function can rapidly vary for a short period after an event. For such cases, the piecewise constant approximation might not work well. The performance of the piecewise constant model would be improved if the approximation accuracy is improved; however, this increases the computational cost. In this experiment, the numbers of the parameters are set to be almost the same between the piecewise constant model and the neural network based model. Moreover, the neural network based model performs better than the piecewise constant model, indicating that the neural network based model is more efficient than the piecewise constant model.

2 Real data

We use the following real datasets for the next experiment.

Finance dataset: This dataset contains the trading records of Nikkei 225 mini, which is the most liquid features contracts in Asia . The timestamps of 182,373 transactions in one day are analyzed, and the first $80\%$ and the last $20\%$ of the data are used for training and testing, respectively.

Emergency call dataset: This dataset contains the records of the police department calls for service in San Francisco . Each record contains the timestamp and address from which the call was made. We prepare 100 separate sequences for the 100 most frequent addresses, which contain a total of 294,865 events. The first $80\%$ and the last $20\%$ of the events in the sequences are used for training and testing, respectively.

Meme dataset: MemeTracker tracks the popular phrases from numerous online resources such as news media and personal blogs . This dataset records the timestamps when the focused phrases appear on the internet. We first extract the 50 most frequent phrases and obtain the corresponding 50 separate sequences. We use 40 sequences out of 50 with 247,579 events for training and the remaining 10 sequences with 61,095 events for testing.

Music dataset: This dataset records the history of music listening of users at https://www.last.fm/ . We prepare the 100 sequences for the 100 most active users in Jan 2009, which contain a total of 299,046 events. The first $80\%$ and the last $20\%$ of the events in the sequences are used for training and testing, respectively.

Figure 4 summarizes the performances for the real datasets. We find that the neural network based model exhibits a competitive or superior score as compared to those of the other models; this demonstrates the practical effectiveness of the proposed model. For the finance dataset, the performances of all the models are close to each other, implying that the constant or exponential hazard function reproduces the event occurrence process of the financial transactions. The neural network based model performs much better than the other models, particularly for the Meme and music datasets: the difference in the scores between the neural network based model and the other models is greater than 0.5. This difference should be significant (e.g., in the case of the right panel in Fig. 3, where the exponential model clearly fails to reproduce the true model, the score difference is about 0.4 between the true and exponential models). For the two datasets, the data contain inter-event intervals that are much longer than the average, and the variability of the inter-event intervals is large. Other than the neural network based model, the models presumably fail to adapt to such a feature.

3 Time prediction experiment

To evaluate the predictive performances based on a metric other than the log-likelihood, we also carry out the time prediction experiments. Specifically, we use the median of the predictive distribution to predict the timing of the coming event and evaluate the prediction by the mean absolute error. The result is summarized in the table below. In the table, the best score is in bold, and it is in red if the difference between the best score and the second best score is statistically significant ($p<0.01$). We find that our model performs better than the other models on average and that the performance of our model is best or close to the best for the most individual datasets. These results demonstrate the effectiveness of our model in the prediction task.

4 Comparison with the continuous-time LSTM model

We here compare our model with the continuous-time LSTM (CT-LSTM) model . Both the two models aim at flexibly estimating the intensity function. The technical advantage of our model over the CT-LSTM model is that the log-likelihood function can be exactly evaluated, therefore the estimation can be carried out efficiently. Our model is also relatively easy to implement. For the CT-LSTM model, the evaluation of the log-likelihood function is based on the Monte Carlo method, which can be computationally expensive and can deteriorate the performance.

We evaluate the predictive performance of the CT-LSTM model in terms of the mean negative log-likelihood (MNLL) and the mean absolute error (MAE). The number of the hidden units in the CT-LSTM model is set to 42, so that the numbers of the free parameters are almost the same between the CT-LSTM model and our model. The following table lists the scores of the CT-LSTM model relative to our model; the positive score means that our model is better than the CT-LSTM model. In the table, the score is in blue if our model is significantly better than the CT-LSTM model ($p<0.01$). The performance of our model is better than the CT-LSTM model on average, demonstrating the effectiveness of our model.

Discussion and Conclusions

In this study, we extended the RNN based point process models such that the time course of the hazard function is represented in a general manner based on a neural network. We then showed the effectiveness of our model by analyzing both synthetic and real datasets. Primary advantages of the proposed model are summarized as follows:

By using a feedforward neural network, our model can reproduce any time course of the hazard function in principle, i.e., the usefulness of fully neural network based modeling for point processes is indicated.

By modeling the cumulative hazard function rather than the hazard function itself, we can avoid the direct evaluation of the integral in the log-likelihood function. The log-likelihood function can be exactly and efficiently evaluated without relying on the numerical approximation because of this approach.

We note that the cumulative hazard function also plays an important role in the diagnostic analysis . We did not consider herein the marks of each event (i.e., the information associated with each event other than the timestamp) because the primary contribution of this work is the development of a general model of the hazard function. However, our approach can be easily extended to marked temporal point processes as well.

This research was partly supported by AMED under Grant Number JP19dm0307009. T. O. and K. A. are supported by Kozo Keikaku Engineering Inc.

References

Supplementary Material for "Fully Neural Network based Model for General Temporal Point Processes"

where $f^{(j)}$, $W^{(j)}$, and $\boldsymbol{b}^{(j)}$ represent the activation function, the weight matrix, and the bias term for the $j$th layer, respectively. The input to the cumulative hazard function network is given as $\boldsymbol{Y}^{(0)}=(\tau,\boldsymbol{h}_{i}^{T})^{T}$, where $\boldsymbol{h}_{i}$ is the RNN output. The $L$th layer is the output layer, and we have $Z_{i}(\tau)=\boldsymbol{Y}^{(L)}$.

To calculate $\partial Z_{i}(\tau)/\partial\tau$, we introduce a notation

and employ a chain rule to obtain a backward operation for $\boldsymbol{y}^{(j)}$ as

where we have $\boldsymbol{y}^{(L)}=1$. By recursively applying the backward operation, we finally obtain the desired derivative as

where $[\boldsymbol{y}^{(0)}]_{1}$ represents the first element of the vector $\boldsymbol{y}^{(0)}$. We can then construct an extended network that produces both $Z_{i}(\tau)$ and $\partial Z_{i}(\tau)/\partial\tau$ by concatenating the original feedforwad network of eq. (15) and the backward network of eq. (18), as shown in Fig. S1.

Evaluating the gradient of the loss function

The loss function of our model, the negative log-likelihood function, depends on both $Z_{i}(\tau)$ and $\partial Z_{i}(\tau)/\partial\tau$. The derivative $\partial Z_{i}(\tau)/\partial\tau$ is computed using backpropagation, as seen in the previous section. Then the computational graph to evaluate the loss function is then obtained as in Fig. S1. The gradient of the loss function with respect to the parameters can be evaluated by applying backpropagation to the computational graph.

Here we use backpropagation twice for training the model, the first for calculating $\partial Z_{i}(\tau)/\partial\tau$ and the second for calculating the gradient of the loss function. This kind of procedure is sometimes called as double backpropagaton .