A neural network based model for multi-dimensional nonlinear Hawkes
  processes
This paper introduces the Neural Network for Nonlinear Hawkes processes
(NNNH), a non-parametric method based on neural networks to fit nonlinear
Hawkes processes. Our method is suitable for analyzing large datasets in which
events exhibit both mutually-exciting and inhibitive patterns. The NNNH
approach models the individual kernels and the base intensity of the nonlinear
Hawkes process using feed forward neural networks and jointly calibrates the
parameters of the networks by maximizing the log-likelihood function. We
utilize Stochastic Gradient Descent to search for the optimal parameters and
propose an unbiased estimator for the gradient, as well as an efficient
computation method. We demonstrate the flexibility and accuracy of our method
through numerical experiments on both simulated and real-world data, and
compare it with state-of-the-art methods. Our results highlight the
effectiveness of the NNNH method in accurately capturing the complexities of
nonlinear Hawkes processes.
Sobin Joseph, Shashi Jain
# A neural network based model for multi-dimensional nonlinear Hawkes processes

###### Abstract

This paper introduces the Neural Network for Nonlinear Hawkes processes (NNNH), a non-parametric method based on neural networks to fit nonlinear Hawkes processes. Our method is suitable for analyzing large datasets in which events exhibit both mutually-exciting and inhibitive patterns. The NNNH approach models the individual kernels and the base intensity of the nonlinear Hawkes process using feed forward neural networks and jointly calibrates the parameters of the networks by maximizing the log-likelihood function. We utilize Stochastic Gradient Descent to search for the optimal parameters and propose an unbiased estimator for the gradient, as well as an efficient computation method. We demonstrate the flexibility and accuracy of our method through numerical experiments on both simulated and real-world data, and compare it with state-of-the-art methods. Our results highlight the effectiveness of the NNNH method in accurately capturing the complexities of nonlinear Hawkes processes.

KEYWORDS: nonlinear Hawkes processes, Hawkes processes with inhibition, neural networks for Hawkes Process, online learning for Hawkes processes

## 1 Introduction

A.G. Hawkes introduced the Hawkes process, which is a type of multivariate point process used to model a stochastic intensity vector that depends linearly on past events Hawkes (1971). This approach has been widely applied in various fields, such as seismology Ogata (1999), Marsan and Lengline (2008), financial analysis Filimonov and Sornette (2012), Bacry et al. (2015), social interaction modeling Crane and Sornette (2008), Blundell et al. (2012), Zhou et al. (2013), criminology Mohler et al. (2011), genome analysis Reynaud-Bouret et al. (2010), Carstensen et al. (2010), and epidemiology Park et al. (2020), Chiang et al. (2022).

The Hawkes process in its original form is linear, i.e., the intensities depend on past events through a linear combination of kernel functions. The primary concern in modelling the Hawkes process is estimating these kernel functions. A common practice has been to assume a parametric form for the kernel function, such as the exponential and power-law decay kernels, and then use maximum likelihood estimation (Ozaki, 1979) to determine the optimal values of the parameters.

While easier to estimate, a parametric kernel function is often not expressive enough for real-world problems. A significant body of work is concerned with the non-parametric estimation of the kernel functions. Lewis and Mohler (2011) propose a non-parametric, Expectation-Maximization (EM) based method to estimate the Hawkes kernel and the non-constant base intensity function. Zhou et al. (2013) uses the method of multipliers and majorization-minimization approach to estimate the multivariate Hawkes kernels. Bacry and Muzy (2014) develop a non-parametric estimation
based on the Wiener-Hopf equations for estimating the Hawkes kernels. Achab et al. (2017) derive an integrated cumulants method for estimating the Hawkes kernel. Xu et al. (2016) proposed a sparse-group-lasso-based algorithm combined with likelihood estimation for estimating the Hawkes kernels. Joseph et al. (2022) use a feed-forward network to model the excitation kernels and fit the network parameters using a maximum likelihood approach. Except for Bacry and Muzy (2014), these non-parametric kernel estimation methods are generally restricted to the linear Hawkes processes.

Modelling the intensities as a linear combination of kernel functions imposes a non-negativity constraint on the kernel functions, which can be interpreted as excitation kernels. Non-negative kernel functions do not allow incorporating inhibitive effects while modelling the intensity process. Modelling inhibition, where the occurrence of an event reduces the intensity of future arrival, has drawn relatively less attention in the literature. However, inhibitory effects are prevalent in various domains. For example, in neuroscience, inhibitory kernel functions can represent the presence of a latency period before the successive activations of a neuron (Reynaud-Bouret et al., 2013).

This paper proposes a non-parametric estimation model for the non-linear Hawkes process. The non-linear Hawkes process allows the inclusion of both excitatory and inhibitory effects to model a broader range of phenomena. The stability condition of the nonlinear Hawkes process is explored in Bremaud and Massoulie (1996). Bonnet et al. (2021) use a negative exponential function to model inhibitive kernels for a univariate nonlinear Hawkes process and use maximum likelihood estimation to determine the optimal parameters. Bonnet et al. (2022) extend this approach to a multivariate nonlinear Hawkes process. Lemonnier and Vayatis (2014) develop the Markovian Estimation of Mutually Interacting Process (MEMIP) method, which utilizes weighted exponential functions to determine kernels for the nonlinear Hawkes process. To extend the MEMIP for large dimensional datasets, Lemonnier et al. (2016) introduce dimensionality reduction features. However, the above approaches require the kernel function to be smooth, which is a drawback. Wang et al. (2016) propose an algorithm that learns the nonlinear Hawkes kernel non-parametrically using isotonic regression and a single-index model. The Isotonic Hawkes Process, however, assumes only a continuous monotonic decreasing excitation kernel to capture the influence of the past arrivals.

In their work, Du et al. (2016) employ a recurrent neural network (RNN) to model event timings and markers simultaneously by leveraging historical data. Conversely, Mei and Eisner (2017) propose a novel continuous-time long short-term memory (LSTM) model to capture the self-modulating Hawkes processes, capable of accounting for the inhibiting and exciting effects of prior events on future occurrences. Moreover, this approach can accommodate negative background intensity values, which correspond to delayed response or inertia of some events. While RNN-based models can capture complex long-term dependencies, inferring causal relationships or extracting causal information from these models can be challenging.

We propose a neural networks based non-parametric model for the nonlinear Hawkes process where feed-forward networks are used to model the kernels and time-varying base intensity functions. The architecture of the neural network is chosen such that the likelihood function and its gradients with respect to the network parameters can be efficiently evaluated. As the likelihood function is non-convex with respect to the parameter space, we use the Stochastic Gradient Descent (SGD) with Adam (Kingma and Ba, 2014) to obtain the optimal network parameters that maximize the log-likelihood. The method is an extension of the Shallow Neural Hawkes model proposed by Joseph et al. (2022), which was designed for linear Hawkes processes and only allows for the modelling of kernels with an excitation feature. We evaluate our model against state-of-the-art methods for non-linear Hawkes processes, using both simulated and real datasets.

The paper is structured as follows: Section 2 provides a definition of the non-linear Hawkes process and formulates the associated log-likelihood maximization problem. In Section 3, we introduce the proposed neural network model for the non-linear Hawkes process and discuss the parameter estimation procedure. Section 4 presents the results obtained by our method on synthetic data and its application to infer the order dynamics of bitcoin and ethereum market orders on the Binance exchange. Additionally, we apply our method to a high-dimensional neuron spike dataset. Finally, in Section 5, we summarize our findings and discuss the limitations of our method.
Preliminary Definitions

Definition of the Hawkes Process:We denote a D-dimensional counting process as \(N(T)\), and its associated discrete events as \(\mathcal{S}=(t_{n},d_{n})_{n\geq 1}\), where \(t_{n}\in[0,T)\) is the timestamp of the \(n\)th event and \(d_{n}\in(1,2,...,D)\) is the dimension in which the \(n\)th event has occurred. Let \(\{t_{1}^{d},\ldots,t_{k}^{d},\ldots,t_{N_{d}(T)}^{d}\}\), be the ordered arrivals for the dimensions \(d=1,\ldots,D\). Given \(t\geq 0\), the count of events in \([0,t)\) for the \(d\)th dimension will be

\[N_{d}(t)=\sum_{n\geq 1}\mathbb{1}_{t_{d}^{\ast}<t}.\]

The conditional intensity function for the counting process at \(t\), \(\lambda_{d}^{\ast}(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}\), is given by,

\[\lambda_{d}^{\ast}(t)=\underset{\Delta t\to 0}{lim}\frac{\mathbb{E} \left[N_{d}(t+\Delta t)-N_{d}(t)\right]\mathcal{H}_{t^{-}}}{\Delta t} \tag{1}\]

