Make Hawkes Processes Explainable by Decomposing Self-Triggering Kernels

TL;DR

This paper introduces a novel framework for decomposing Hawkes Processes into interpretable kernels, enhancing understanding of event dependencies and improving prediction accuracy in real-world datasets.

Contribution

It presents the first multiplicative kernel composition method for Hawkes Processes, enabling automatic decomposition and better interpretability of event triggering mechanisms.

Findings

Outperforms existing methods in event prediction accuracy

Provides interpretable decomposition of self-triggering kernels

First to introduce multiplicative kernel composition for Hawkes Processes

Abstract

Hawkes Processes capture self-excitation and mutual-excitation between events when the arrival of an event makes future events more likely to happen. Identification of such temporal covariance can reveal the underlying structure to better predict future events. In this paper, we present a new framework to decompose discrete events with a composition of multiple self-triggering kernels. The composition scheme allows us to decompose empirical covariance densities into the sum or the product of base kernels which are easily interpretable. Here, we present the first multiplicative kernel composition methods for Hawkes Processes. We demonstrate that the new automatic kernel decomposition procedure outperforms the existing methods on the prediction of discrete events in real-world data.

Figures

Tables

Table 1. Table 1: Base kernels and their equations.

Type	Equation
Exponential (EXP($\alpha$,$\beta$))	$\alpha e^{-\beta t}$
Power-Law (PWL(K,c,p))	$\frac{K}{{\left(c+t\right)}^p}$, $\left(p>1\right)$
Pulse (SQR(B,L))	$B\left(u\left(t\right)-u\left(t-L\right)\right)$
Sinusoidal (SNS(A,$\omega$))	$Asin\left(\omega t\right)$, $t\in [0,\frac{\pi }{\omega }]$

Table 2. Table 2: Stationarity conditions of the base kernels: the expressions in the stationarity condition should be set to less than 1.

Type	Stationarity Condition
EXP($\alpha$,$\beta$)	$\alpha /\beta$
PWL(K,c,p)	$K{c}^{1-p}/(p-1$), $\left(p>1\right)$
SQR(B,L)	$BL$
SNS(A,$\omega$)	$2A/\omega$

Table 3. Table 3: Stationarity Condition for Multiplicative Combination of the four Base Kernels.

Base Kernel	Base Kernel	Condition
EXP($\alpha$,$\beta$)	EXP($\alpha$,$\beta$)	${\alpha }_1{\alpha }_2/\left({\beta }_1+{\beta }_2\right)$ (closed under multiplication)
EXP($\alpha$,$\beta$)	PWL(K,c,p)	$\alpha K{\beta }^{p-1}{e}^{\beta c}\Gamma (1-p,\beta c)$
EXP($\alpha$,$\beta$)	SQR(B,L)	$\left(\alpha B\left(1-e^{-\beta L}\right)\right)/\beta$
EXP($\alpha$,$\beta$)	SNS(A,$\omega$)	$\left(A\alpha \omega \left(1+e^{\frac{-\beta \pi }{\omega }}\right)\right)/\left({\omega }^2+{\beta }^2\right)$
PWL$(K_1,c_1,p_1)$	PWL$(K_2,c_2,p_2)$	$\le \left(K_1{K}_2\right)/\left((p_1+p_2-1)\min {(c_1,c_2)}^{\left(p_1+p_2-1\right)}\right)$ (upper bound)
PWL(K,c,p)	SQR(B,L)	$\left(KB\left(c^{-\left(p-1\right)}-{\left(c+L\right)}^{-\left(p-1\right)}\right)\right)/\left(p-1\right)$
PWL(K,c,p)	SNS(A,$\omega$)	$\le KA\left({\left(c+\frac{\pi }{\omega }\right)}^{1-p}-c^{1-p}\right)/\left(1-p\right)$ (upper bound)
SQR(B,L)	SQR(B,L)	$BL$
SQR(B,L)	SNS(A,$\omega$)	$2AB/\omega$
SNS(A,$\omega$)	SNS(A,$\omega$)	$\pi A/\left(2\omega \right)$

Table 4. Table 4: Average log-likelihood, over each of the five financial data sets, for the two levels of decomposition.

Dataset	l(K1)	l(K2)
Technology	-2097.0	-1894.6
Healthcare	-2677.7	-2446.2
Industrial	-2309.7	-2127.9
Services	-2368.9	-2218.7
Utilities	-2427.9	-2266.3

Table 5. Table 5: Aggregate comparison, among the gradient descent based HP model and the first- and second-level decompositions of the proposed algorithm, for each of the five financial data sets.

