Hawkes processes for credit indices time series analysis: How random are trades arrival times?

TL;DR

This paper applies Hawkes processes to model and analyze the timing of trades in credit indices, capturing self-excitement and interactions, and introduces a new fitting method validated on simulated and real data.

Contribution

It presents a simple, effective maximum likelihood method for fitting multidimensional Hawkes processes with exponential kernels to credit trading data.

Findings

Successfully fitted Hawkes models to real credit index data

Quantified self-excitement and cross-influence among indices

Validated the method on simulated data

Abstract

Targeting a better understanding of credit market dynamics, the authors have studied a stochastic model named the Hawkes process. Describing trades arrival times, this kind of model allows for the capture of self-excitement and mutual interactions phenomena. The authors propose here a simple yet conclusive method for fitting multidimensional Hawkes processes with exponential kernels, based on a maximum likelihood non-convex optimization. The method was successfully tested on simulated data, then used on new publicly available real trading data for three European credit indices, thus enabling quantification of self-excitement as well as volume impacts or cross indices influences.

Figures

Tables

Table 1. Table 1: Statistics on inter-arrival times Δ ​ t Δ 𝑡 \Delta t (period: 03/01/2017 to 14/12/2017)

	$𝔼\left(\Delta t\right)$	$\sigma \left(\Delta t\right)$	$Q_1\left(\Delta t\right)$	$Q_2\left(\Delta t\right)$	$Q_3\left(\Delta t\right)$	Number of trades
itxeb	7 min	10 min	1 min	3 min	8 min	22 227
itxex	8 min	12 min	2 min	4 min	10 min	18 152
itxes	17 min	24 min	3 min	9 min	21 min	7 194

Table 2. Table 2: Estimated parameters for ITXEB , ITXES and ITXEX over the period (2017-01-15 to 2017-12-15)

	Estimated ${\mu }^{-1}$	Estimated $\alpha$	Estimated ${\beta }^{-1}$
ITXEB	22 min	0.62	20 min
ITXEX	24 min	0.60	23 min
ITXES	58 min	0.65	36 min

Equations

$$\lambda(t|\mathscr{F}_{t})=\lim_{h\rightarrow 0}\mathbb{E}\bigg{[}\frac{N(t+h)-N(t)}{h}\bigg{|}\mathscr{F}_{t}\bigg{]}$$

$$\lambda (t)=\mu (t)+{\int }_{-\infty }^t\varphi (t-\tau )dN(\tau )=\mu (t)+\sum _{t_i<t}\varphi (t-t_i)$$

$${\lambda }_i(t)={\mu }_i+\sum _{j=1}^M{\int }_0^t{\varphi }_{i,j}(t-\tau )d{N}_j(\tau )={\mu }_i+\sum _{j=1}^M\sum _{n=1}^{N_j(t)}{\varphi }_{i,j}(t-t_{j,n}),$$

$${\mathcal{L}}_t(\mathbf{\lambda }):=-\sum _{i=1}^M\left({\int }_0^t\log {\lambda }_i(\tau )d{N}_i(\tau )-{\int }_0^t{\lambda }_i(\tau )d\tau \right).$$

$$\mathcal{L}_{t}(\bm{\mu},\bm{\alpha},\bm{\beta})=\sum_{m=1}^{M}\bigg{(}\mu_{m}T\ +\ \sum_{n=1}^{M}\alpha_{mn}\sum_{\{k:t_{k}^{n}<T\}}[1-e^{-\beta_{mn}(T-t_{k}^{n})}]\\ -\sum_{\{k:t_{k}^{n}<T\}}log[\mu_{m}+\sum_{n=1}^{M}\alpha_{mn}\beta_{mn}\sum_{\{k:t_{k}^{n}<t_{i}^{m}\}}e^{-\beta_{mn}(t_{i}^{m}-t_{k}^{n})}]\bigg{)}$$

$${\mathcal{L}}_t^{\ast }(\mathbf{\beta })={\mathcal{L}}_t({\mathbf{\mu }}^{\ast }(\mathbf{\beta }),{\mathbf{\alpha }}^{\ast }(\mathbf{\beta }),\mathbf{\beta })=mi{n}_{(\mathbf{\mu },\mathbf{\alpha })}{\mathcal{L}}_t(\mathbf{\mu },\mathbf{\alpha },\mathbf{\beta })$$