where, \(\mathcal{H}_{t^{-}}\), denotes the history of the counting process upto \(t\).

For a \(D\)-dimensional Hawkes process the conditional intensity \(\lambda_{d}(t),\) for the \(d\)-th dimension is expressed as,

\[\lambda_{d}(t)=\mu_{d}(t)+\sum_{j=1}^{D}\sum_{\{\forall k|(t_{k}^{j}<t)\}} \phi_{dj}(t-t_{k}^{j}), \tag{2}\]

where \(\mu_{d}(t)\), the exogenous base intensity for the \(d\)-th dimension, does not depend upon the history of the past arrivals. \(\phi_{dj}(t-t_{k}^{j}),\)\(1\leq d,j\leq D,\) are called the excitation kernels that quantify the magnitude of excitation of the conditional intensity for the \(d\)-th dimension at time \(t\) due to the past arrival at \(t_{k}^{j},\)\(\{\forall k|t_{k}^{j}<t\}\) in the \(j\)th dimension. These kernel functions are positive and causal, i.e., their support is in \(\mathbb{R}^{+}\). Inferring a Hawkes process requires estimating the base intensity function \(\mu_{d}\) and its kernels functions \(\phi_{dj},\) either by assuming a parametric form for the kernels or in a non-parametric fashion. A more generic form of the Hawkes process conditional intensity \(\lambda_{d}^{\ast}(t),\) is given by,

\[\lambda_{d}^{\ast}(t)=\Psi_{d}(\lambda_{d}(t)),\ \forall d=1,...,D, \tag{3}\]

where \(\Psi_{d}:\mathbb{R}\rightarrow\mathbb{R}^{+}\).

If \(\Psi_{d}\) is an identity function then, Equation 3 is equivalent to the linear Hawkes process expressed in Equation 2. However, if \(\Psi_{d}\) is a non-linear function, then the process is called a Non-linear Hawkes Process. For the stability and stationarity of the nonlinear Hawkes process, \(\Psi_{d}\) should be Lipschitz continuous, as given by Theorem 7 of Bremaud and Massoulie (1996). The advantage of the nonlinear Hawkes process is that it allows the kernel output to take negative values to model inhibitory effects as well as positive values to model excitatory effects, i.e., \(\phi_{dj}:\mathbb{R}^{+}\rightarrow\mathbb{R}\). Wang et al. (2016) uses a sigmoid function \((1+e^{-x})^{-1}\) and a decreasing function \(1-(1+e^{-x})^{-1},\) while Bonnet et al. (2021) and Costa et al. (2020) use a \(\max\) function as \(\Psi\).

The Log-Likelihood Function of non-linear Hawkes ProcessThe categorization of methods used for estimating Hawkes processes, as discussed in Cartea et al. (2021), can be broadly classified into three groups:

* _Maximum likelihood estimation (MLE)_ It is the most commonly used approach for estimating the kernels of the Hawkes process, as employed by many methods such as Bonnet et al. (2021), Du et al. (2016), and Lemonnier et al. (2016). However, MLE methods have high computational costs as their time complexity increases quadratically with the number of arrivals. Expectation maximization methods, such as those used by Lewis and Mohler (2011) and Zhou et al. (2013), also belong to this category.
* _Method of moments_Wang et al. (2016) and Achab et al. (2017) use a moment matching method to estimate the Hawkes process. These methods typically rely on the spectral properties of the Hawkes process. The WH method (Bacry and Muzy, 2014), which is used as a benchmark model in numerical experiments, estimates the Hawkes process as a system of Wiener-Hopf equations and also falls under this category.
* _Least squares estimation (LSE)_ This method involves minimizing the L2 error of the integrated kernel. Cartea et al. (2021) propose an efficient stochastic gradient-based estimation method for linear multivariate Hawkes processes, which applies to large datasets and general kernels. However, the method is not applicable suitable for nonlinear Hawkes process.

In this paper the optimal kernels are estimated through the maximum likelihood approach.

For a \(D\)-multidimensional nonlinear Hawkes process with \(\mathbf{\mu}=[\mu_{d}(t)]_{D\times 1}\) and \(\mathbf{\Phi}=[\phi_{dj}(t)]_{D\times D}\) ; \(1\leq d,j\leq D\), let \(\mathbf{\theta}=[\theta_{1},\ldots\theta_{p}\ldots\theta_{P}]\) be the set of parameters used to model \(\mathbf{\mu}\) and \(\mathbf{\Phi}\). The parameters can be estimated by maximizing the log-likelihood function over the sampled events from the process. The log-likelihood (LL) function corresponding to the dataset \(\mathcal{S}\) for the nonlinear Hawkes process is given by (see for instance Rubin (1972); Daley and Vere-Jones (2007)),

\[\mathcal{L}(\mathbf{\theta},\mathcal{S})=\sum_{d=1}^{D}\left[\sum_{t_ {n},d_{n}\in\mathcal{S}}\left\{\log(\lambda_{d}^{*}(t_{n}))-\int_{t_{n-1}^{d }}^{t_{n}^{d}}\lambda_{d}^{*}(s)ds\right\}\mathbb{1}_{d_{n}=d}\right] \tag{4}\]

Depending on the parametric form of the kernel,\(\phi_{dj}(t)\), the LL function may or may not be concave. For instance, even for the exponential kernel, \(\phi_{dj}(t)=\alpha_{dj}\exp(-\beta_{dj})(t)\), the LL function is not concave in the parameter space. We use the Stochastic Gradient Descent (SGD) method to search the local optima for the LL function in the parameter space. In order to estimate the parameters \(\mathbf{\theta}\) using the SGD we need an unbiased estimator of the gradient of the LL, as given by Equation 4, with respect to \(\theta\).

\[\nabla_{\theta_{p}}\left(\mathcal{L}(\mathbf{\theta},\mathcal{S})\right) = \nabla_{\theta_{p}}\left(\sum_{d=1}^{D}\left[\sum_{t_{n},d_{n}\in \mathcal{S}}\left\{\log(\lambda_{d}^{*}(t_{n}))-\int_{t_{n-1}^{d}}^{t_{n}^{d}} \lambda_{d}^{*}(s)ds\right\}\mathbb{1}_{d_{n}=d}\right]\right) \tag{5}\] \[= \sum_{d=1}^{D}\left[\sum_{t_{n},d_{n}\in\mathcal{S}}\nabla_{ \theta_{p}}\left(\left\{\log(\lambda_{d}^{*}(t_{n}))-\int_{t_{n-1}^{d}}^{t_{n} ^{d}}\lambda_{d}^{*}(s)ds\right\}\right)\mathbb{1}_{d_{n}=d}\right]\]

Following Equation 5, an unbiased estimator of the gradient of \(\mathcal{L}\) with respect to parameter \(\theta_{p}\) will be,

\[\nabla_{\theta_{p}}\widehat{\mathcal{L}}(\mathbf{\theta},t_{n}^{d}) :=\nabla_{\theta_{p}}\left(\log(\lambda_{d}^{*}(t_{n}^{d}))-\int_{t_{n-1}^{d}} ^{t_{n}^{d}}\left(\lambda_{d}^{*}(s)\right)ds\right), \tag{6}\]

where \(t_{n}^{d}\) is randomly sampled from \(\mathcal{S}\). The proposed model does not assume a specific parametric form for the kernels and the base intensities. The model and its estimation is discussed in Section 3.

## 3 Proposed Model

Hornik et al. (1989) show that a multilayer feed-forward networks with as few as one hidden layer are capable of universal approximation in a very precise and satisfactory sense. We here propose a neural network-based system called the Neural Network for Non-Linear Hawkes (NNNH) to estimate the kernels and the base intensity of the nonlinear Hawkes process. Each kernel, \(\phi_{dj},\,1\leq d,j\leq D,\) and base intensity function \(\mu_{d}\), is separately modelled as a feed-forward network. However, the parameters of all the networks are jointly estimated by maximizing the log-likelihood of the training dataset.

In the NNNH, each kernel \(\phi_{dj}\), \(1\leq d,j\leq D\) is modelled as a feed-forward neural network \(\hat{\phi}_{dj}:\mathbb{R}^{+}\rightarrow\mathbb{R}\),

\[\hat{\phi}_{dj}(t)=A_{2}\circ\varphi\circ A_{1}, \tag{7}\]

where \(A_{1}\): \(\mathbb{R}\rightarrow\mathbb{R}^{p}\) and \(A_{2}\): \(\mathbb{R}^{p}\rightarrow\mathbb{R}\). More precisely,