Dataset	# of Seq.	$l(\text{K1},\text{K2})>l\left(\text{EXP}\right)$	$l\left(\text{K2}\right)>l\left(\text{K1}\right)$
Technology	70	67.14%	98.57%
Healthcare	117	62.39%	92.31%
Industrial	53	64.15%	94.34%
Services	61	54.09%	85.25%
Utilities	48	77.08%	93.75%

Equations

$$\lambda (t)=\underset{h\to 0}{\lim }\frac{\mathbb{E}[N(t+h)-N(t)]}{h}$$

$$\underset{h\to 0}{\lim }\frac{\mathbb{E}[N(t+h)-N(t)|\mathcal{H}(t)]}{h}=\mu +{\int }_{-\infty }^t\varphi (t-u)dN(u),$$

$$||\varphi ||:={\int }_0^{\infty }\varphi (t)dt\le 1,$$

$${\nu }_{\tau }^{(h)}=\frac{1}{h}E\left(({\int }_0^{h}d{N}_s-\Lambda h)({\int }_{\tau }^{\tau +h}d{N}_s-\Lambda h)\right)$$

$${\nu }_{\tau ,\delta }^{(h)}=\frac{1}{T}\sum _{i=1}^{\lfloor T/\delta \rfloor }(d{N}_{i\delta }^{(h)}-d{N}_{(i-1)\delta }^{(h)})(d{N}_{i\delta +\tau }^{(h)}-d{N}_{(i-1)\delta +\tau }^{(h)}),$$

$${\hat{\psi }}_z=\sum _{n=1}^{+\infty }{\hat{\varphi }}_z^n=\frac{{\hat{\varphi }}_z}{(1-{\hat{\varphi }}_z)}.$$

$$(1+{\hat{\psi }}_z^{\star })\Lambda (1+{\hat{\psi }}_z^{\star }{)}^{†}=\frac{\widehat{{\nu }_z^{(h)}}}{{\widehat{g_z}}^{(h)}},$$

$${\hat{\varphi }}_{i\omega }=1-e^{-\log |1+{\hat{\psi }}_{i\omega }|+iH(\log |1+{\hat{\psi }}_{i\omega }|)},$$

$${\left[\text{mbox}EXP\right]}^{k_1}\times {\left[\text{mbox}PWL\right]}^{k_2}\times {\left[\text{mbox}SQR\right]}^{k_3}\times {\left[\text{mbox}SNS\right]}^{k_4}$$

$${\left[EXP\right]}^k=\underset{\text{k terms}}{\underbrace{{\alpha }_1{e}^{-{\beta }_{1}x}\ast {\alpha }_2{e}^{-{\beta }_{2}x}\ast ...\ast {\alpha }_k{e}^{-{\beta }_{k}x}}}.$$

$${\left[\text{mbox}EXP\right]}^{k_1}\times {\left[\text{mbox}PWL\right]}^{k_2}\times {\left[\text{mbox}SQR\right]}^{k_3}\times {\left[\text{mbox}SNS\right]}^{k_4}$$

$$0\le {\int }_0^{\infty }EXP({\alpha }_1,{\beta }_1)EXP({\alpha }_2,{\beta }_2)dx<1$$

$$0\le {\int }_0^{\infty }{\alpha }_1{e}^{{\alpha }_{1}x}{\alpha }_2{e}^{{\beta }_{2}x}dx<1$$

$${\int }_0^{\infty }{\alpha }_1{e}^{-{\beta }_{1}x}{\alpha }_2{e}^{-{\beta }_{2}x}dx$$

$$0\le {\int }_0^{\infty }EXP(\alpha ,\beta )PWL(K,c,p)dx<1$$

$$0\le {\int }_0^{\infty }\alpha e^{-\beta x}\frac{K}{(x+c{)}^p}dx<1$$

$${\int }_0^{\infty }\alpha e^{-\beta x}\frac{K}{(x+c{)}^p}dx$$

$$0\le {\int }_0^{\infty }EXP(\alpha ,\beta )SQR(B,L)dx<1$$

$$0\le {\int }_0^L\alpha B{e}^{-\beta x}dx<1$$

$${\int }_0^L\alpha B{e}^{-\beta x}dx={\left[\frac{\alpha B{e}^{-\beta x}}{\beta }\right]}_0^L=\frac{\alpha B(1-e^{-\beta L})}{\beta }$$