$$mi{n}_{(\mathbf{\mu },\mathbf{\alpha },\mathbf{\beta })}{\mathcal{L}}_t(\mathbf{\mu },\mathbf{\alpha },\mathbf{\beta })=mi{n}_{(\mathbf{\beta })}mi{n}_{(\mathbf{\mu },\mathbf{\alpha })}{\mathcal{L}}_t(\mathbf{\mu },\mathbf{\alpha },\mathbf{\beta })$$

$$\mu=\bigl{(}\begin{smallmatrix}0.1\\ 0.2\end{smallmatrix}\bigr{)}\quad\alpha=\bigl{(}\begin{smallmatrix}0.5&0.00\\ 0.4&0.3\end{smallmatrix}\bigr{)}\quad\beta=\bigl{(}\begin{smallmatrix}0.3&0.00\\ 0.2&0.2\end{smallmatrix}\bigr{)}$$

$$\frac{1}{\mu}=\bigl{(}\begin{smallmatrix}40\\ 42\\ 38\end{smallmatrix}\bigr{)}\quad\alpha=\bigl{(}\begin{smallmatrix}0.40&0.07&0.57\\ 0.00&0.53&0.21\\ 0.00&0.16&0.59\end{smallmatrix}\bigr{)}\quad\frac{1}{\beta}=\bigl{(}\begin{smallmatrix}12&21&53\\ 71&18&12\\ 53&14&34\end{smallmatrix}\bigr{)}$$

$$\frac{1}{\mu}=\bigl{(}\begin{smallmatrix}29\\ 30\\ 125\end{smallmatrix}\bigr{)}\quad\alpha=\bigl{(}\begin{smallmatrix}0.44&0.25&0.20\\ 0.23&0.37&0.18\\ 0.07&0.08&0.45\end{smallmatrix}\bigr{)}\quad\frac{1}{\beta}=\bigl{(}\begin{smallmatrix}23&9&21\\ 14&13&31\\ 39&17&25\end{smallmatrix}\bigr{)}$$

$$(0,\infty )=(0,{\tau }_1]\cup ({\tau }_1,{\tau }_2]\cup \cdots \cup ({\tau }_{d-1},{\tau }_d]\cup \dots$$

$$\forall t\ge 0,\lambda (t)=\mu +\tilde{\lambda }(t)$$

$$\forall t\in ]{\tau }_{d-1},{\tau }_d]\tilde{\lambda }(t)=\pi \tilde{\lambda }({\tau }_{d-1})e^{-\rho (t-{\tau }_{d-1})}+{\int }_{[{\tau }_{d-1},t)}\alpha e^{-\beta (t-u)}dN(u)$$

$${\mathcal{L}}_T(\lambda )=-\left({\int }_0^T\log \lambda (t)d{N}_i(t)-{\int }_0^T\lambda (t)dt\right)$$

$${\mathcal{L}}_T(\lambda )=\mu T+\sum _{d\in days}\pi \tilde{\lambda }({\tau }_{d-1})(1-e^{-\rho ({\tau }_d-{\tau }_{d-1})})$$

$$+\alpha \sum _{t\in ]{\tau }_{d-1},{\tau }_d]}(1-e^{-\rho (t-{\tau }_{d-1})})-\sum _{t_i}\log \lambda (t_i)$$

$$\forall i\in [1\dots D],{\lambda }_i(t)={\mu }_i+\sum _{j=1}^D\int {\varphi }_{ij}(t-s)d{N}_j(s)$$

$${\varphi }_{ij}(t)={\alpha }_{ij}{\beta }_{ij}\exp (-{\beta }_{ij}t)1_{t>0}$$

$$\forall t\text{ between two consecutive days,}{\lambda }_i(t)=0$$

$${\mathcal{L}}_t(\mathbf{\lambda }):=-\sum _{i=1}^M\left({\int }_0^t\log {\lambda }_i(\tau )d{N}_i(\tau )-{\int }_0^t{\lambda }_i(\tau )d\tau \right).$$