\[A_{1}(x)=\mathbf{W}_{1}x+\mathbf{b_{1}}\ \ \text{for}\ x\in\mathbb{R},\ \mathbf{W}_{1}\in\mathbb{R}^{p\times 1}, \mathbf{b_{1}}\in\mathbb{R}^{p},\]
\[A_{2}(\mathbf{x})=\mathbf{W}_{2}\mathbf{x}+b_{2}\,\,\,\text{for}\,\,\mathbf{x} \in\mathbb{R}^{p},\,\,\mathbf{W}_{2}\in\mathbb{R}^{1\times p},b_{2}\in\mathbb{R},\]

and \(\varphi:\mathbb{R}^{j}\rightarrow\mathbb{R}^{j},j\in\mathbb{N}\) is the component-wise Rectified Linear Unit (ReLU) activation function given by:

\[\varphi(x_{1},\ldots,x_{j}):=\left(\max(x_{1},0),\ldots,\max(x_{j},0)\right).\]

We take \(\mathbf{W}_{1}=[a_{1}^{1},a_{1}^{2},....,a_{1}^{p}]^{\top}\), \(\mathbf{W}_{2}=[a_{2}^{1},a_{2}^{2},....,a_{2}^{p}]\) and \(\mathbf{b_{1}}=[b_{1}^{1},b_{1}^{2},....,b_{1}^{p}]^{\top}.\) Therefore, the kernel function can be written as:

\[\widehat{\phi}_{dj}(x)=b_{2}+\sum_{i=1}^{p}a_{2}^{i}\max\left(a_{1}^{i}x+b_{1 }^{i},0\right) \tag{8}\]

As the kernel, \(\widehat{\phi}_{dj},\) maps from positive real numbers to real numbers, it is versatile in its ability to model both inhibitory and excitatory effects.

With a choice of \(p\) neurons for the hidden layer, the dimension of the parameter space for the above network will be \(3p+1.\) For a \(D\)-dimensional nonlinear Hawkes process there will be \(D^{2}\) kernels and the total number of parameters to be estimated would be \((3p+1)D^{2}.\) In case the base intensity is not constant, we also model it as a feed-forward neural network, with an architecture similar to the ones used for the kernels. We then need to approximate in total \((3p+1)\left(D^{2}+D\right)\) parameters. If the base intensity is assumed to be a constant, the total number of parameters to be estimated would be \((3p+1)D^{2}+D.\) The Figure 2 is a schematic representation of the NNNH model, with each \(\phi_{dj},\) and \(\mu_{d}\) being a network (as expressed in Figure 1). In order to estimate \(\lambda_{d}^{*}(t),\) we need in addition to \(t,\) the history of arrivals until \(t,\) i.e., \(\mathcal{H}_{t}^{-}.\) The NNNH therefore estimates \(\lambda_{d}^{*}(t)\) as:

\[\widehat{\lambda}_{d}^{*}(t)=\max\left(\widehat{\mu}_{d}(t)+\sum_{j=1}^{D} \sum_{\{\forall k\mid(t_{k}^{j}<t)\}}\widehat{\phi}_{dj}(t-t_{k}^{j}),0\right). \tag{9}\]

The outer \(\max\) function ensures that the estimated \(\lambda_{d}^{*}(t)\) is positive and mimics the \(\Psi\) function in the non-linear Hawkes process, as given in Equation 3.

Figure 1: The architecture of the feed-forward neural network used for modeling kernels and base intensity functions of the nonlinear Hawkes process in NNNH. Here \(\varphi,\)\(f\) are the activation functions. In the NNNH we take \(\varphi\) as ReLU and \(f\) as the identity function.
This non-parametric method employs a feed-forward neural network that can precisely estimate any continuous kernel function and base intensity function defined on a compact set, with a desired level of accuracy (Leshno et al. (1993)). However, the challenge is to jointly estimate the parameters of the \(D^{2}+D\) neural networks used to model the base intensity and the kernel functions. To do this, the set of parameter values that maximize the log-likelihood over the observed dataset \(\mathcal{S}\) is found. Since the log-likelihood function is non-convex in the parameter space, the batch stochastic gradient descent is used to estimate the parameter values that provide a local maxima of the log-likelihood function The gradient of the log-likelihood function as given by Equation 6 with respect to each network parameter is computed for updating the parameter values in the batch stochastic gradient method.

### Estimating the parameters of the NNNH

The challenge for the NNNH model is, given the dataset \(\mathcal{S}\) to efficiently obtain the parameters of the networks that locally maximize the log-likelihood function. More precisely, let \(\boldsymbol{\theta}=[\theta_{1},\ldots\theta_{p}\ldots\theta_{P}]^{\top}\) denote the parameters of the network and let \(\boldsymbol{\Theta}\) denote the parameter space; then we want to find, for given dataset \(\mathcal{S},\) the values of the parameters that maximize the log-likelihood function over the parameter space, i.e.,

\[\widehat{\boldsymbol{\theta}}=\operatorname*{arg\,max}_{\boldsymbol{\theta} \in\boldsymbol{\Theta}}\mathcal{L}\left(\boldsymbol{\theta},\mathcal{S}\right). \tag{10}\]

Gradient descent is an effective optimization technique when the log-likelihood function can be differentiated with respect to its parameters. This is because computing the first-order partial derivatives of all parameters has the same computational complexity as evaluating the log-likelihood function, making it a relatively efficient approach. With neural networks, stochastic gradient descent, has been shown to be more effective and efficient method for optimization. The NNNH parameter estimation algorithm utilizes Adam (Kingma and Ba (2014)), a stochastic gradient descent technique that adjusts learning rates adaptively and relies solely on first-order gradients. Equation 6 provides

Figure 2: NNNH model for multivariate nonlinear Hawkes process.
the necessary unbiased estimator of the log-likelihood function for the NNNH method to be used with the SGD. Algorithm 1 gives the pseudo-code for the NNNH parameter estimation using Adam.

```
0: dataset \(\mathcal{S}\) Input: learning rate \(\alpha=0.01\), decay rates \(\beta_{1}=0.9,\beta_{2}=0.999\), small constant \(\epsilon=10^{-8}\) Output: Optimal parameter value \(\widehat{\mathbf{\theta}}\)  initialize parameter vector \(\mathbf{\theta}_{\mathbf{0}}\); \(m_{0}\gets 0\) (initialize first moment vector); \(v_{0}\gets 0\) (initialize second moment vector);  initialize iteration step \(i=0\); whilestopping criterion not metdo \(i\gets i+1\)  randomly sample \((t_{i},d_{i})\in\mathcal{S}\) \(g_{i}\leftarrow\nabla_{\mathbf{\theta}}\widehat{\mathcal{L}}(\mathbf{\theta}_{i-1},t_{ i}^{d})\) compute unbiased gradient of the likelihood function \(m_{i}\leftarrow\beta_{1}m_{i-1}+(1-\beta_{1})g_{i}\) update biased first moment estimate \(v_{i}\leftarrow\beta_{2}v_{i-1}+(1-\beta_{2})g_{i}^{2}\) update biased second moment estimate \(\widehat{m}_{i}\leftarrow\frac{m_{i}}{1-\beta_{1}^{2}}\) compute bias-corrected first moment estimate \(\hat{v}_{i}\leftarrow\frac{v_{i}}{1-\beta_{2}^{2}}\) compute bias-corrected second moment estimate \(\mathbf{\theta}_{i}\leftarrow\mathbf{\theta}_{i-1}-\alpha\frac{\widehat{m}_{i}}{\sqrt{ v_{i}+\epsilon}}\) update parameters  end while
```

**Algorithm 1**Estimation of parameters of NNNH using Adam

The essential step while estimating the parameters for the NNNH model is computing the unbiased gradient of the likelihood function, i.e.,

\[\nabla_{\theta_{p}}\widehat{\mathcal{L}}(\mathbf{\theta},t_{n}^{d}):=\nabla_{ \theta_{p}}\left(\log(\widehat{\lambda}_{d}^{*}(t_{n}^{d}))-\int_{t_{n-1}^{d}} ^{t_{n}^{d}}\left(\widehat{\lambda}_{d}^{*}(s)\right)ds\right).\]

This in turn involves computing the following gradients,

\[\nabla_{\theta_{p}}\log(\widehat{\lambda}_{d}^{*}(t_{n}^{d}))\text{ and }\nabla_{\theta_{p}}\int_{t_{n-1}^{d}}^{t_{n}^{d}}\left(\widehat{\lambda}_{d}^{*} (s)\right)ds.\]

While it is possible to compute the former gradient efficiently, the challenge lies in accurately estimating the latter gradient. The former gradient calculation involves the following steps:

\[\nabla_{\theta_{p}}\log(\widehat{\lambda}_{d}^{*}(t_{n}^{d})) = \frac{1}{\widehat{\lambda}^{*}(t_{n}^{d})}\nabla_{\theta_{p}} \max\left(\widehat{\mu}_{d}(t_{n}^{d})+\sum_{j=1}^{D}\sum_{\{\forall k|(t_{k}^ {j}<t_{n}^{d})\}}\widehat{\phi}_{dj}(t_{n}^{d}-t_{k}^{j}),0\right),\] \[= \frac{1}{\widehat{\lambda}^{*}(t_{n}^{d})}\left(\nabla_{\theta_{ p}}\left(\widehat{\mu}_{d}(t_{n}^{d})\right)+\sum_{j=1}^{D}\sum_{\{\forall k|(t_{k}^ {j}<t_{n}^{d})\}}\nabla_{\theta_{p}}\left(\widehat{\phi}_{dj}(t_{n}^{d}-t_{k}^ {j})\right)\right)\mathbbm{1}_{\widehat{\lambda}_{d}(t_{d}^{d})>0}.\]

Here,

\[\widehat{\lambda}_{d}^{*}(t)=\max\left(\widehat{\lambda}_{d}(t),0\right).\]

In Appendix A, you can find the expressions for the gradients of the neural networks \(\widehat{\mu}_{d}(t)\) and \(\widehat{\phi}_{dj}(t).\) In Section 3.2 we provide an approach for efficient estimation of \(\nabla_{\theta_{p}}\int_{t_{n-1}^{d}}^{t_{n}^{d}}\left(\widehat{\lambda}_{d}^ {*}(s)\right)ds.\)
### Gradient of the integrated Hawkes intensity function

The computation of the first order derivatives of,

\[\widehat{\Lambda}_{d}(t_{n}^{d}):=\int_{t_{n-1}^{d}}^{t_{n}^{d}}\widehat{\lambda }_{d}^{*}(s)ds,\]

with respect to all the parameters will have the same computational complexity as evaluating the function itself. In the SGD method, it is crucial to perform efficient computations of \(\widehat{\Lambda}_{d}(t_{n}^{d})\) since these gradients must be computed multiple times. Gaussian Quadrature is a numerical technique that can be used to estimate the definite integral, but it is limited to well-behaved integrands that can be approximated by a certain degree of polynomial. Unfortunately, the non-linear Hawkes intensity function lacks smoothness, making the use of Gaussian Quadrature prone to substantial error. Moreover, Gaussian Quadrature entails computing the quadrature points and weights, which can be time-consuming and resource-intensive. The number of quadrature points required for accurate approximation also grows quickly with the polynomial degree, resulting in increased computational cost.

Our selection of network architecture can enable efficient evaluation of the integral \(\widehat{\Lambda}_{d}(t_{n}^{d})\) in the NNNH. The above integral can be expanded as:

\[\int_{t_{n-1}^{d}}^{t_{n}^{d}}\widehat{\lambda}_{d}^{*}(s)ds = \int_{t_{n-1}^{d}}^{t_{n}^{d}}\max\left(\widehat{\lambda}_{d}(s), 0\right)ds\] \[= \int_{t_{n-1}^{d}}^{t_{n}^{d}}\max\left(\left[\widehat{\mu}_{d}(s )+\sum_{j=1}^{D}\sum_{\{t_{k}^{j}<s\}}\widehat{\phi}_{dj}(s-t_{k}^{j})\right], 0\right)ds\] \[= \int_{t_{n-1}^{d}}^{t_{n}^{d}}\max\left(\left[\widehat{\mu}+\sum_ {j=1}^{D}\sum_{\{t_{k}^{j}<s\}}\sum_{i=1}^{p}a_{2}^{i}\max\left(a_{1}^{i}(s-t_{ k}^{j})+b_{1}^{i},0\right)\right],0\right)ds,\]

with \(\widehat{\phi}_{dj}(s)\) substituted using Equation 8 and \(\widehat{\mu}_{d}(s)\) set as a constant for simplicity in the third line. In addition, we set the bias \(b_{2}\) to zero for ease of exposition. Therefore, the above integral takes the following form:

\[\int_{t_{n-1}^{d}}^{t_{n}^{d}}\widehat{\lambda}_{d}^{*}(s)ds=\int_{t_{n-1}^{d }}^{t_{n}^{d}}\max\left(\left[\widehat{\mu}+\sum_{\{k,i,j\}}a_{2}^{i}\max\left( a_{1}^{i}(s-t_{k}^{j})+b_{1}^{i},0\right)\right],0\right)ds \tag{11}\]

Let \(\mathbf{X}\equiv\{x_{1},\ldots,x_{P^{*}}\}\) be the sequence of _zero-crossings_ obtained for all the \(P^{*}\) neurons, where the zero-crossing for the \(i\)th neuron is obtained by solving:

\[a_{1}^{i}(x_{i}-t_{k}^{j})+b_{1}^{i} = 0\] \[x_{i} = t_{k}^{j}-\frac{b_{1}^{i}}{a_{1}^{i}}.\]

\(P^{*}\) represents the total number of neurons used in the neural networks for modeling the kernels, \(\widehat{\phi}_{dj}\) with \(1\leq d,j\leq D,\) of a \(D\)-dimensional nonlinear Hawkes process. For example, if each kernel in a \(D\)-dimensional nonlinear Hawkes process is modelled using a network with \(P\) neurons, then \(P^{*}=P(D^{2})\).

We re-index the sequence \(\mathbf{X}\) in monotone increasing order as \(\mathbf{S}\equiv\{s_{1}\leq s_{2}\leq\cdots\leq s_{P^{*}}\}.\) Let \(\mathbf{S}^{*}\equiv\{s_{l}\leq\cdots\leq s_{u}\},\) where \(t_{n-1}^{d}\leq s_{l}\leq\cdots\leq s_{u}\leq t_{n}^{d},\) and \(1\leq l\leq u\leq P^{*},\) be the largest subsequence of the sorted sequence \(\mathbf{S}.\) Therefore, \(\mathbf{S}^{*}\) is the set of all zero-crossings that lie in the range \([t_{n-1}^{d},t_{n}^{d}].\)
We then write,

\[\int_{t_{n-1}^{d}}^{t_{n}^{d}}\widehat{\lambda}_{d}^{*}(s)ds=\int_{t_{n-1}^{d}}^{s _{1}}\max\left(\widehat{\lambda}_{d}(s),0\right)ds+\ldots\int_{s_{q-1}}^{s_{q}} \max\left(\widehat{\lambda}_{d}(s),0\right)ds\ldots+\int_{s_{u}}^{t_{n}^{d}} \max\left(\widehat{\lambda}_{d}(s),0\right)ds,\]

and exploit the fact that \(\widehat{\lambda}_{d}(s)\) will be linear in \(s\) in the sub-intervals, \([t_{n-1}^{d},s_{l}],\ldots,[s_{u},t_{n}^{d}]\). If \(\widehat{\lambda}_{d}(s)\) is linear in \(s\) in the interval \([s_{q-1},s_{q}]\), then there will be at the max one zero crossing for the function \(\max(\widehat{\lambda}_{d}(s),0)\) that lies within this interval. If \(s_{q}^{\text{out}}\) is the zero-crossing for \(\max(\widehat{\lambda}_{d}(s),0)\) and \(s_{q}\leq s_{q}^{\text{out}}\leq s_{q}\), we can write:

\[\int_{s_{q-1}}^{s_{q}}\max\left(\widehat{\lambda}_{d}(s),0\right)ds=\left( \int_{s_{q-1}}^{s_{q}^{\text{out}}}\widehat{\lambda}_{d}(s)ds\right)1_{ \widehat{\lambda}_{d}(s_{q-1})>0}+\left(\int_{s_{q}^{\text{out}}}^{s_{q}} \widehat{\lambda}_{d}(s)ds\right)1_{\widehat{\lambda}_{d}(s_{q})>0}.\]

In case \(s_{q}^{\text{out}}\notin[s_{q-1},s_{q}]\) the above can be integrated as:

\[\int_{s_{q-1}}^{s_{q}}\max\left(\widehat{\lambda}_{d}(s),0\right)ds=\left( \int_{s_{q-1}}^{s_{q}}\widehat{\lambda}_{d}(s)ds\right)1_{\widehat{\lambda}_{ d}(s_{q-1})>0}.\]

Therefore, we evaluate \(\widehat{\Lambda}_{d}(t_{n}^{d})\), by splitting the integral into intervals within which \(\max(\widehat{\lambda}_{d}(s),0)\) is linear in \(s\). The integral,