$$0\le {\int }_0^{\infty }EXP(\alpha ,\beta )SNS(A,\omega )dx<1$$

$$0\le {\int }_0^{\frac{\pi }{\omega }}A\alpha e^{-\beta x}sin(\omega x)dx<1$$

$${\int }_0^{\frac{\pi }{\omega }}A\alpha e^{-\beta x}sin(\omega x)dx$$

$$0\le {\int }_0^{\infty }PWL(K_1,c_1,p_1)PWL(K_2,c_2,p_2)dx<1$$

$$0\le {\int }_0^{\infty }\frac{K_1}{(x+c_1{)}^{p_1}}\frac{K_2}{(x+c_2{)}^{p_2}}dx<1$$

$${\int }_0^{\infty }\frac{K_1}{(x+c_1{)}^{p_1}}\frac{K_2}{(x+c_2{)}^{p_2}}dx$$

$$0\le {\int }_0^{\infty }PWL(K,c,p)SQR(B,L)dx<1$$

$$0\le {\int }_0^L\frac{KB}{(x+c{)}^p}dx<1$$

$${\int }_0^L\frac{KB}{(x+c{)}^p}dx={\left[\frac{KB}{(1-p)(x+c{)}^{(p-1)}}]\right]}_0^L=\frac{KB(c^{-(p-1)}-(c+L{)}^{-(p-1)})}{p-1}$$

$$0\le {\int }_0^{\infty }PWL(K,c,p)SNS(A,\omega )dx<1$$

Full text

Introduction

Hawkes Processes (HPs) (?) model self-exciting behavior, i.e., when the arrival of one event makes future events more likely to happen. This type of behavior has been observed in various domains, such as earthquakes, financial markets, web traffic patterns, crime rates (?; ?) and social media (?).

As an example, in high-frequency finance, buyers and sellers of stocks demonstrate herding behavior (?; ?). After the main earthquake, several aftershocks follow according to a time-clustered pattern (?). In web data, hyperlink proliferation across pages exhibit self- and mutual-excitation (?). In criminology, gang-related retaliatory crime patterns are grouped in time (?). In social media, the ‘infectiousness’ of posts can be shown to be modeled through a self-excitement and mutual-excitement assumption (?).

In HPs analysis, parametric kernels capture intra-domain typical behaviors: quick time-decaying exponential excitation in the case of finance and web data (?; ?); slower power-law decay in earthquake-related data (?); and periodicity-inducing sinusoidal kernel in TV-watching data (?).

When an appropriate kernel is selected, the kernel parameters are fitted to predict future events. The parameters may be fitted to the data through the gradient descent (GD) method over a likelihood function penalized by a regularization criterion (e.g., Akaike Information Criterion) on the number of parameters (?). Another method of kernel estimation is through the use of the power spectrum of the second order statistics of the process: covariance density and normalized covariance (?). These are well defined when the self-triggering function induces what is called stationary behaviour.

However, kernel selection in HPs analysis is a challenging problem, since an appropriate kernel should be manually selected in practice. In this paper, we present a kernel structure search algorithm for HPs. Given base kernels, our algorithm finds the best fitting one, considering composition (sum and product) of base kernels. For verifying the stationarity property of each composite kernel, we also derived analytical expressions for the stationarity conditions. To our best knowledge, our method is the first multi-type kernel composition framework for HPs.

The main steps of the automatic framework, which will be thoroughly explained in the following sections, are discretized kernel estimation and greedy search in kernel composition space.

Related Work

Automatic analysis frameworks for Gaussian Processes (GPs) are proposed in (?) and (?). However, due to fundamental distinctions between GPs and HPs (such as stationarity conditions and causality assumptions for the latter), the techniques proposed for GPs can not be extended to HPs in a straightforward manner.

(?) uses exponential kernels for modeling quick-decay in finance or web data. (?) models slow decay influence with power-law kernels in earthquake, while (?) performs power law modeling experiments with social media-related data. (?) uses sinusoidal kernels for modeling periodicity-inducing influence in TV watching-related data (IPTV) , in which watching one episode of a TV program makes the viewer more likely to watch further ones. Since these shows are usually broadcasted weekly, the TV-watching behavior is likely to demonstrate a weekly self-excitement. In addition, according to (?), homicide rates show a pronounced seasonal effect, peaking in the summer and tapering in the winter.

More recent works, such as the neural network-based Hawkes processes in (?; ?) and the time-dependent Hawkes process (TiDeH) (?), allow for learning very flexible Hawkes processes with highly complicated intensity functions, while depending on the size and the quality of data. In this work, however, we focus on the interpretability, or explainability, of said functions and their corresponding typical behaviours, which are core factors on the Hawkes kernel selection and optimization.