$${\mathcal{L}}_t(\mathbf{\lambda }):=D\ast \delta \sum _{m=1}^M{\mu }_m+\sum _{m=1}^M\sum _{n=1}^M{\alpha }_{mn}\sum _{d\in days}\sum _{t_i^n\in \text{day d}}(1-e^{-{\beta }_{mn}({\tau }_d-t_i^n)})$$

$$-\sum _{m=1}^M\sum _{t_i^m}log({\mu }_m+\sum _{n=1}^M{\alpha }_{mn}{\beta }_{mn}{R}_{mn}(i))$$

$$R_{mn}(i)$$

$$=e^{-{\beta }_{mn}(t_i^m-t_{i-1}^m)}{R}_{mn}(i-1)+\sum _{\{k:t_{i-1}^m\le t_k^n<t_i^m\}}{e}^{-{\beta }_{mn}(t_i^m-t_k^n)}$$

$$\frac{\partial {\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\mu }_m}$$

$$\frac{\partial {\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\alpha }_{mn}}$$

$$\frac{{\partial }^2{\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\mu }_m\partial {\mu }_m}$$

$$\frac{{\partial }^2{\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\mu }_m\partial {\alpha }_{ml}}$$

$$\frac{{\partial }^2{\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\alpha }_{ml}\partial {\mu }_m}$$

$$\frac{{\partial }^2{\mathcal{L}}_t(\mathbf{\lambda })}{\partial {\alpha }_{mk}\partial {\alpha }_{ml}}$$

Full text

1 Introduction

1.1 From credit derivative indices to Hawkes processes

Credit indices are financial instruments comprised of a set of credit securities, mainly used to hedge credit default risk. Each new index series is issued every 6 months and expires after a defined maturity. Though liquid, those indices are traded at a rather ”mid-frequency” rate, with trades occurring at a minute scale. For European indices, the market is continuously open from 7:00 to 17:00, with a total of five to ten billion euros reported each day on average. Fig 1 gives an idea of the activity on a sampled day of trading. If recent regulations have led to a greater public reporting of trading activities, thus releasing more amount of traded data, this data is by essence quite sparse. Picturing such market behaviour over time represents a challenging opportunity: it originally motivated the work presented in this article.

This study focuses on three principal European credit indices, of the most traded 5-year maturity111Data provided by otcstreaming.com, and available on demand. : the Main Index, noted here itxeb, a combination of 125 equally weighted investment grade entities; the Crossover Index, noted itxex, with 75 sub-investment grade names; and the Senior Financial, noted itxes, a subset of 30 financial entities from Main Index, referencing senior debt.

When trying to capture the underlying dynamics of this market, one of the first questions that comes to mind is: how random is the timing of trades? A very naive approach is to model trades arrival times as purely random and uncorrelated processes, for example a Poisson process by trying to fit an exponential law over the distribution of inter-arrival times. Yet, such a model fails to fit the data: from Fig 2, we can notice a high density of very short inter-arrival times, suggesting that some self-excitation phenomenon occurs, with trades triggering more trades; an intuition that every practitioner of the field would confirm.

Literature review reveals that one model has recently driven considerable attention and a growing interest from both the scientific and quantitative communities: Hawkes processes. This model, introduced in 1971 by Hawkes ([6], [7]) to characterize earthquake tremors, can be used to describe timing of trades and their cross influences from a point process perspective.

After summarizing related work for Hawkes processes, we will recall the theoretical framework around this stochastic model. We will then focus on optimization methods for fitting such a process and describe Two Stage Hawkes Likelihood Optimization (2SHLO), a maximum likelihood based algorithm to fit multidimensional Hawkes processes with a parametric exponential kernel. Lastly, we will present our results, firstly on simulated data, then on credit derivative indices mentioned above. Hawkes processes in finance have mainly been applied to high frequency data. Is such a model appropriate for ”mid-frequency” credit data? Can we specify the impact of volumes in trades or describe mutual influences between different indices?

1.2 Related work