\[\int_{t_{n-1}^{d}}^{t_{n}^{d}}\widehat{\lambda}_{d}^{*}(s)ds\]

as well as its gradients with respect to the network parameters can then be exactly evaluated. The expression of the above integral in a sub-interval \([s_{q-1},s_{q}],\) where \(\widehat{\lambda}_{d}^{*}(s)\) is linear is provided in the Appendix B.

## 4 Experiments and Results

In this section, we evaluate the effectiveness of the NNNH method on synthetic and real-world datasets1.. We assess the performance of the NNNH method in estimating the conditional intensity function of a non-homogeneous Poisson Processes (NHPP). We then showcase the adaptability of the NNNH method by utilizing one-dimensional non-linear Hawkes processes. Specifically, we explore the NNNH's capacity to estimate Hawkes processes with smooth, non-smooth, and negative kernels, and examine its robustness in handling different variants of non-linear Hawkes models. Additionally, we investigate the NNNH's ability to estimate kernels and base intensity function for multidimensional non-linear Hawkes processes, where we evaluate a simultaneous combination of smooth, non-smooth, and inhibitive kernels.

Footnote 1: The source code and the dataset used in the experiments described here are available at [https://github.com/sobin-joseph/NNNH](https://github.com/sobin-joseph/NNNH)

We investigate the practical applications of non-linear Hawkes models through two case studies. In the first case, we apply the NNNH method to analyze tick-by-tick cryptocurrency trading data on the Binance exchange, seeking to identify any causal relationships between the buy and sell market orders of the BTC/USD and the ETH/USD pairs. For the second case, we utilize the NNNH to analyze neuronal spike trains recorded from the motor cortex of a monkey, revealing the interdependence between individual neurons and how they function in tandem.

Prior to delving into the outcomes of our numerical experiments, we provide a brief overview of the data preprocessing steps and the hyper-parameters selected for fitting the NNNH model to a dataset. Specifically we discuss the scaling approach adopted for the dataset, the choice of the number of neurons opted for the network, the initial learning rates employed for the Adam optimizer, and the stopping criteria for the optimizer.
### Data preprocessing and choice of hyper-parameters:

For all our experiments, we divide the dataset into training, validation, and test set. The dataset is partitioned in a 60:20:20 ratio for analysis. Due to the dependence on history in the Hawkes process, we do not alter the chronology of the events. Before splitting the dataset, we first scale the timestamps. Scaling the dataset is critical step in the preprocessing stage of building neural networks as it can enhance their performance, stability, and convergence. Additionally, scaling helps to ensure that all input features have a similar range of variance, which can lead to improved initialization of the network parameters and overall performance during training. Given dataset \(\mathcal{S}=(t_{n},d_{n}),\) where \(t_{n}\in[0,T),\) let \(\{t_{1}^{d},\ldots,t_{k}^{d},\ldots,t_{N_{d}(T)}^{d}\},\) be the ordered arrivals for the dimensions \(d=1,\ldots,D.\) We define \(T_{\max}\) as

\[T_{\max}=\max\left(t_{N_{1}(T)}^{1},\ldots,t_{N_{d}(T)}^{d},\ldots,t_{N_{D}(T)} ^{D}\right),\]

and \(N(T_{\max})\) as:

\[N(T_{\max})=\sum_{d=1}^{D}N_{d}(T_{\max}).\]

We scale the original timestamps as:

\[\widehat{t}_{n}^{d}=t_{n}^{d}\frac{N(T_{\max})}{T_{\max}}.\]

Initialization of network parameters and choice of batch sizeIn the NNNH method as discussed in Section 3 we model each kernel and base intensity function as a feed-forward neural network. To model the base intensity functions, \(\widehat{\mu}_{d}(t),\) we use a feed-forward network with fifty neurons in the hidden layer. The following initialization is used for the parameters of the network,

\[a_{1}^{i} \sim \mathcal{U}(-10^{-3},10^{-3})\] \[a_{2}^{i} \sim \mathcal{U}(0,0.2)\] \[b_{1}^{i} \sim \mathcal{U}(-1,1),\]

where \(\mathcal{U}(a,b),\) denotes uniform distribution between \(a\) and \(b\).

For the kernels, \(\widehat{\phi}_{dj}(t),\)\(1\leq d,j\leq D,\) of the nonlinear Hawkes process we use feed-forward networks with thirty two neurons in the hidden layer and the following scheme for initialising the network parameters,

\[a_{1}^{i} \sim \mathcal{U}(0,-0.3)\] \[a_{2}^{i} \sim \mathcal{U}(0,0.2)\] \[b_{1}^{i} \sim \mathcal{U}(0,0.3).\]

If we use a constant baseline intensity we initialize \(\widehat{\mu}=1.\)

We use a batch size of hundred for calculating the stochastic gradient. We use different learning rates for updating the parameters in the hidden layer and the output layer, as we observe faster convergence with this choice. The following learning rates were used for all the experiments: For fitting the networks used to model the base intensity function we used a learning rate of \(10^{-3}\) for updating the network parameters in the output layer and \(10^{-6}\) for updating the parameters in the hidden layer. For the networks used to model the kernel function the we used corresponding learning rates of \(10^{-2}\) and \(10^{-3}\) for the output and hidden layers, respectively.

Stopping criteria:In this study, for all the experiments, we utilize a stopping criterion based on the negative log-likelihood value computed on the validation dataset. The negative log-likelihood values are computed at the end of each iteration of the batch stochastic gradient descent algorithm. The NNNH parameter estimation algorithm stops when the updated parameters fail to improve the best
recorded validation error for a specified number of iterations, following the approach proposed by Goodfellow et al. (2016). This early stopping criterion helps prevent over fitting, as demonstrated in Figure 3, where we observe a reduction in training negative log-likelihood value while the validation negative log-likelihood value starts increasing after a certain number of iterations, indicating that the model is starting to over-fit the training data.

Through numerical experiments, Section 4.5 presents our investigation of the sensitivity of NNNH estimation to various hyper-parameter choices.

### Estimation of Non Homogenous Poisson Process Intensity

A Non-Homogenous poisson Process (NHPP) \(N(t),t\geq 0\) is a generalized counting process in which the rate at which events occurs varies with time. The intensity function, denoted by \(\mu(t):\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}\), captures the varying rate of the events over time, which could be influenced by factors such as external events, seasonal patterns, or other underlying phenomena. To comprehend this counting process, estimating the intensity function is necessary. The intensity function of NHPP can be estimated parametrically or non-parametrically. Parametric estimation assumes that the parametric form of the underlying intensity function is known. For instance, Rigdon and Basu (1989) estimates the parameters by assuming the intensity function as power law function, while Lee et al. (1991) uses a general exponential-polynomial-trigonometric intensity function. If one does not assume the parameteric form of the intensity function, they would resort to non-parametric techniques to estimate the conditional intensity function. In Leemis (1991), a technique is proposed that utilizes a piecewise linear function to estimate the cumulative intensity of the NHPP, while Xiao and Dohi (2013) presents another non-parametric approach that employs a wavelet-based estimation method to estimate the intensity function.

In this section, we show that the NNNH technique can be used to model the intensity function, \(\mu(t),\) of a non-homogeneous Poisson process. A single-layered feed-forward neural network can be used to model the intensity function as,

\[\hat{\mu}(t)=\max\left(b_{2}+\sum_{i=1}^{p}a_{2}^{i}\max(a_{1}^{i}t+b_{1}^{i},0 ),0\right), \tag{13}\]

Figure 3: Plot of Negative log-likelihood for train and validation datasets in one dimensional case, indicating the stopping criteria.
where, as explained in Section3, \([a_{1}^{i},b_{1}^{i},a_{2}^{i},b_{2}],\,1\leq i\leq p,\) are the parameters of the network that are estimated by maximizing the log-likelihood value,

\[\mathcal{L}(\mathbf{\Theta})=\sum_{t_{n}\in\mathcal{S}}log(\widehat{\mu}(t_{n}))- \int_{0}^{T}\widehat{\mu}(s)ds, \tag{14}\]

for the training dataset \(\mathcal{S}\). We use the SGD with Adam, for adaptive learning rates, to find the optimal network parameters, as discussed in Section 3.1.

We consider an exponential, linear, polynomial, and a trigonometric function, as given in Table 1, to model the conditional intensity. The arrival times are simulated using the thinning algorithm of Lewis and Shedler (1979). Table 1 compares the true negative log-likelihood of the simulated NHPP with the negative log-likelihood values obtained using the fitted NNNH model 2.