Hawkes Processes

A point process with a sequence of n time-events is expressed by a vector of the form ($t_{1}$,$t_{2}$, … , $t_{n}$). Treating the real line as a time axis, the vector can be intuitively associated with a counting process $N(t)$, such that $dN(t)=1$, if there is an event at time t; and $dN(t)=0$, otherwise.

A point process can be described through its intensity function ($\lambda(t)$), which can be understood as the instantaneous expected rate of arrival of events, or the expectation of derivative of the counting process $N(t)$:

[TABLE]

This intensity function uniquely characterizes the finite-dimensional distributions of the point process (?). A simple example of this function would be the constant mean rate of arrival, $\mu$, in the case of a homogeneous Poisson process.

HPs model the intensity function in terms of self-excitation: the arrival of an event makes subsequent arrivals more likely to happen (?). HPs can be described through the following conditional intensity function $\lambda(t)$:

[TABLE]

where

• •

$\mathcal{H}(t)$ is the history of the process, the set containing all the events up to time t;

• •

$\mu$ is called background rate, or exogenous intensity, which is fixed as the mean rate of a homogeneous Poisson process;

• •

$\phi(t)$ is denominated self-triggering kernel, or excitation function.

From this function, one may notice that the intensity at time t will likely be affected by events which happened before the time t, described by the history of the process. From (?), we have that, if:

[TABLE]

then the corresponding process will show wide-sense stationary behavior, from which the asymptotic steady arrival rate, or first-order statistics, $\Lambda=\tfrac{\mu}{(1-||\phi||)}$, can be obtained, along with its covariance function, or second-order statistics, which is independent of t, $\nu(\tau)=\mathbb{E}[dN(t)dN(t+\tau)]$.

Estimating $\Lambda$ and $\nu(\tau)$ requires wide-sense stationarity assumptions which, besides being analytically convenient, are also connected to the fact that, in real data, the chain of self-excitedly induced further events will always be of finite type, or without ‘blowing up.’ This corroborates the practicality of the estimated model.

Discretized Kernel Estimation

Being one possible way of recovering the triggering kernel of a HP, this step is fully described in (?), and basically consists of building an estimator of $\phi(t)$ from empirical measurements of $\nu(\tau)$, the stationary covariance.

Given a finite sequence of ordered time-events in $[0,T]$, we fix a window size of h, and estimate $\nu(\tau)$ as:

[TABLE]

In practice, this estimation is done in discrete time steps $\delta$, up to a maximum value of $\tau$. 111In our case, we used a carefully designed heuristics explained in the section Experimental Results.:

[TABLE]

where $dN_{i\delta}^{(h)}$ is the total number of events happening between $t=i\delta$ and $t=i\delta+h$.

From (?), we have that, given $g_{t}^{(h)}=(1-\frac{|h|}{t})^{+}$, i.e., a triangular kernel density estimator with bandwidth h, we have the following relation in Laplace domain: $\hat{\nu_{z}^{(h)}}=\hat{g_{z}}^{(h)}(1+\hat{\psi}_{z}^{\star})\Lambda(1+\hat{\psi}_{z}^{\star})^{{\dagger}},$ where222Given a function $f_{t}$, $\hat{f}_{z}$ is its Laplace Transform, and the “$\star$” symbol corresponds to its conjugate.:

[TABLE]

Working with the Fourier transform restriction, i.e., ($z=i\omega$, with $\omega\in\mathbb{R}$) and given that $\hat{g}_{i\omega}^{(h)}=\tfrac{4}{\omega^{2}h}\sin^{2}(\tfrac{\omega h}{2}),$ we get to

[TABLE]

where we fix $h=\delta$, so we do not bother with cancellations of $\hat{g}_{z}^{(h)}$. Then, from: $|1+\hat{\psi}_{i\omega}|^{2}{=}\tfrac{\hat{\nu}_{z}^{(h)}}{\Lambda\hat{g_{z}}^{(h)}},$ we get to the discretized estimation of $\phi_{t}$ by taking the inverse Fourier transform of:

[TABLE]

in which the operator $H(\cdot)$ refers to the Hilbert transform.

Automatic Kernel Decomposition for HPs

This section presents the second step of the automatic kernel identification: a parametric kernel search through our new kernel decomposition scheme.

Self-Exciting Kernels