Hawkes processes were originally designed for seismology analysis and deeply studied for that purpose by Ogota [12], [10], [11]. The concept has been utilized to model other effects where self and cross ignitions happen, such as social media tweets cascading [14], crime occurrences due to gang retaliations [9], or financial market events. As regards to financial applications, Bowsher [4] proposed in 2002 a generalized Hawkes process model taking into account night gaps, inter-day dependencies as well as intraday seasonality (with a piece-wise baseline intensity), which was used to model interactions between trades and price changes and also estimate of the price volatility based on mid-quote intensity for NYSE stock. Market price microstructure & market impact (influence of market orders on forthcoming prices) was modelled in [3], price impact was also studied in [1]. [13] focused on the impact of volume and order types on the limit order book. [2] (2015) exposes a full review of Hawkes process applications to finance, such as price & volatility modeling, market reflexivity measurement, order book modeling or risk contagion modeling.

2 Modeling Self and Cross Excitement with Hawkes Processes

This section re-frames the required formalism, as very clearly stated in [2], [8] or [15], starting from point and counting process definitions, through the key concept of intensity function. Hawkes model formulation is detailed for the reader, with a peculiar focus on exponential kernel structure.

2.1 Core concepts: from counting & point processes to intensity functions

Definition 1 (Point Process)

Let $(\Omega,\mathscr{F},P)$ be a probability space.

Let $(t_{k})_{k\in\mathbb{N}^{*}}$ be a sequence of non-negative random variables such that $\forall\ k\in\mathbb{N}^{*}$, $t_{k}<t_{k+1}$. $(t_{k})_{k\in\mathbb{N}^{*}}$ is called a (simple) point process on $\mathbb{R}_{+}$.

Definition 2 (Counting Process)

Let $(t_{k})_{k\in\mathbb{N}}$ be a point process. The right-continuous stochastic process defined for all $t\in\mathbb{R}_{+}$ as $N(t)=\sum_{k\in\mathbb{N}^{*}}1_{t_{k}\leq t}$ is called the counting process associated with $(t_{k})_{k\in\mathbb{N}^{*}}$.

The study of point processes goes through a single mathematical object, the intensity function, which is defined as the conditional probability density of occurrence of an event in the immediate future.

Definition 3 (Intensity Function)

Let $N$ be a counting process adapted to a filtration $\mathscr{F}_{t}$. The left-continuous intensity is heuristically defined as

[TABLE]

The intensity function depends on the choice of filtration $\mathscr{F}_{t}$, which represents the amount of information available until time $t$ (See [5] for a rigorous definition of the intensity function). In the context of Hawkes processes, we will simply use the natural filtration, with all previous information being available, and consider $\lambda(t)$.

Homogeneous Poisson process

One of the simplest point processes is the homogeneous Poisson process, for which the intensity function is constant over time: $\forall\ t\geq 0,\ \lambda(t)=\lambda$. In that case, durations (or inter-event waiting times) are independent and identically distributed (following an exponential distribution of hazard rate $\lambda$). As presented in the introduction, credit trades cannot be modeled by this memory-less model. Let’s introduce Hawkes processes, a peculiar kind of non-homogeneous Poisson processes with linear dependencies over functions of past events.

2.2 Self-excitement: one-dimensional Hawkes processes

Definition 4 (One-dimensional Hawkes Processes)

Let $(t_{k})_{k\in\mathbb{N}}$ be a point process and $N$ the associated counting process, such that its intensity function $\lambda$ is defined for each time $t\geq 0$ as

[TABLE]

where $\mu:\mathbb{R}\mapsto\mathbb{R}_{+}$ is an exogenous base intensity and $\phi:\mathbb{R}_{+}\mapsto\mathbb{R}^{+}$ is a non-negative, measurable function such that $||\phi||_{1}=\int_{0}^{\infty}\phi(s)ds<1$.

$N$ is called a Hawkes process with baseline $\mu$ and kernel $\phi$.

The kernel $\phi$ expresses the positive influence of past events on the current value of the intensity. Each jump $dN(\tau)\neq 0$ increases the probability of future events through the kernel $\phi$. Clustering effects and branching structure are well depicted by Hawkes processes, with the baseline activity generating immigrant events and descendant events enhancing the intensity.