Footnote 2: The source code and the dataset used in the experiments described here are available at [https://github.com/sobin-joseph/NNNH](https://github.com/sobin-joseph/NNNH)

We see that for all the cases considered, the estimated negative log-likelihood values are reasonably close to the true negative log-likelihood values. Figure 4 compares the NNNH estimates of the intensity function with the true intensity function. We see for all the cases considered the NNNH is able to recover the intensity function reasonably well.

### Univariate Hawkes Estimation

This section focuses on examining the effectiveness of the NNNH method in estimating one-dimensional non-linear Hawkes processes using the following criteria:

* Estimation of non-smooth kernels
* Estimation of negative kernels
* Estimation of non-linear Hawkes processes for different variations of \(\Psi\).
* Estimation of Hawkes processes with varying base intensity function.

We simulate the arrival times for the different variants of the Hawkes processes using the _Ogata's thinning algorithm_ proposed in Ogata (1981). We compare the performance of our algorithm to the following state of the art methods:

* WH: An algorithm proposed in Bacry and Muzy (2014), a non-parametric estimation method which solves a Wiener-Hopf system derived from the auto-covariance of the multivariate Hawkes processes.
* Bonnet: A maximum likelihood based estimation method for Hawkes processes with self excitation or inhibition as proposed in Bonnet et al. (2021).

Estimation of non-smooth kernels:We first consider a linear Hawkes process with a non-smooth kernel, specifically a rectangular kernel of the form,

\[\text{Rectangular, }\phi(t)=\left\{\begin{array}{l}\alpha\beta,\text{if}\, \delta\leq t\leq\delta+\frac{1}{\beta},\\ 0,\text{otherwise}\end{array}\right.\]

\begin{table}
\begin{tabular}{|c|l|c|c|} \hline Function & Underlying equation & TNLL & NNLL \\ \hline \hline Exponential & \(\mu(t)=0.5e^{-0.001t}\) & 1177 & 1173 \\ \hline Linear & \(\mu(t)=0.5\) & 2528 & 2518 \\ \hline Parabola & \(\mu(t)=2\times 10^{-7}(1.5t-2000)\) & 1703 & 1699 \\ \hline Sin Curve & \(\mu(t)=0.4(\text{sin}(2\pi t\times 0.0004-1000)+1.1)\) & 3482 & 3538 \\ \hline \end{tabular}
\end{table}
Table 1: The NHPP functions and their corresponding True Negative Log-Likelihood(TNLL) and NNNH Negative Log-Likelihood(NNLL) obtained from the NNNH method
The smoothness assumption of the kernel is a prerequisite for some non-parametric methods used in estimating the kernels of Hawkes processes, such as the Markovian Estimation of Mutually Interacting Process (MEMIP) proposed in Lemonnier and Vayatis (2014). As a result, it is necessary to compare the performance of the NNNH method using a non-smooth kernel in order to assess its effectiveness. We simulate the process using the following parameter values for the rectangular kernel,

\[\mu=0.05,\,\alpha=0.7,\,\beta=0.4,\,\delta=1.\]

Simulating the process until \(T=60000\) yields \(N(T)=10001\) arrivals. As a first step, we visually compare the estimated kernels obtained using the NNNH and the benchmark methods with the true kernel. Figure 5 illustrates the kernels fitted by the three estimation methods. As expected we see that the non-parametric models, i.e., the NNNH and the WH are able to better recover the ground truth. This is also reflected in the negative log-likelihood values obtained for WH, Bonnet, and NNNH methods, i.e., \(4228\), \(12838\), and \(2800\) respectively.

Figure 4: NNNH estimated kernel and base intensity for a one-Dimensional Hawkes process with sin base intensity and exponential kernel
Estimation of negative kernels:We next consider a non-linear Hawkes process,

\[\lambda^{*}(t)=\max\left(\mu+\sum_{\tau<t_{n}}\phi(t-\tau),0\right),\]

where, \(\Psi\) in Equation 3 is a \(\max\) function, i.e. \(\Psi(\lambda(t))=\max(\lambda(t),0)\). The inhibitive kernel, \(\phi,\) is specified as an exponential function given by,

\[\phi(t)=-0.5e^{-2t}.\]

Given that non-parametric methods, such as the Expectation Maximization (EM) approach introduced in Lewis and Mohler (2011), are capable of estimating only positive kernels, it is of interest to evaluate the capacity of the NNNH method in estimating the negative kernels.

We simulate the above process for \(T=14000\), with a base intensity, \(\mu=0.9\). This resulted in \(N(T)=10001\) events. We first compare the kernel obtained using the NNNH, WH, and Bonnet method with the ground truth in Figure 6. All three methods are able to recover the kernel, but visually the NNNH and Bonnet methods appear to provide better estimates. Furthermore, this conclusion is supported by the estimated negative log-likelihood values, which are \(9833\), \(9559\), and \(9564\) for the WH, Bonnet, and the NNNH, respectively.

Estimation of non-linear Hawkes processes for different variations of \(\Psi:\) Diverse variations of the non-linear Hawkes process can be derived by making distinct selections of the function \(\Psi\) in

Figure 5: Fitted Rectangular kernel: A comparison of the NNNH, WH, and Bonnet estimates with the actual (real) kernel.

Figure 6: Fitted negative exponential kernel: A comparison of teh NNNH, WH, and Bonnet estimates with the actual kernel.
Equation 3. In contrast to the preceding example where \(\Psi\) was chosen as a \(\max\) function, in this instance, we adopt a sigmoid function for \(\Psi\) given by

\[\lambda^{*}(t)=\frac{1}{1+e^{-(\lambda(t)-2)}}\]

. We consider the following kernel,

\[\phi(t)=e^{-2t},\]

for this example.

The nonlinear Hawkes process was simulated for a duration of \(T=10000\), resulting in \(N(T)=3028\) events. According to Equation 9, the NNNH model converts a non-linear Hawkes process with any \(\Psi\) into a non-linear Hawkes process with \(\Psi\) as a \(\max\) function. As a result, the recovered kernels from the NNNH method might not correspond to the kernels of the original process. However, we can compare the actual simulated intensity process with the intensities recovered from the NNNH method and the WH method. Figure 7 demonstrates that the recovered intensities match the simulated intensities.

The accuracy of NNNH estimates of quantiles is evident from the QQ plot obtained from the fitted WH and the NNNH method on the test dataset, as shown in Figure 8. The negative log-likelihood values for NNNH and WH are \(2954\), and \(3813\) respectively.

Figure 8: QQ plot comparing estimated and true intensities for sigmoid nonlinear Hawkes process using WH and NNNH

Figure 7: Comparison of real and estimated intensity function by WH and NNNH for nonlinear Hawkes process with sigmoid as \(\psi\) function.
Estimation of Hawkes processes with varying base intensity functionTo complete the univariate case test for the NNNH, we examine a Hawkes process with a base intensity that varies over time. We adopt a trigonometric sin function for the time-varying base intensity \(\mu(t)\), which is specified in Table 1. The associated kernel for the Hawkes process is expressed as:

\[\phi(t)=e^{-2t}.\]

We simulate the process till \(T=3000,\) which results in \(N(T)=2482\) events. The comparison between the recovered base intensity function and the excitation kernel using the NNNH method with the true values is depicted in Figure 9. The results demonstrate that the NNNH approach can accurately estimate both the base intensity function and the kernel concurrently.

### Multivariate Hawkes Estimation

In this section, we evaluate the NNNH approach for estimating multivariate non-linear Hawkes processes using a simulated dataset. The multivariate non-linear Hawkes process is simulated using the _Ogata's thinning algorithm_. The parameters of the neural networks utilized to model the kernels in the NNNH method are optimized by maximizing the log-likelihood through the stochastic gradient descent method, as described in Section 3.1. We utilize the WH method (Bacry and Muzy, 2014) and the Bonnet Multivariate method (Bonnet et al., 2022) as our comparative models. While the WH method is non-parametric, the Bonnet Multivariate approach assumes a parametric structure for the kernels in the Hawkes process. Both reference models are capable of estimating kernels that exhibit inhibitory effects.

We consider two examples for the simulated case. We first consider a multivariate Hawkes process, with both positive and negative exponential kernels. Specifically the parameter chosen for the model are:

\[\alpha=\begin{bmatrix}-0.9&3\\ 1.2&1.5\end{bmatrix},\]

\[\beta=\begin{bmatrix}5&5\\ 8&8\end{bmatrix}\]

, and a constant base intensity

\[\mu=\begin{bmatrix}0.5\\ 1.0\end{bmatrix}.\]

The kernels of the process are defined as: \(\phi_{dj}(t)=\alpha_{dj}e^{-\beta_{dj}t},\)\(1\leq d,j\leq 2.\)

We simulate the process till \(T=1000,\) which results in \(N(T)=1002\) events. Figure 10 plots the true kernels and the kernels recovered using the NNNH, the WH, and the Bonnet Multivariate

Figure 9: Estimated kernel and base intensity by applying NNNH to a one-dimensional Hawkes process with sin base intensity and an exponential kernel
methods. Through a visual examination, it can be observed that the Bonnet Multivariate approach closely approximates the actual kernels. This outcome is in line with expectations as the simulated data employed a parametric form of the Hawkes process that is specifically suitable for the Bonnet Multivariate model. Regarding the non-parametric models, the WH estimates exhibit greater variability compared to the NNNH estimates. The respective negative log-likelihoods for the NNNH, the WH, and the Bonnet Multivariate methods are \(1480\), \(1967\), and \(1460\).

Subsequently, we examine a bivariate non-linear Hawkes process that encompasses a combination of diverse types of kernels. Precisely, we select the kernels as follows:

\[\text{exponential kernel}\,\phi_{11}(t)=0.3e^{-3t},\] \[\text{rectangular kernel}\,\phi_{12}(t)=\left\{\begin{array}{l}0.7\times 0.4,\text{if}\,1\leq t\leq 1+\frac{1}{0.4},\\ 0,\text{otherwise}\end{array}\right.,\] \[\text{negative exponential}\,\phi_{21}(t)=-0.2e^{-t},\] \[\text{exponential}\,\phi_{22}(t)=0.4e^{-2t},\]

Figure 10: The kernels fitted using the NNNH, WH, and Bonnet Multivariate methods, for a bivariate Hawkes process with both positive and negative kernels.
and use the following base intensity values:

\[\mu=\begin{bmatrix}0.1\\ 0.2\end{bmatrix}\]

We perform simulations up to time \(T=10000\), resulting in the occurrence of a total \(8489\) events in the two dimensions. Figure 11 presents the findings obtained from the NNNH approach in comparison with other techniques. Notably, the advantages of employing a non-parametric estimation method for Hawkes processes are evident in this context. Both the NNNH and the WH methods successfully capture all the kernels, whereas the Bonnet Multivariate approach does not accurately depict the rectangular kernel. It is worth mentioning that the estimated kernels obtained using the WH method appear to be more erratic than those of the NNNH. The corresponding negative log-likelihood values for the WH, the Bonnet Multivariate, and the NNNH are \(3051\), \(3062\), and \(2899\) respectively.

### Sensitivity Analysis

In the context of the NNNH approach, we have made deliberate decisions concerning the quantity of neurons deployed in the network's hidden layer to emulate the kernels and the learning rate employed to regulate the parameter updates in the Adam algorithm. To gauge the sensitivity of the NNNH method, we have conducted a numerical investigation with regard to these decisions. Specifically, we have examined the influence of variations in the number of neurons and learning rates on the negative log-likelihood value, while holding all other factors constant. Note that the likelihood values reported are for the validation dataset which was not used in the training of the network

The sensitivity analysis is conducted for the negative exponential kernel, which is discussed in Section 4.3. The study focuses on the variation in the number of neurons employed in the neural networks, and the corresponding impact on the negative log-likelihood values. The results are presented in

Figure 11: The fitted kernels, using the NNNH, WH, and Bonnet Multivariate method, for a bivariate non-linear Hawkes process with a mix of exponential and rectangular kernels.
Table 2 presents the percentage variation in the negative log-likelihood values linked to the optimal parameters obtained by varying the learning rate for the output layer of the neural network. The learning rate for the hidden layer remains unchanged, and is taken as ten percent of the learning rate for the output layer.

We find that the optimal network parameters obtained using varying learning rates results in similar negative log-likelihood values and therefore use a learning rate of 0.01 for all our experiments.

### Real Data

#### 4.6.1 Financial Dataset

In this study, we evaluate the efficacy of the NNNH method on high-frequency order book data pertaining to two of the most frequently traded cryptocurrencies, namely bitcoin and Ethereum. The data comprises of buy and sell trade records for the BTC-USD (bitcoin- US dollar) and ETH-USD (Ethereum- US Dollar) pairs. We streamed the Binance exchange order book data, as several popular cryptocurrencies are traded in this exchange, and the exchange has high trade volumes. We obtain the tick-by-tick arrival times for both buy and sell trades from the exchange. The market is composed of makers and takers, with makers generating buy or sell orders that are not carried out immediately, thereby creating liquidity for that cryptocurrency. In contrast, takers place market orders that are executed instantly, thereby taking away the liquidity.

The Binance exchange furnishes two streams of tick-by-tick data: trade arrival data and trade stream data. The trade arrival data comprises of limit orders, which are orders placed with a specified price limit, and market orders, which are orders executed at the prevailing market price. Limit orders are executed when the best available market price reaches the set limit, while market orders are executed instantly at the current best limit order. The Binance trade-stream data provides the timestamp of these order arrivals, along with price and volume features, and a unique identifier for the buyer/seller. Since a single market order may necessitate multiple limit orders to fulfill the requested volume, several trades are recorded with a common identifier. Therefore, we cleaned the dataset by filtering out the data with common IDs and retaining only the unique trade events. Finally, based on whether the buyer was a market maker or taker, the trades were marked as either a buy or a sell market order.

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline No of Neurons & NLL & \(\%\) change & \(\%\) time change \\ \hline
2 & 2057.78 & -0.061 & +164.1 \\
4 & 2058.01 & -0.057 & +95.4 \\
8 & 2059.06 & -0.004 & +89.9 \\
16 & 2059.14 & -0.001 & +90.9 \\
32 & 2059.16 & 0.00 & 0.00 \\
64 & 2061.29 & +0.10 & +6.7 \\
128 & 2061.63 & +0.12 & +77.3 \\
256 & 2061.62 & +0.12 & +76.1 \\ \hline \end{tabular}
\end{table}
Table 2: Percentage change in negative log-likelihood (NLL) values and corresponding computational time for different numbers of neurons in the neural network

\begin{table}
\begin{tabular}{|c|c|c|} \hline Learning Rate & NLL & \(\%\) change \\ \hline
0.01 & 2064.54 & +0.00 \\
0.005 & 2064.19 & -0.02 \\
0.001 & 2065.49 & +0.05 \\
0.0005 & 2067.87 & +0.12 \\ \hline \end{tabular}
\end{table}
Table 3: Percentage change in negative log-likelihood (NLL) values for different learning rates.
Our analysis focuses exclusively on the trades conducted for the BTC-USD and ETH-USD pairs, which account for the predominant trading volumes in the cryptocurrency exchange. The arrival time for market orders for these two pairs can be modelled as a four dimensional non-linear Hawkes process, i.e.,

First Dimension: Intensity process for the sell market orders for the BTC-USD pair

Second Dimension: Intensity process for the buy market orders for the BTC-USD pair

Third Dimension: Intensity process for the sell market orders for the ETH-USD pair

Fourth Dimension: Intensity process for the buy market orders for the ETH-USD pair.

The analysis was performed for the market orders arrival times on December 7 2021, between 12:00 to 12:10 (UTC) on the Binance exchange. Table 4 summarises the data after classifying the trades as buy or sell market order events.

Our aim is to employ the NNNH method to jointly model the intensity process of buy and sell market orders for BTC-USD and ETH-USD currency pairs. This modelling approach enables us to gain insights into the causal relationships between the two pairs, such as whether the arrival of a buy BTC-USD order affects the arrival of a sell ETH-USD market order. Furthermore, identifying the functional form of the self and cross modulation due to the arrival of market orders is of interest. This highlights the necessity of non-parametric estimation techniques, as the true form of the modulation function remains unknown. To achieve our objective, we partition the dataset into train, validation,

\begin{table}
\begin{tabular}{|c|c|c|c|} \hline Market Order Type & BTC-USD & ETH-USD & Aggregate \\ \hline Sell & 4219 & 3480 & 7699 \\ Buy & 4374 & 2999 & 7373 \\ \hline Total & 8593 & 6479 & 15072 \\ \hline \end{tabular}
\end{table}
Table 4: Summary of the crypto-currency market orders recorded on December 7 2021 between 12:00 and 12:10 UTC

Figure 12: Histogram for the inter event time arrivals of the BTC-USD sell trade
and test sets. We optimize the parameters of the model using stochastic gradients computed on the train set and use the validation set for applying the early stopping criteria. Specifically, if there is no improvement in negative log-likelihood values computed on the validation set for ten consecutive iterations, we stop training. Finally, we evaluate the goodness of fit using the negative log-likelihood reported on the test set.

We adopt the widely used non-parametric methods, namely the WH and the EM method (Lewis and Mohler, 2011), as reference models. In Figure 13, we compare the obtained kernels using the three methods. The corresponding negative log-likelihood values are \(9630\), \(6732\), and \(6400\) for the EM, WH, and NNNH methods, respectively. The estimated base intensities for BTC-USD sell, BTC-USD buy, ETH-USD sell, and ETH-USD buy, in the order of appearance, using the NNNH method are \(0.0028\), \(0.0032\), \(0.0021\), and \(0.0022\), respectively.

By examining Figure 13, we can draw the following observations for the market orders that arrived during the training window:

* Significant self-excitation is observed for both BTC-USD and ETH-USD buy and sell market orders.
* BTC-USD sell orders result in excitation of ETH-USD buy and sell orders.
* The self-excitation of ETH-USD buy and sell orders is higher compared to their respective BTC-USD orders.
* A certain level of inhibition is observed in BTC-USD buy orders due to sell BTC-USD orders.

Notably, all three methods effectively capture a lagged cross-excitation effect in ETH-USD buy orders caused by sell BTC-USD orders.

Figure 13: Fitted kernels for BTC-USD and ETH-USD sell-buy market orders estimated using the NNNH, the WH, and the EM methods.
Predictive capacity of the NNNH methodAccurately inferring the kernels is critical in developing an improved predictive model for the arrival process. For example, a trading algorithm relies on predicting the time at which the next market order will arrive with 90% confidence. With an accurate prediction model, approximately 90% of the predictions are expected to be correct, and in roughly 10% of the cases, the order will arrive after the predicted time. If the model consistently produces accurate predictions (i.e., does not fail for 10% of the cases), then the predictions are overly conservative, indicating that the predicted time is set too far in the future. Conversely, if more than 10% of the predictions fail, the predicted time is closer than anticipated. Therefore, a prediction model that produces predictions with Q% certainty should have Q% correct outcomes. The accuracy of such a prediction model can be evaluated using a QQ plot.

Figure 14 presents the QQ plot for the BTC-USD and ETH-USD order arrivals based on the EM, WH, and NNNH models. The sample quantiles for the QQ plot are obtained from the test dataset. The results indicate that all three models make reasonably accurate predictions. Notably, the NNNH method appears to have more precise predictions, particularly for BTC-USD and ETH-USD sell orders, as evidenced by the visual assessment.

#### 4.6.2 Neuron Spike Train Dataset

In this study, we utilized the NNNH method to analyze a dataset of neuron spikes, which was obtained from an experiment conducted by Engelhard et al. (2013). The primary objective of the experiment was to examine the relationship between unit recordings from the motor cortex of monkeys and the position and grip strength of their hands as they utilized a joystick to manipulate a robotic arm. Consistent with the approach taken by Aljadeff et al. (2016), we analyzed the data obtained from the nine simultaneous electrode recordings. Although the dataset included additional information

Figure 14: QQ plot comparing order arrivals of BTC-USD and ETH-USD using EM, WH, and NNNH methods
regarding hand motion, such as cursor trajectory and grip force, our analysis was focused exclusively on the neuron spike train.

To provide a brief biological background about neuron spikes, we refer to Duval et al. (2022). Neurons employ an electrical wave of short duration, known as the action potential or spike, to reliably transmit signals from one end to the other (Castelfranco and Hartline, 2016). Action potentials are all-or-none events and are triggered when the neuron's membrane potential, which is the electrical potential difference between the inside and the outside of the neuron, is sufficiently large. Between two successive action potentials, the neuron adds up its inputs, causing changes in its membrane potential. When this potential reaches a sharp threshold, the neuron fires a spike that propagates (Luo et al., 2015).

The histogram in Figure 15 presents the inter-arrival times of spikes generated by a single neuron. Notably, the peak of the histogram is marginally displaced from the origin, in contrast to the histogram of a conventional exponential distribution. This observation suggests the presence of an inhibitory process. To model the neuron spike train, we employed a non-linear Hawkes process that is nine-dimensional. Each dimension of the process is associated with the intensity of neuron spike arrivals at an individual electrode. We use the NNNH method to estimate the kernel functions and the base intensities for this process.

Table5 gives the number of spike events for each neurons used for the analysis.

Figure 16 displays the kernels that were fitted using the NNNH method to model the neuron spike train. Additionally, Table 6 reports the base intensity values for the arrival of neuron spikes at each electrode. The results indicate that a majority of the neurons exhibit a self-inhibitory characteristic. Specifically, after firing, neurons typically undergo a refractory period during which they cannot generate additional spikes. The inhibitory pattern evident in the plots likely corresponds to this refractory period. Based on the figures, it is reasonable to assume that the refractory period for the recorded motor neurons falls within the 20-60 ms range.

Figure 15: Histogram of inter-spike time intervals for neuron 1 in grip-and-reach task

\begin{table}
\begin{tabular}{|c|c|c|c|c||c|c|c|c|} \hline Neuron & \(1\) & \(2\) & \(3\) & \(4\) & \(5\) & \(6\) & \(7\) & \(8\) & \(9\) \\ \hline \hline No. of events & 4121 & 1895 & 912 & 2509 & 221 & 2784 & 149 & 3854 & 1360 \\ \hline \end{tabular}
\end{table}
Table 5: Summary statistics of neuron spike data
## 5 Conclusion

This paper proposes a non-parametric method for estimating non-linear Hakes processes. We call the method Neural Network for Non-linear Hawkes processes or the NNNH method. The NNNH method involves modelling the kernels and the base intensity functions of a non-linear Hawkes process as individual neural networks. The parameters of the neural networks involved are jointly estimated by maximizing the log-likelihood function using the batch stochastic gradient descent with Adam for adaptive learning. To apply the SGD method, an unbiased estimator for the gradients with respect to the network parameters is proposed. The paper provides an efficient scheme to evaluate the integrated intensity and its gradients, which can be a computational bottleneck.

\begin{table}
\begin{tabular}{|c|c|c|c||c|c|c|c|c|} \hline \(\mu_{1}\) & \(\mu_{2}\) & \(\mu_{3}\) & \(\mu_{4}\) & \(\mu_{5}\) & \(\mu_{6}\) & \(\mu_{7}\) & \(\mu_{8}\) & \(\mu_{9}\) \\ \hline \hline
0.009 & 0.035 & 0.001 & 0.046 & 0.001 & 0.043 & 0.001 & 0.076 & 0.022 \\ \hline \end{tabular}
\end{table}
Table 6: Base intensities of neuron spike data obtained using the NNNH for the grip-and-reach task

Figure 16: Kernels of motor neurons demonstrating the self-dependency and mutual dependency of neurons for the grip-and-reach task estimated by NNNH.
The effectiveness of the NNNH method in learning the underlying process is investigated through numerical experiments. Unlike many recent neural network-based models that predict conditional intensities, the NNNH model retains the interpretability of the Hawkes process by enabling the inference of the kernels rather than just direct prediction of the conditional intensities. The ability to recover the kernels is desirable for understanding causal relationships between the arrivals in different dimensions, for instance. The NNNH method performed satisfactorily for the diverse set of problems considered. We also see that the NNNH can be used to estimate non-homogeneous Poisson process. The NNNH method is able to recover non-smooth kernels, kernels with negative intensities, and estimate other variants of the non-linear Hawkes processes.

The NNNH approach is applied to examine the process of order arrivals for buy and sell transactions on the Binance exchange in relation to the BTC-USD and ETH-USD currency pairs. We find evidence of both self and cross excitation in the order arrivals for the two currency pairs. Notably, self excitation is observed to be significant for both currency pairs, while the cross dependence exhibits an asymmetry whereby ETH-USD order arrivals are more strongly influenced by BTC-USD order arrivals than vice versa.

We applied the NNNH model to a publicly accessible dataset of neuron spikes obtained from the motor cortex of monkeys engaged in a grip-and-reach task. As anticipated, the kernels obtained from the analysis exhibit self-inhibition, which could be associated with the refractory period of the neurons. Our findings provide evidence of the NNNH method's efficacy in analyzing high-dimensional datasets.

The NNNH method has a disadvantage in that the estimation of model parameters relies on the use of stochastic gradient descent (SGD), which updates the network parameters iteratively using first-order derivatives. Convergence of SGD-based methods can be slow, often requiring several iterations before the stopping criterion is met. However, the use of SGD makes the NNNH method well-suited for online learning, enabling the model to be trained on new data points that are continuously streaming in, as opposed to a static dataset. Moreover, the NNNH method is amenable to parallel computing of the gradients for batch SGD. A potential avenue for future research is to extend the NNNH method to model marked non-linear Hawkes processes.