From the definition of the conditional intensity function, the self-excitation of the process is expressed through the kernel function $\phi(t)$. For the kernel decomposition, four base kernels will be used for identifying and estimating typical triggering behaviors as shown in Table 1:

• •

EXP($\alpha$,$\beta$): The decay exponential kernel is parameterized by the amplitude $\alpha$ and decay rate $\beta$, and is useful for modeling quick influence decay, in which initial transactions/hyperlinks have a lot of impact initially but rapidly reduce their influence over time;

• •

PWL(K,c,p): The power law kernel is parameterized by the amplitude K, the exponent p, and the constant c, modeling a slower decaying trend than the exponential;

• •

SQR(B,L): The pulse kernel is described by the amplitude B and the length L. Being a trivial, steady, and self-exciting dynamics on its own, it may also work as an offset level for the combined triggering with other kernel types, in the case of addition, and as a horizon truncation, in the case of multiplication333$u(t)$ is the step function.;

• •

SNS($A$,$\omega$): A truncated sinusoidal kernel, parameterized by the amplitude A and the angular velocity $\omega$. This type of kernel base function captures well the self-excitement of periodic events.

Here, the discretized kernel estimation is optional when a direct optimization of kernel structure is possible. Unfortunately, discontinuous functions (SQR, SNS) do not allow such optimization (e.g., Gradient Descent, Nelder-Mead). In this paper, we use the discretized kernel estimation as a unified method for both continuous (EXP, PWL) and discontinuous (SQR, SNS) kernels; and, most importantly, their combinations.

Furthermore, another great advantage of this step, compared with traditional sequential methods, is the fact that the value of $\nu$ for each value of $\tau$ can be calculated independently, while, in Gradient Descent, the value of the parameters at step t must be obtained before the values for step $t+1$. When combined with the parallelization of loops in our algorithm, this step significantly improves the speed of obtaining the most likely parametric representations of the sample processes.

Kernel Decomposition

For expressing the discretized estimation in terms of the four base kernels, the following steps are executed:

Calculate residues ($L^{1}$-error) w.r.t. the four basic kernels $\{$EXP, PWL, SQR, SNS$\}$; 2. 2.

Select the kernel with the minimum residue $MR_{1}$, denominated $K_{1}$; 3. 3.

Check whether the estimated parameters of the kernel satisfy the stationarity condition, by using the closed-form expressions from Table 2; 4. 4.

Calculate residues w.r.t. a total of 8 kernel expansions, resulting from 2 operations (addition and multiplication) per base kernel $\{$+EXP, $\times$EXP, +PWL, $\times$PWL, +SQR, $\times$SQR, +SNS, $\times$SNS$\}$, while fixing the optimized parameters for $K_{1}$, in the case of Additive Combination, and recalculating all the parameters, in the case of Multiplicative Combination; 5. 5.

Select the kernel with minimum residue $MR_{2}$, denominated $K_{2}$, and check the spectral radius condition (calculated in closed-form from Table 3); 6. 6.

If both K1 and K2 are stable, and $MR_{1}<MR_{2}/\eta$ ($\eta$ would act as a regularization parameter), pick $K_{1}$ . Else, pick $K_{2}$. 7. 7.

If likelihood (llh) of direct optimization (GD, Nelder-Mead) is greater than likelihood of kernel decomposition, output GD model. Else, output the decomposition model.

Regarding the computational efficiency of the decomposition algorithm, two strategies yielded results at a much lower computational cost, without altering the results of the decomposition:

• •

Selecting the best kernel through the error, instead of likelihood;

• •

Greedy search of K2 based on the selected K1, instead of doing a brute-force search over all the $4\times 8=32$ possible combinations for K2.

Figure 5 explains the algorithm up to the depth two, for illustration purposes. Our kernel decomposition scheme can be expanded into multiple depths, as explained in Section Higher-order Kernel Decomposition. Our algorithm is presented in Algorithm 1.

Stationarity Conditions

Verifying the stationarity condition is one of the most important steps in the kernel search. When we end up with a non-stationary kernel, estimating future events can not be accurate.

To solve this issue, we developed closed-form expressions, either in the form of equality or as an upper bound, which are shown in Table 2, for the case of a single kernel, and Table 3, for multiplicative combinations of two kernels 444 $\Gamma(\cdot,\cdot)$ is the well-known incomplete Gamma function: $\Gamma(a,y)=\int_{y}^{\infty}t^{a-1}e^{-t}dt$. The conditions for additive combinations can be derived from the conditions for single kernels in a straightforward manner.