For this study, we focus on exponential kernels, defined as $\phi(t)=\alpha\beta\ e^{-\beta t}\ \mathds{1}_{t\geq 0}$, with parameters $\alpha,\beta\geq 0$, which allows easy interpretation: $\alpha$ or adjacency represents the weight of previous events while $\beta$ is the decay / typical duration of influence for a past event.

The branching ratio $n=||\phi||_{1}=\int_{0}^{\infty}\phi(s)ds$ represents the average number of descendants for any event. For exponential kernel, $n=\alpha$.

2.3 Including mutual-excitement with multidimensional Hawkes processes

Definition 5 (Multidimensional Hawkes Processes)

Let $M\in\mathbb{N}^{*}$ and $\{(t_{k}^{i})_{k}\}_{i=1,...,M}$ be a M-dimensional point process.

We denote by $N_{t}=(N_{t}^{1},...,N_{t}^{M})$ the associated counting process such that its vector intensity function $\lambda:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}^{M}$, is defined as, for all $t\geq 0$, $i\in{1,...,M}$:

[TABLE]

with $\mu=(\mu_{i})_{i=1,...,M}$ an exogenous base intensity vector

and $\phi(t)=(\phi_{i,j})_{i,j=1,...,M}$ a matrix-valued kernel that is component-wise positive and causal (null values when $t<0$) and with each component belonging to the space of $L^{1}$-integrable functions.

$N_{t}$ is called a multidimensional (or multivariate) Hawkes process with baseline $\mu$ and kernel $\phi$.

The choice of an exponential kernel can be generalized for multivariate Hawkes processes with kernel components being expressed as $\phi_{i,j}(t)=\alpha_{i,j}\beta_{i,j}e^{-\beta_{i,j}t}\mathds{1}_{t\geq 0}$.

This article focuses on Hawkes processes with exponential kernel and constant baselines, to which we propose a maximum likelihood based fitting method.

3 Hawkes Processes: Model Calibration

The most commonly used technique for parametric inference of Hawkes processes is a direct numerical optimization of the Maximum Likelihood, which was first introduced in the work of Ogata (1978) [10]. He proved the asymptotic consistency and efficiency of the Maximum Likelihood Estimator (MLE) under the assumption of stationary condition of the underlying point process. The negative log-likelihood function of a multidimensional Hawkes process over the time interval $[0,t]$ is given in Daley & Vere-Jones [5] Proposition 13.1.VI by :

[TABLE]

Depending on the shape of kernels, $\mathcal{L}_{t}$ is generally non-convex. Classic non-convex numerical optimization pitfalls are to be feared: a direct numerical optimization could converge to a merely local minimum. On the other hand, it may take too many iterations to converge as the negative log-likelihood function can be flat on some regions of the space of parameters. This problem is also faced when using an Expectation Maximization (EM) algorithm to find the MLE.

On the computational side, the repeated evaluation of the negative log-likelihood can be highly time consuming, essentially due to the nested sum in the first part of 4 where the conditional intensity is also expressed as a sum over the history.

3.1 Maximum Likelihood Estimation (MLE) for exponential multidimensional Hawkes processes

Let’s consider a M-dimensional multivariate Hawkes model with constant baselines $\mu_{i}$ and kernel functions $\phi_{i,j}$ as defined in 5. The choice of an exponential kernel has proved to be very interesting as it allows intuitive and meaningful interpretations of its parameters and also reduces the computational cost of the evaluation of the negative log-likelihood function as noted by Ogata (1981) [11] who exhibited a recursive formula that eliminated the nested sum evaluation problem.

The existing methods to fit exponential parameters through MLE directly use non-linear optimization algorithms such as Nelder-Mead (also called downhill simplex method) or BFGS, with performance decreasing as the number of parameters increases. This represents quite an issue as $M(2M+1)$ parameters describe a M-dimensional Hawkes process with exponential kernel.

In addition, other existing methods often make the strong assumption that the decays $\beta_{i,j}$ of the exponential kernel are given and fixed a priori. Consequently, the estimation of each kernel function is equivalent to the estimation of its adjacency coefficient $\alpha_{i,j}$ and since the conditional intensity $\lambda_{i}(t)$ is linear with respect to kernel functions $\phi_{i,1}(t),\ldots,\phi_{i,p}(t)$, the negative log-likelihood function $\mathcal{L}_{t}$ is convex with respect to all parameters. Therefore we can use the widely available convex optimization machinery to find the global minimum efficiently. The main drawback of this method is the possible inaccuracy of the decays parameters $\beta_{i,j}$ which can lead to model mismatch.

3.2 Introducing Two Stage Hawkes Likelihood Optimization (2SHLO):

In order to combine the benefits of both approaches described above, we propose the following estimation method to fit exponential Hawkes processes.

The negative log-likelihood function for exponential multivariate Hawkes processes is given by (see Ogata[11] for the recursive formulation):

[TABLE]

As we mentioned earlier, if $\beta$ is fixed, then $\mathcal{L}_{t}(\cdot,\cdot,\beta)$ is convex with respect to parameters $(\bm{\mu},\bm{\alpha})$. As a result it can be minimized using any convex optimization algorithm, for example an Accelerated Gradient Descent (AGD) which converges to a global minimum $(\bm{\mu}^{*}(\bm{\beta}),\bm{\alpha}^{*}(\bm{\beta}))$ such that:

[TABLE]

$\mathcal{L}_{t}^{*}(\bm{\beta})$ is still a non-convex function but defined on a space with a lower dimension than the initial parameter space. Any non-linear heuristic algorithms can now be used to minimize $\mathcal{L}_{t}^{*}(\bm{\beta})$ leading to the MLE as222if we call $m=min_{(\bm{\mu},\bm{\alpha},\bm{\beta})}\mathcal{L}_{t}(\bm{\mu},\bm{\alpha},\bm{\beta})$, the following inequality holds for all $(\bm{\mu},\bm{\alpha},\bm{\beta})$: $m\leq\mathcal{L}_{t}(\bm{\mu}^{*}(\bm{\beta}),\bm{\alpha}^{*}(\bm{\beta}),\bm{\beta})\leq\mathcal{L}_{t}(\bm{\mu},\bm{\alpha},\bm{\beta})$ 3.2 follows immediately by taking the minimum over $(\bm{\mu},\bm{\alpha},\bm{\beta})$ in each side of the inequality. :

[TABLE]

We summarize the fitting method as follow :

3.3 Goodness of fit

The compensator function $\Lambda$ of a point process with intensity $\lambda$ is defined as $\forall\ t\geq 0,\ \Lambda(t)=\int_{0}^{t}\lambda(s|\mathscr{F}_{t})ds$. To assess the goodness of fit of our estimation, we use the residual point process analysis theorem, as stated in [5]:

The transformed sequence $\{t^{*}_{k}\}=\{\Lambda(t_{k})\}$ is a realization of a unit rate Poisson process if and only if the original sequence $\{t_{k}\}$ is a realization from the point process defined by $\lambda$.

The goodness of fit can be tested by comparing the transformed estimated times to a standard Poisson process using Q–Q plots and comparison of density functions.

4 Results - 2SHLO in practice

We have been running our experiments in Python 3, relying on tick[16] library for Hawkes process simulations, AGD optimization and analytic tools. 2SHLO is also benefiting of well-known scipy scientific package.

To provide an idea of the computational performance of the fitting method, we mention that, on average, it takes 4.69 seconds to fit 2 dimensional exponential hawkes process with training length of T = 10000.

4.1 Validating 2SHLO performances over simulated data

To generate Hawkes processes, a few methods are available (as summarized in [2]), either based on thinning, time-change or cluster algorithm. We have here used tick functions for simulations which are based on the thinning algorithm described in (Ogata, 1981, p.25, Algorithm 2) [11].

To validate the performance of the fitting algorithm, we ran a series of 100 simulation and fitting procedure on the simulated 2D Hawkes time-stamps using different end times T to have different training sizes and using the following simulation parameters :

[TABLE]