The kernel is said to induce stationarity if the result of the expression calculated using the estimated parameters belongs to the interval [0,1). This can be justified both from the point-of-view of HP as a branching process, also called immigrant-birth representation (?), and of the boundedness of the spectral radius (largest absolute value among the eigenvalues) of the excitation matrix. 555For the univariate HP case, the excitation matrix has dimension one, being only the excitation function, $\phi(t)$.

Scale-Independence Criterion

For an automatic time series analysis, scale-independence is indispensable, as time sequences of disjoint datasets may occur in time scales differing by several orders of magnitude. As an example, earthquake events’ occurrences in a sequence are spaced by intervals of monthly and yearly scales. Thus, setting a horizon of a few months as the maximum value of $\tau$ in Equation (4) might result in a satisfactory discrete estimation grid. However, using the same time length for estimating the triggering behavior of a finance-related sequence would require an impractically large grid resolution.

A histogram of all the time intervals between events in a sequence may be readily generated, and is an indicator of the overall magnitude of the spacing among the events. Thus, as a rule of thumb, the horizon length for $\tau$ may be set as the smaller time interval strictly larger than a percentage of the sequence’s intervals. The values of 50 % and 95 % were used. In practice, this value of horizon length is obtained with the help of a histogram composed by 100 bins.

Higher-order Kernel Decomposition

A sequential additive decomposition of the discretized estimation vector is rather straightforward, since one may just set the residual vector from the previous stages as the input of the next ones.

In the case of multiplicative decomposition, it is nontrivial to find the result of intraclass decomposition. To the best of our knowledge, no analysis on multiplicative HP kernel decomposition is reported yet.

In this paper, we provide a new upper bound over an interclass kernel product of unknown degree, as in:

[TABLE]

for $k_{i}\in\mathbb{Z}^{*}$, where the operator “$\left[\cdot\right]^{k}$” corresponds to the set of functions which can be decomposed into a k-th order product of kernels, e.g:

[TABLE]

By deriving the four possible intraclass kernel products, one may observe that the typical self-exciting behavior features of each kernel type are preserved, as in the following:

• •

$\left[\mbox{EXP}\right]^{k_{1}}$ reduces to the case of a single exponential with $\alpha=\prod_{i=1}^{k_{1}}\alpha_{i}$ and $\beta=\sum_{i=1}^{k_{1}}\beta_{i}$, thus still accounting for its ‘quick-decay’ behavior: $\left[\mbox{EXP}\right]^{k_{1}}\subset\left[\mbox{EXP}\right]$

• •

$\left[\mbox{PWL}\right]^{k_{2}}$ is lower bounded by a single PWL kernel with $K=\prod_{i=1}^{k_{2}}K_{i}$, $c=max(c_{1},...,c_{k_{2})}$ and $p=\sum_{i=1}^{k_{2}}p_{i}$, thus still accounting for its ‘slow-decay’ behavior

• •

$\left[\mbox{SQR}\right]^{k_{3}}$ reduces to a single SQR kernel with $B=\prod_{i=1}^{k_{4}}B_{i}$ and $L=min(L_{1},...,L_{k_{4}})$, thus still accounting for its ‘steady-triggering’ behavior: $\left[\mbox{SQR}\right]^{k_{3}}\subset\left[\mbox{SQR}\right]$

• •

$\left[\mbox{SNS}\right]^{k_{4}}$ has $A=\prod_{i=1}^{k_{4}}A_{i}$ and a ‘spikier’ aspect (higher bandwidth), thus still accounting for its ‘periodicity-inducing’ behavior

Thus, on deepening the decomposition algorithm by overly increasing the number of levels above 2, we may be, in fact, adding little information on the qualitative aspect of the self-exciting behavior analysis of the data while making it more prone to overfitting to the noisiness of the discretized estimation vectors.

Upper Bound

Furthermore, regarding the boundedness of the higher-order decompositions, from the exact results for EXP and SQR intraclass decompositions and the upper bounds for the PWL and SNS ones, we have that:

[TABLE]

for $0\leq x\leq min(L,\dfrac{\pi}{\omega})$, and 0 otherwise.

Experimental Results

To demonstrate the benefits of the kernel decomposition framework, we conducted experiments with synthetic, financial and earthquake data.

For real-world data sets, no prior information about the kernel (type and parameters) is available. Thus, we use the log-likelihood of the kernel function over the time sequence as a quality criterion.