Fig. 3 confirms that the estimated parameters converge to their true optimal values for increasing N number of ticks in the training set. This reassures that when the Hawkes process is well specified and indeed recovers the optimal parameters asymptotically.

4.2 Calibrating univariate Hawkes processes on ”mid-frequency” credit derivative trades

We have fitted a modified one-dimensional Hawkes process on the period from 15/01/2017 to 15/12/2017 of reported trades data333The trading data have been cleaned by discarding transactions that do not reflect market signals such as roll and switch trades on 3 indices ITXEB, ITXES and ITXEX. The modification consists of imposing a null value of the intensity in the gap between two days and the standard exponential kernel format during trading days so as to account for the absence of trading activity overnight. The model is thus slightly modified and described in the appendix section in a more detailed way. Before making this modeling choice, we experimented with a more elaborate model that accounts for the overnight spillover effect as described in the work of C. G.Bowsher [4] where the author defines the intensity as recursively depending on the level of the stochastic parts of the intensity at the end of the previous trading day and thus extending the effect of trading days over time. The fitting procedure resulted into a null parameter estimation of the spillover effect which motivated the choice of the model modification described above.

To validate the ’universality’ of our estimates, we tested the goodness of fit of the Hawkes model with the estimated parameters on out of sample period from 15/01/2018-15/07/2018. The model fitted well for the 3 indices time series (see Fig. 4) suggesting that each series has stable intrinsic parameters that characterize the pattern of trades arrival times.

4.3 Measuring the influence of traded volumes

When considering the distribution of trade volumes for credit indices, volume clusters are clearly outlined. For example, for index itxeb, three bins can be distinguished: ”small” trades, ”medium” trades and ”big” trades (See Fig. 5). Using multivariate Hawkes processes model, volume impacts have been modeled by considering 3 counting processes, splitting trades per volume size.

By using the same model as the last section to take overnight gaps into consideration, the fitting algorithm trained on 2018 data resulted in the following parameters estimates ($\mu$ and $\beta$ expressed in minute*-1*) :

[TABLE]

Visualizing in Fig. 5 the branching ratio which corresponds the the values of $\alpha$ allows some intuitive interpretation of the fitted model as it shows that large trades are purely self exciting process and they have a major influence in triggering small trades, we also notice that all trade categories have a non negligible self excitation component.

4.4 Measuring cross interactions between traded indices

We have adopted again a multivariate Hawkes processes with null intensity overnight model to study the cross interaction between indices trades times. The fitting algorithm trained on 2017 data resulted in the following parameters estimates ($\mu$ and $\beta$ expressed in minute*-1*) :

[TABLE]

From Fig. 6 we can notice that itxes index is the least influenced by the other indices with the very low cross excitation components and a high self excitation component which was expected due the its low liquidity. We also notice that the estimated baselines of all indices are greater for this multivariate model compared to the uni-dimensional model which can be explained by the fact that we decreased the ’Poissonnian’ behaviour of the arrival times by including more information about the influence of other indices.

5 Conclusion

Through point process analysis, Hawkes processes allows a relatively simple and interpretable representation of time dependant events with self and mutual excitements. With 2SHLO; a maximum likelihood based algorithm mixing convex and non-convex optimization, we have been able to fit fastly and properly exponential kernel multivariate Hawkes processes. The algorithm has been applied to credit derivative trading data, an unexplored ”mid-frequency” market for such processes. We have been able to emphasize self excitement for three index trades and to measure impacts of traded volumes. With the same framework, we were also able to quantify cross influences between indices. Next prospects turns towards testing 2SHLO performances on high dimensional datasets, for instance studying both price & volume trade influences as well as market bid/ask dynamics on credit derivatives in a forecasting perspective.

Appendix: Hawkes Processes Formulas

5.1 The Univariate Bowsher Hawkes process model

5.1.1 Definition:

As described in [4], the intensity of the univariate Bowsher process is defined recursively depending on the level of the stochastic parts of the intensity function at the end of the $(d-1)$th trading day and the contributions of the events occurring on day $d$. Thus it is necessary to perform the following data transformation: the real half-line is partitioned into intervals of length $l$ corresponding to the different trading days. This partition is written as :

[TABLE]

The model is thus defined by the (scalar) stochastic intensity :

[TABLE]

such that :

[TABLE]

where the parameters used are :

• •

$]\tau_{d-1},\tau_{d}]$ is the interval defining the day number d

• •

$\mu$ is the baseline intensity

• •

$\pi$ is the spillover adjacency coefficient

• •

$\rho$ is the spillover decay

• •

$\alpha$ is the self excitement adjacency coefficient

• •

$\beta$ is the self excitement decay

• •

$dN$ is the process differentiate

5.1.2 The Loglikelihood:

The minus loglikelihood of a univariate point process is defined as :

[TABLE]

Computing the minus loglikelihood in the Hawkes Bowsher model results into the following closed form :

[TABLE]

5.1.3 The Fitting procedure:

As the task of computing the gradient of each parameter of the loglikelihood is a fastidious one because of the recursive definition of the intensity, and also because we have only 5 parameters to fit $(\mu,\pi,\rho,\alpha,\beta)$, we choose to use the L-BFGS-B algorithm to perform the minimization of $\mathcal{L}_{T}(\lambda)$ which had delivered satisfying performance on fitting simulated data.

5.1.4 Data simulation:

We use the thinning algorithm to simulate artificial data described by the following algorithm :

5.2 The Day Gaps Multivariate Exponential Hawkes process model

5.2.1 Definition:

The intensity of the Day Gaps Multivariate Exponential Hawkes process is defined just like the standard multivariate exponential Hawkes process model described in [6] given by :

[TABLE]

Where :

• •

$D$ is the number of nodes

• •

$\mu_{i}$ are the baseline intensities

• •

$\phi_{ij}$ are the kernels

• •

$dN_{j}$ are the processes differentiates

with an exponential parametrisation of the kernels:

[TABLE]

• •

$\alpha_{ij}$ are the adjacency of the kernel

• •

$\beta_{ij}$ are the decays of the kernel

With the modification consisting on considering the intensity as a null function between two consecutive days to take the day gaps into consideration :

[TABLE]

5.2.2 The Loglikelihood:

The minus loglikelihood of a Multivariate Point Process is defined as :

[TABLE]

Computing the minus loglikelihood in The Day Gaps Multivariate Exponential Hawkes process model results into the following closed form :

[TABLE]

with $R_{mn}(i)$ are defined using the Ogoata [11] recursive formula :

[TABLE]

Where :

• •

$D$ is the number of days in the data

• •

$\delta$ is the length of one day

• •

$\mu_{i}$ are the baseline intensities

• •

$\alpha_{ij}$ are the adjacency coefficients

• •

$\beta_{ij}$ are the decays coefficients

• •

$\tau_{d}$ is the time defining the last time of the day number d

5.2.3 The Fitting procedure:

To fit this model, we use the algorithm described in this paper, using the projected Newton Descent for the convex optimization part and Nelder Mead for the non convex optimization step. The choice of the newton descent is justified by the easy computation of the gradient and hessian matrix which is sparse and because this algorithm has a quadratic convergence rate compared to gradient descent methods.

The gradient of $\mathcal{L}_{t}(\bm{\lambda})$ is thus given by the following close formulas :

[TABLE]

The hessian of $\mathcal{L}_{t}(\bm{\lambda})$ is given by the following close formulas :

[TABLE]

and all other coefficients are null. The fitting algorithm is described as follow:

And the projected newton descent is described as follow :

5.3 Some empirical results:

5.3.1 Loglikelihood landscape:

Visualizing the minus loglikelihood function for 2D exponential Hawkes process with regards to two non convex variables $(\beta_{0,0},\beta_{1,1})$ resulted into the following figure, emphasizing the flatness of the function to minimize and the existence of valleys that pose problem to gradient based optimization :

5.3.2 Goodness of fit size bins study:

Assessing the goodness of fit of the multivariate Hawkes process fitted on 2017 itxeb reported trades split into three series with regards to the size of trades:

5.3.3 Goodness of fit Multi index study:

Assessing the goodness of fit of the multivariate Hawkes process fitted on 2017 itxeb, itxes and itxex reported trades