Given a realization $(t_{1},t_{2},...,t_{k})$ of some regular point process on [0,T], its log-likelihood (l) is expressed as: $l(t_{1},...,t_{k})=\sum_{i=1}^{k}\log(\lambda(t_{i}))-\int_{0}^{T}\lambda(u)du$.

For an automatic time series analysis, scale-independence is indispensable, as time sequences of disjoint datasets may occur in time scales differing by several orders of magnitude. As an example, earthquake events’ occurrences in a sequence are spaced by intervals of monthly and yearly scales. Thus, setting a horizon of a few months as the maximum value of $\tau$ in Equation (4) might result in a satisfactory discrete estimation grid. However, using the same time length for estimating the triggering behavior of a finance-related sequence would require an impractically large grid resolution.

A histogram of all the time intervals between events in a sequence may be readily generated, and is an indicator of the overall magnitude of the spacing among the events. Thus, as a rule of thumb, the horizon length for $\tau$ may be set as the smaller time interval strictly larger than a percentage of the sequence’s intervals. The values of 50 % and 95 % were used. In practice, this value of horizon length is obtained with the help of a histogram composed by 100 bins.

Financial Data

In the finance domain, HPs have become more prevalent, due to its structure being naturally adapted to model systems in which the discrete nature of the jumps in $N_{t}$ is relevant, making the model remarkably suited to modeling high-frequency data (?).

Here, we picked the 19 top-varying companies of the Technology, Healthcare, Industrial, Services and Utilities categories of Yahoo Finance. We extracted tick data from every two minutes of 30 business days (02/02/2017 to 02/23/2017 for Technology and 04/07/2017 to 05/18/2017 for the other ones). Whenever a stock price changed by some magnitude higher than some threshold, an event was logged in the corresponding time sequence. Ten different percentual thresholds, increasing at equally spaced intervals from 0.03% to 0.3%, were applied. This procedure resulted in a number of valid sequences, for each category, indicated in Table 5, since the remaining ones did not contain enough points for the splitting between training and validation subsequences.

As in extrapolation tasks, the 80% of the first elements for each sequence were then used as training data, and the remaining 20% were used for validation, i.e., we estimated the parameters of the kernel using the first 24 days and then calculated the log-likelihood on the last 6 days of each sequence. The kernel was then normalized to 2 min = 120 sec. When comparing the log-likelihoods of first and second level decompositions, we observed that the second level, with composite kernels, resulted in a higher log-likelihood in the majority of sequences from each category, as indicated in Table 5, what corroborates that a more flexible model of the kernel provides a more accurate description of the underlying dynamics of the process. The average log-likelihood for each level is shown in Table 4. The comparison of each sequence is shown in Figure 7.

When comparing the performance of the best estimation among the two levels and the usual exponential HP model used in financial analysis, fitted through the gradient-based method from (?), it is possible to see that the kernel composition exhibited a much more robust performance. Although the exponential HP performed well in some sequences, it tended to get stuck in local maxima with very poor performance, usually leading to unstable or negative combinations of parameters, for which the likelihood is null. The kernel composition performs better in the majority of sequences, as shown in Table 5. We provide the comparison of each individual sequence in the supplementary material.

Earthquake Data

The data considered for the earthquake experiment was a set of 100 time sequences extracted from the USGS NCSN Catalog (NCEDC database), from the day of 01/01/1966 to 01/01/2015. The latitude range was [30,55], and the longitude range was [-140,-110]. Different length intervals and resulting areas were considered. Whenever the magnitude of an event exceeded some threshold, its time coordinate was added to the corresponding input time sequence. The magnitude thresholds were varied among 2.5, 3.0, 3.5 and 4.0; and the grid resolution was set to 20 and 100 points.

Seeking scale-independent search, we use the aforementioned histogram heuristics: earthquakes events are separated by time intervals of monthly or yearly scales. Thus, the estimation horizon for financial data, for example, lasting usually only a few seconds, would hardly capture the overall aspect of the triggering behavior in this case.

The results indicate a strong agreement with the long standing assumption of a power-law shaped kernel for the intensity of aftershocks’ occurrences (‘Omori’s Law’ (1894)). For 20-point grid resolution, the relative frequency of each kernel was $\mbox{(EXP, PWL, SQR, SNS)}=\mbox{(0, 97, 2, 1)}$. For the 100-point grid resolution, the relative frequency was (0, 99, 1, 0). Q-Q plots from the estimated models are shown in Figure 8, in which comparisons to the original sequence are made among sequences generated by our kernel composition, disconsidering the stabiliy check, the discretized estimate and the usual Power-Law kernel model, fitted through the gradient descent based method (GD) (?). Both our method and the discretized estimate perform very close to the original sequence, while the GD method tended to get stuck in local optima with poor performance. A timescale-based initialization of $\mu$ was used.

Verifying Scale-Independence Criterion

To verify the histogram criteria (explained in Section Scale-Independence Criterion), used the data sets ( Stocks and Earthquake). As shown in Figure 9, The histogram criteria allows us to find a good resolution of kernels in highly different scales.

Conclusion

Hawkes processes are point processes which capture self-exciting discrete events in time series data. To predict future events with HPs, an appropriate kernel is selected by hands, previously. In this paper, we proposed a new temporal covariance-based kernel decomposition method to represent various self-exciting behaviors. We also presented a model (structure/parameter) learning algorithm to select the best HP kernel given the temporal discrete events. The stationarity conditions are derived to guarantee the validity of the kernel learning algorithm. In experiments, we demonstrated that the proposed algorithm performs better than existing methods to predict future events by automatically selecting kernels.

Appendix A Derivations of Stationarity Criteria for Multiplicative Combinations of Kernels

This appendix introduces the full derivations of stationarity criteria for the second order multiplicative compositions of the four base kernels.

EXP x EXP

For the combination “EXPxEXP”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Thus:

[TABLE]

So, this case reduces to the case of a single exponential.

EXP x PWL

For the combination “EXPxPWL”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Thus:

[TABLE]

where $\Gamma(\cdot,\cdot)$ is the well-known Incomplete Gamma Function: $\Gamma(a,y)=\int_{y}^{\infty}t^{a-1}e^{-t}dt$.

EXP x SQR

For the combination “EXPxSQR”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Thus:

[TABLE]

So, in the case of a multiplicative combination, the SQR kernel acts as a truncation horizon.

EXP x SNS

For the combination “EXPxSNS”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Where:

[TABLE]

PWL x PWL

In the case of the combination “PWLxPWL”, an upper bound is derived as follows:

[TABLE]

[TABLE]

Then:

[TABLE]

PWL x SQR

For the combination “PWLxSQR”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Where:

[TABLE]

So, once again, the SQR kernel acts as a truncation horizon.

PWL x SNS

In the case of the combination “PWLxSNS”, an upper bound is derived as follows:

[TABLE]

[TABLE]

Where:

[TABLE]

SQR x SQR

For the combination “SQRxSQR”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Where:

[TABLE]

So, the multiplicative combination of two SQR kernels may be reduced to the case of a single SQR kernel.

SQR x SNS

In the case of combinations of discontinuous kernels (SQR and SNS), we assume they have the same starting and ending points, i.e., $L=\dfrac{\pi}{\omega}$. So, for the combination “SQRxSNS”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Where:

[TABLE]

SNS x SNS

In the case of combinations of discontinuous kernels (SQR and SNS), we assume they have the same starting and ending points. So, for the combination “SNSxSNS”, we have that, for stationarity to be achieved:

[TABLE]

[TABLE]

Where:

[TABLE]

Appendix B Derivation of the Log-likelihood formula for HPs

This derivation follows the steps on (?). Given a realization $(t_{1},t_{2},...,t_{k})$ of some regular point process observed over the interval [0,T], the log-likelihood is expressed as:

[TABLE]

Proof.

Let be the joint probability density of the realization:

[TABLE]

It can be written in terms of the Conditional Intensity Function. We can then find f in terms of $\lambda$:

[TABLE]

where, given the history up to last arrival u, $\mathcal{H}(u)$, F(t) is then defined as the conditional cumulative probability distribution of the next arrival time $T_{k+1}$:

[TABLE]

Integrating both sides of Equation (B) over $(t_{k},t)$:

[TABLE]

Given that the realization is assumed to have come from a so-called simple process, i.e., a process in which multiple arrivals cannot occur at the same time, we have that $F(t_{k})$ = 0 as $T_{k+1}>t_{k}$, which simplifies equation (B) to:

[TABLE]

Further rearranging the expression:

[TABLE]

and

[TABLE]

Thus, the likelihood becomes:

[TABLE]

Given that the process is observed on [0,T], the likelihood must include the probability of seeing no arrivals in $(t_{k},T]$:

[TABLE]

Through using the formulation of F(t), we have that:

[TABLE]

Finally, getting the logarithm of the expression, we have the formula for l:

[TABLE]

∎

Appendix C Automatic Report

Appendix D Comparison between Gradient-based and Discretized Estimation steps for the financial datasets