Multi-kernel property in high-frequency price dynamics under Hawkes model

TL;DR

This paper introduces a multi-kernel Hawkes model to capture different responsive speeds of market participants in high-frequency trading, validated with empirical US stock data.

Contribution

It proposes a novel multi-kernel Hawkes framework for modeling high-frequency price dynamics and assesses the optimizer's effectiveness in complex likelihood landscapes.

Findings

Identification of UHF, VHF, and HF kernels in US stock data

Estimation of arrival times and contribution levels for each kernel

Validation of the model's ability to capture diverse market participant behaviors

Abstract

This study investigates and uses multi-kernel Hawkes models to describe a high-frequency mid-price process. Each kernel represents a different responsive speed of market participants. Using the conditional Hessian, we examine whether the numerical optimizer effectively finds the global maximum of the log-likelihood function under complicated modeling. Empirical studies that use stock prices in the US equity market show the existence of multi-kernels classified as ultra-high-frequency (UHF), very-high-frequency (VHF), and high-frequency (HF). We estimate the conditional expectations of arrival times and the degree of contribution to the high-frequency activities for each kernel.

Tables

Table 1. Table 1 : Summary statistics for estimates from the NBBO of IBM, 2018 under the one, two, and three-kernel Hawkes model for top to bottom

			UHF
stock	statistic	$\mu$	${\alpha }_s$	${\alpha }_c$	$\beta$
	mean	0.2768	322.8	128.6	883.5
IBM	median	0.2239	318.2	123.1	871.8
	SD	0.1595	79.86	34.68	191.4
			UHF			VHF
stock	statistic	$\mu$	${\alpha }_s$	${\alpha }_c$	$\beta$	${\alpha }_s$	${\alpha }_c$	$\beta$
	mean	0.2029	615.2	193.5	1915	2.909	4.605	35.75
IBM	median	0.1728	619.8	188.2	1922	2.786	4.344	34.47
	SD	0.1055	124.5	51.69	294.8	1.376	1.619	12.63
			UHF			VHF			HF
stock	statistic	$\mu$	${\alpha }_s$	${\alpha }_c$	$\beta$	${\alpha }_s$	${\alpha }_c$	$\beta$	${\alpha }_s$	${\alpha }_c$	$\beta$
	mean	0.1382	677.8	196.6	2191	6.010	8.397	78.13	0.0940	0.1271	1.449
IBM	median	0.1190	662.4	192.2	2132	4.862	7.457	66.77	0.0803	0.0661	1.180
	SD	0.0755	146.3	55.16	388.1	5.502	5.398	54.32	0.1065	0.2501	2.055

Table 2. Table 2 : Conditional expected time of the next arrival by each kernel, IBM, 2018

kernel	UHF
one	self	cross
	$319.4$$\mu s$	$145.7$$\mu s$
	UHF		VHF
two	self	cross	self	cross
	$132.1$$\mu s$	$59.09$$\mu s$	$2.207$$\mathrm{ms}$	$3.327$$\mathrm{ms}$
	UHF		VHF		HF
three	self	cross	self	cross	self	cross
	$115.8$$\mu s$	$55.61$$\mu s$	$1.032$$\mathrm{ms}$	$1.530$$\mathrm{ms}$	$72.74$$\mathrm{ms}$	$62.12$$\mathrm{ms}$

Table 3. Table 3 : The proportion of the causes of mid price changing events under two-kernel model and three-kernel model

two-kernel model
stock	year	base	UHF	VHF	HF
IBM	2016	38.40%	39.81%	21.79%
	2017	37.03%	43.96%	19.01%
	2018	35.52%	42.79%	21.69%
	2019	37.35%	42.41%	20.24%
AAPL	2016	35.45%	40.44%	24.11%
	2017	32.39%	43.10%	24.50%
	2018	31.34%	45.91%	22.76%
	2019	33.64%	43.31%	23.04%
three-kernel model
IBM	2016	26.54%	35.67%	20.29%	17.49%
	2017	24.54%	41.04%	18.51%	15.91%
	2018	24.22%	40.81%	19.70%	15.27%
	2019	25.17%	40.65%	16.91%	17.27%
AAPL	2016	26.24%	34.31%	16.53%	22.92%
	2017	18.70%	38.28%	20.90%	22.12%
	2018	16.30%	43.83%	20.61%	19.26%
	2019	16.12%	40.46%	20.42%	23.01%

Table 4. Table 4: Estimates and standard errors under the two-kernel model for January 2018 based on the NBBO of IBM

date	$\mu$	${\alpha }_{1s}$	${\alpha }_{1c}$	${\alpha }_{2s}$	${\alpha }_{2c}$	${\beta }_1$	${\beta }_2$
2018-01-02	0.1073	683.0	229.2	3.429	4.353	2041	46.34
	(0.0017)	(0.0596)	(0.1365)	(0.0347)	(0.0771)	(0.0500)	(0.1181)
2018-01-03	0.1697	637.0	215.3	4.773	6.473	1707	56.46
	(0.0021)	(0.0111)	(0.0392)	(0.0261)	(0.0464)	(0.0049)	(0.0352)
2018-01-04	0.1295	797.5	260.7	7.451	9.918	2423	81.54
	(0.0018)	(0.0850)	(0.0762)	(0.0421)	(0.0964)	(0.0386)	(0.1340)
2018-01-05	0.1193	804.3	200.2	4.314	5.815	2434	46.56
	(0.0018)	(0.1918)	(0.1414)	(0.2975)	(0.1434)	(0.2182)	(0.0872)
2018-01-08	0.1223	806.1	218.1	4.633	6.917	2430	51.03
	(0.0018)	(0.0605)	(0.0879)	(0.0394)	(0.0560)	(0.0153)	(0.0693)
2018-01-09	0.1129	707.7	179.5	3.519	3.869	2182	38.15
	(0.0017)	(0.0851)	(0.0237)	(0.1310)	(0.3077)	(0.0891)	(0.0997)
2018-01-10	0.0747	575.1	186.1	3.315	4.344	2006	40.49
	(0.0014)	(0.3245)	(0.3192)	(0.1709)	(0.2196)	(0.0404)	(0.1087)
2018-01-11	0.0918	737.8	200.5	3.537	5.074	2094	45.51
	(0.0015)	(0.0893)	(0.0980)	(0.0330)	(0.1472)	(0.0844)	(0.0306)
2018-01-12	0.1033	608.9	176.5	5.976	7.822	1755	64.47
	(0.0016)	(0.0746)	(0.1005)	(0.0481)	(0.2354)	(0.0312)	(0.0941)
2018-01-16	0.1778	484.1	186.6	4.066	6.259	1582	46.20
	(0.0021)	(0.0530)	(0.0377)	(0.0327)	(0.0576)	(0.0181)	(0.0488)
2018-01-17	0.1815	552.5	192.1	5.098	8.305	1777	60.78
	(0.0023)	(0.0381)	(0.1950)	(0.1334)	(0.1306)	(0.0661)	(0.1500)
2018-01-18	0.2074	660.1	207.4	5.144	5.520	1922	54.13
	(0.0023)	(0.0288)	(0.0242)	(0.0554)	(0.0470)	(0.0137)	(0.0209)
2018-01-19	0.2851	602.5	310.7	2.597	3.306	7.392	36.94
	(0.0027)	(0.2681)	(0.2784)	(0.2357)	(0.0352)	(0.2185 )	(0.0155)
2018-01-22	0.1443	662.8	290.4	2.333	4.713	1791	44.51
	(0.0019)	(0.0343)	(0.0358)	(0.0424)	(0.0164)	(0.0210)	(0.0648)
2018-01-23	0.1250	474.6	172.4	2.835	5.142	1657	38.61
	(0.0018)	(0.1237)	(0.1455)	(0.2019)	(0.1499)	(0.0533)	(0.1033)
2018-01-24	0.2082	576.3	178.0	2.506	6.525	2183	33.79
	(0.0024)	(0.0528)	(0.0775)	(0.0588)	(0.0717)	(0.0404)	(0.1154)
2018-01-25	0.1656	485.5	116.1	2.658	8.710	2130	36.02
	(0.0021)	(0.0863)	(0.0879)	(0.0960)	(0.0717)	(0.0537)	(0.0372)
2018-01-26	0.0873	585.7	140.9	4.259	7.836	2093	54.11
	(0.0015)	(0.1618)	(0.1657)	(0.0982)	(0.2668)	(0.1149)	(0.1861)
2018-01-29	0.1543	522.4	105.8	1.974	8.045	2209	31.83
	(0.0021)	(0.0347)	(0.1528)	(0.1201)	(0.1694)	(0.2561)	(0.0945)
2018-01-30	0.1676	566.0	138.7	2.614	7.365	2085	35.06
	(0.0021)	(0.0369)	(0.0625)	(0.0331)	(0.0574)	(0.0223)	(0.0281)
2018-01-31	0.1274	651.3	141.3	3.326	6.150	2073	36.84
	(0.0018)	(0.0635)	(0.0241)	(0.0723)	(0.0443)	(0.0900)	(0.0579)

Equations

$${\mathbf{N}}_t=\left[\begin{array}{c}{N}_1(t)\\ N_2(t)\end{array}\right],{\mathbf{\lambda }}_t=\left[\begin{array}{c}{\lambda }_1(t)\\ {\lambda }_2(t)\end{array}\right]$$

$${\mathbf{\lambda }}_t=\mathbf{\mu }+{\int }_{-\infty }^t\mathbf{h}(t-u)\mathrm{d}{\mathbf{N}}_u$$

$$\mathbf{h}=\sum _{k=1}^K{\mathbf{h}}_k$$

$${\mathbf{h}}_k(t)={\mathbf{\alpha }}_k\circ \left[\begin{array}{cc}{\mathrm{e}}^{-{\beta }_{k11}t}& {\mathrm{e}}^{-{\beta }_{k12}t}\\ {\mathrm{e}}^{-{\beta }_{k21}t}& {\mathrm{e}}^{-{\beta }_{k22}t}\end{array}\right].$$

$${\mathbf{\lambda }}_t=\mathbf{\mu }+\sum _{k=1}^K{\mathbf{\lambda }}_k(t)$$

$${\mathbf{\lambda }}_k(t)=\left[\begin{array}{c}{\lambda }_{k1}(t)\\ {\lambda }_{k2}(t)\end{array}\right]={\int }_{-\infty }^t{\mathbf{h}}_k(t-u)\mathrm{d}{\mathbf{N}}_u.$$

$${\beta }_{k1}:={\beta }_{k11}={\beta }_{k12},{\beta }_{k2}:={\beta }_{k21}={\beta }_{k22},$$

$${\mathbf{h}}_k(t)={\mathbf{\alpha }}_k\circ \left[\begin{array}{cc}{\mathrm{e}}^{-{\beta }_{k1}t}& {\mathrm{e}}^{-{\beta }_{k1}t}\\ {\mathrm{e}}^{-{\beta }_{k2}t}& {\mathrm{e}}^{-{\beta }_{k2}t}\end{array}\right].$$

$$\mathrm{d}{\mathbf{\lambda }}_t=\sum _{k=1}^K\mathrm{d}{\mathbf{\lambda }}_k(t),$$

$$\mathrm{d}{\mathbf{\lambda }}_k(t)=-{\mathbf{\beta }}_k{\mathbf{\lambda }}_k(t)\mathrm{d}t+{\mathbf{\alpha }}_k\mathrm{d}{\mathbf{N}}_t$$

$${\mathbf{\beta }}_k=\left[\begin{array}{cc}{\beta }_{k1}& 0\\ 0& {\beta }_{k2}\end{array}\right].$$

$$\mathbb{E}[{\mathbf{\lambda }}_t]={\left(\mathrm{I}-\sum _{k=1}^K{\mathbf{\beta }}_k^{-1}{\mathbf{\alpha }}_k\right)}^{-1}\mathbf{\mu },$$

$$\mathbb{E}[{\mathbf{\lambda }}_k(t)]={\mathbf{\beta }}_k^{-1}{\mathbf{\alpha }}_k\mathbb{E}[{\mathbf{\lambda }}_t].$$

$${\mathbf{\Lambda }}_t=\left[\begin{array}{c}{\mathbf{\lambda }}_1(t)\\ ⋮\\ {\mathbf{\lambda }}_K(t)\end{array}\right],\mathbf{\alpha }=\left[\begin{array}{c}{\mathbf{\alpha }}_1\\ ⋮\\ {\mathbf{\alpha }}_K\end{array}\right],\mathbf{\beta }=\mathrm{Dg}([{\beta }_{11}{\beta }_{12}\cdots {\beta }_{K1}{\beta }_{K2}])$$

$$\mathrm{d}{\mathbf{\Lambda }}_t=\mathbf{\beta }{\mathbf{\Lambda }}_t\mathrm{d}t+\mathbf{\alpha }\mathrm{d}{\mathbf{N}}_t.$$

$$\mathrm{d}{\mathbf{X}}_t={\mathbf{a}}_t\mathrm{d}t+{\mathbf{f}}_x(t)\mathrm{d}{\mathbf{N}}_t,\mathrm{d}{\mathbf{Y}}_t={\mathbf{b}}_t\mathrm{d}t+{\mathbf{f}}_y(t)\mathrm{d}{\mathbf{N}}_t.$$

$$\mathrm{d}[\mathbf{X}{\mathbf{Y}}^{\top }{]}_t={\mathbf{f}}_x(t)\mathrm{Dg}(\mathrm{d}{\mathbf{N}}_t){\mathbf{f}}_y^{\top }(t).$$

$$\mathrm{d}({\mathbf{X}}_t{\mathbf{Y}}_t^{\top })$$

$$\frac{\mathrm{d}\mathbb{E}[{\mathbf{X}}_t{\mathbf{Y}}_t^{\top }]}{\mathrm{d}t}=\mathbb{E}[{\mathbf{X}}_{t-}({\mathbf{b}}_t^{\top }+{\mathbf{\lambda }}_t^{\top }{\mathbf{f}}_y^{\top }(t))]+\mathbb{E}[({\mathbf{a}}_t+{\mathbf{f}}_x(t){\mathbf{\lambda }}_t){\mathbf{Y}}_{t-}^{\top }]+{\mathbf{f}}_x(t)\mathrm{Dg}(\mathbb{E}[{\mathbf{\lambda }}_t]){\mathbf{f}}_y^{\top }(t).$$

$$(\mathbf{\beta }-\mathbf{\alpha }\mathrm{J})\mathbb{E}[{\mathbf{\Lambda }}_t{\mathbf{\Lambda }}_t^{\top }]+\mathbb{E}[{\mathbf{\Lambda }}_t{\mathbf{\Lambda }}_t^{\top }](\mathbf{\beta }-{\mathrm{J}}^{\top }{\mathbf{\alpha }}^{\top })=\mathbf{\alpha }\mathbf{\mu }\mathbb{E}[{\mathbf{\Lambda }}_t^{\top }]+\mathbb{E}[{\mathbf{\Lambda }}_t]{\mathbf{\mu }}^{\top }{\mathbf{\alpha }}^{\top }+\mathbf{\alpha }\mathrm{Dg}(\mathbb{E}[{\mathbf{\lambda }}_t]){\mathbf{\alpha }}^{\top }$$

$$\mathrm{J}=\left[\begin{array}{ccccc}1& 0& \cdots & 1& 0\\ 0& 1& \cdots & 0& 1\end{array}\right].$$

$$(\mathrm{I}\otimes (\mathbf{\beta }-\mathbf{\alpha }\mathrm{J})+(\mathbf{\beta }-\mathbf{\alpha }\mathrm{J})\otimes \mathrm{I})\mathrm{vec}\left(\mathbb{E}[{\mathbf{\Lambda }}_t{\mathbf{\Lambda }}_t^{\top }]\right)=\mathrm{vec}\left(\mathbf{\alpha }\mathbf{\mu }\mathbb{E}[{\mathbf{\Lambda }}_t^{\top }]+\mathbb{E}[{\mathbf{\Lambda }}_t]{\mathbf{\mu }}^{\top }{\mathbf{\alpha }}^{\top }+\mathbf{\alpha }\mathrm{Dg}(\mathbb{E}[{\mathbf{\lambda }}_t]){\mathbf{\alpha }}^{\top }\right)$$

$$\mathbb{E}[{\mathbf{\lambda }}_t{\mathbf{N}}_t^{\top }]\approx \mathrm{A}t+\mathrm{B}$$

$$\mathrm{A}$$

$$\mathrm{B}$$

$${\mathbf{\lambda }}_k(t){\mathbf{\lambda }}_t^{\top }={\mathrm{J}}_k{\mathbf{\Lambda }}_t(\mathbf{\mu }+\mathrm{J}{\mathbf{\Lambda }}_t{)}^{\top }$$

$$[\bm{\mathrm{J}}_{k}]_{ij}=\left\{\begin{array}[]{@{}ll@{}}1,&\text{if}\ i=1\ \mathrm{and}\ j=2k-1\ \mathrm{or}\ i=2\ \mathrm{and}\ j=2k\\ 0,&\text{otherwise.}\end{array}\right.$$

$$\mathbb{E}[{\mathbf{N}}_t{\mathbf{N}}_t^{\top }]=\mathrm{A}{t}^2+\left(\mathrm{B}+{\mathrm{B}}^{\top }+\mathrm{Dg}(\mathbb{E}[{\mathbf{\lambda }}_t])\right)t$$

$$\mathrm{Var}(N_1(t)-N_2(t))=2\left[\begin{array}{cc}1& -1\end{array}\right]\mathrm{B}\left[\begin{array}{c}1\\ -1\end{array}\right]+\left[\begin{array}{cc}1& 1\end{array}\right]\mathbb{E}[{\mathbf{\lambda }}_t].$$

$${\mathbf{N}}_t=\left[\begin{array}{c}{N}_1(t)\\ ⋮\\ N_m(t)\end{array}\right],{\mathbf{\lambda }}_t=\left[\begin{array}{c}{\lambda }_1(t)\\ ⋮\\ {\lambda }_m(t)\end{array}\right]$$

$$h_{ij}(t)=\sum _{k=1}^K{\alpha }_{ijk}{\mathrm{e}}^{-{\beta }_{ijk}t}.$$

$$L(\mu ,\alpha ,\beta )=\sum _{i=1}^m\left(\sum _{n=1}^{N_i(T)}\log {\lambda }_i(t_{i,n})-{\int }_0^T{\lambda }_i(u)\mathrm{d}u\right)$$

Full text

1 Introduction

In recent years, interest in financial activities under the high-frequency level has increased as time records for financial market transactions, and quotes are achieved at higher time resolutions. Intraday events in major stock markets were recorded in milliseconds previously, and then the resolution gradually increases to microseconds and nanoseconds. High-frequency financial data comprise the actions of ultra-high-frequency and high-frequency traders as well as low-frequency participants in financial markets.

The Hawkes process (Hawkes, 1971), which has been applied to the natural and social sciences, captures the jump increment, called excitement, in intensities originating from past events and decays of the effect. Activities in financial markets under a tick structure, where prices vary tick-by-tick, such as transactions, quote revisions, and cancellations are activated by past events; their aftermath effects diminish over time. Due to this commonality, the Hawkes model has been actively used as a financial model in recent years. For the earliest studies, Bowsher (2007) proposed a generalized Hawkes model for multivariate financial market events, which incorporates various dependencies in market features, and Large (2007) utilized a mutually exciting Hawkes model to describe the resiliency in the limit order book.

Although the Hawkes model can be used in various fields of finance, such as credit risk analysis (Errais et al., 2010; Aït-Sahalia et al., 2015; Ma and Xu, 2016; Ketelbuters and Hainaut, 2022) and optimal execution (Choi et al., 2021; Da Fonseca and Malevergne, 2021; Jusselin, 2021; Gašperov and Kostanjčar, 2022), studies focusing on price dynamics, such as our study, are also abundant. To mention a few, Bacry et al. (2013a) used the multivariate Hawkes model for the tick dynamics of asset prices with the applications of microstructure noise analysis. Under settings similar to the standard Hawkes models, studies by Bacry and Muzy (2014); Da Fonseca and Zaatour (2014b, 2017); Arouri et al. (2019); Ji et al. (2020) contributed to the literature in aspects such as stylized facts of market impact, moment analysis, lead-lag relationship, capturing correlated jumps, extreme risk, and so on.

In addition, Hawkes models of a modified and extended form suitable for more complex structures have also been proposed. For example, to describe the bid and ask price dynamics, Zheng et al. (2014) introduced a multivariate constraint Hawkes-like point process and Lee and Seo (2022) proposed an extended Hawkes model to capture the spread dependent intensities. Jang et al. (2020) suggested the Hawkes flocking model to assess the systemic risk in crude oil and gasoline futures markets. Morariu-Patrichi and Pakkanen (2022) examined state-dependent Hawkes processes with an application to limit order book dynamics. Swishchuk (2021) proposed general regime-switching Hawkes models to describe the price processes in the limit order books. Zhang et al. (2021) suggested a new model by combining a Hawkes process and a finite-range contact process for stock price movements.

For reviews of other interesting studies on the Hawkes model with financial applications that are been mentioned here, please refer to Law and Viens (2015), Bacry et al. (2015), and Hawkes (2018).

We focus on the multi-kernel property in the high-frequency mid-price process in stock and futures markets. The mid-price, the mean of the best bid and ask prices, moves in multiples of the minimum tick in market microstructure. As mentioned before, many studies have modeled price movement in the tick structure using the Hawkes model. Among these, a single kernel is mainly used, as in Bacry et al. (2013b); Da Fonseca and Zaatour (2014a); Lee and Seo (2017b) or power-law kernel, as in Hardiman et al. (2013); Filimonov and Sornette (2015); Bacry et al. (2016). However, relatively few studies have been conducted on multi-kernel exponential models.

The intraday movements of the stock price are significantly influenced by ultra-high-frequency traders. However, there are also numerous activities that are performed in relatively lower frequencies. The responsiveness to previous events differs depending on their computing power, network speed, and available services such as colocation. The multi-kernel Hawkes model can capture the different responsive speeds and categorize the different types of activities in financial activities. We examine the differences in responsiveness among traders verified by multiple kernels.

As a basic property, we derive the first and second moments under the bi-variate multi-kernel Hawkes model, which is an extension of the results in a single kernel model, such as in Da Fonseca and Zaatour (2014b), Lee and Seo (2017b), and Cui et al. (2020). The moment expression is a closed-form formula and it helps verify the accuracy of the simulation of the multi-kernel Hawkes model.

We carefully examine the estimation procedure because the number of parameters increases as the number of kernels increases. We discuss the property of the global maximum of the log-likelihood based on a conditional Hessian to examine the reliability of the multivariate and multi-kernel model estimation. Various simulation examples are also presented.

The remainder of this paper is organized as follows. Section 2 proposes the multi-kernel Hawkes model for mid-price processes and discusses the global maximum property and moments. Section 3 present an empirical study based on high-frequency stock prices and the filtering effect. Finally, Section 4 concludes the paper.

2 Model

2.1 Multi-kernel model in high-frequency data

The mid-price processes of financial assets exhibit typical signs of high-frequency trading. The price reacts rapidly to previous events and the impact dissipates rapidly. However, market participants may have different processing speeds. We separate the kernels based on the frequency spectrum by introducing multi-kernels in the Hawkes model. It has often been argued that ultra-high-frequency financial raw data recorded at a high-resolution are not adequately suited to the single-kernel exponential Hawkes model. Note that stock movements have been recorded in nanoseconds in major stock exchanges since approximately 2015. Given that there are few in-depth studies on the financial application of the multi-kernel exponential Hawkes model, through this study, we intend to identify the advantages of the multi-kernel model as well as its disadvantages.

Consider a bi-variate counting process and its intensities:

[TABLE]

where $N_{1}(t)$ and $N_{2}(t)$ count the number of the up and down movements of a mid-price process up to $t$, respectively. Thus, the mid-price process can be represented as the difference between the two counting processes. In this study, bold face letters indicate vectors or matrices, as we deal with the bi-variate model. Let the vector of the intensity process be

[TABLE]

where $\bm{\mu}=[\mu_{1},\mu_{2}]^{\top}$ is constant.

This study focuses on the multi-kernel property, which is different from other studies on the single kernel model; that is

[TABLE]

where $k$ denotes the kernel number. Each kernel is

[TABLE]

where $\bm{\alpha}_{k}$ is a $2\times 2$ matrix whose element in $i$-th row and $j$-th column are represented by $\alpha_{ijk}.$ The smaller the $k$, the higher the kernel frequency, hence the larger the $\beta_{k}$.

As the number of parameters increases proportionally with the number of kernels, for the sake of model parsimony, we can assume that $\beta$s are the same for each kernel or for each raw of a kernel depending on the context. Multiple kernels are designed to capture different responsiveness speeds for high-frequency activities in the financial market.

2.2 Moment property

In this section, we computes the moments using the bi-variate multi-kernel Hawkes model. The intensity process of the multi-kernel Hawkes process can be rewritten as

[TABLE]

where the components are defined by

[TABLE]

To derive these formulas, we have additional restrictions on the Markov property:

[TABLE]

i.e.,

[TABLE]

Following this, the differential form of the intensity process can be represented by

[TABLE]

where

[TABLE]

For the simplicity of the formula, we assume that the processes are in a steady state at time 0, that is, by considering the intensity process to begin at $-\infty$, we always assume that the unconditional distribution of the intensity processes at time 0 are in a steady state. Under this assumption, the expectation of the intensity processes and their components are constant.

Proposition 1**.**

Under the steady state assumption, for $t>0$,

[TABLE]

To proceed further, define

[TABLE]

where $\mathrm{Dg}(\cdot)$ denotes a diagonal matrix whose diagonal entries are elements in the argument. Thus, $\Lambda_{t}$ is a $2K\times 1$ matrix, and $\bm{\alpha}$ is $2K\times 2$ and $\bm{\beta}$ is a diagonal $2K\times 2K$ matrix. Following this,

[TABLE]

Lemma 2**.**

Consider $2\times 1$ vector processes $\bm{X}$ and $\bm{Y}$ such that

[TABLE]

Subsequently, the differential form of the quadratic variation matrix process of $\bm{X}\bm{Y}^{\top}$ is represented by

[TABLE]

Thus,

[TABLE]

In addition,

[TABLE]

Proposition 3**.**

Under the steady state assumption, $\mathbb{E}[\bm{\Lambda}_{t}\bm{\Lambda}_{t}^{\top}]$ satisfies the following Sylvester equation:

[TABLE]

where $\bm{\mathrm{J}}$ is a $2\times 2K$ matrix such that

[TABLE]

Remark 4**.**

Solving Eq. (5) is equivalent to solving

[TABLE]

where $\otimes$ denotes the Kronecker product, and $\mathrm{vec}$ is the vectorization operator.

Proposition 5**.**

We have

[TABLE]

where

[TABLE]

Remark 6**.**

Note that

[TABLE]

where $\bm{\mathrm{J}}_{k}$ is a 2 by $2K$ matrix such that

[TABLE]

Proposition 7**.**

Under the steady state condition,

[TABLE]

where $\bm{\mathrm{A}}$ and $\bm{\mathrm{B}}$ are defined by Eqs. (6) and (7), respectively.

Proposition 8**.**

Under simple calculation, we have

[TABLE]

2.3 Conditional concavity

Multi-kernel Hawkes models are inferred by the maximum likelihood estimation (MLE). It is worthwhile to verify whether the numerical procedure of MLE is performed properly because of the complexity of the model. We discuss the issue on finding the global maximum of the log-likelihood function in terms of the conditional concavity introduced by Lee and Seo (2017a) in a single-kernel framework.

Consider a general case of the Hawkes model:

[TABLE]

where the intensities follow the $m$-dimensional version of Eq. (1). The component of $\bm{h}$ is

[TABLE]

The log-likelihood function is

[TABLE]

where $\{t_{i,n}\}$ are the $n$-th jump times of $N_{i}$. We rewrite

[TABLE]

where

[TABLE]

which is the remaining impact on $\lambda_{i}$ by the events which occur before time zero and this term can be ignored. If $\beta_{i}:=\beta_{ijk}$ for all $j$ and $k$, then,

[TABLE]

If the system is Markovian, then the assumption that $\varepsilon$ can be ignored is unnecessary.

By integrating Eq. (10),

[TABLE]

and by substituting the above in Eq. (8),

[TABLE]

where

[TABLE]

and we use

[TABLE]

We examine the conditional concavity of $L(\mu,\alpha|\beta)$ with fixed $\beta$s with the conditional Hessian. Let

[TABLE]

then

[TABLE]

Since the second derivatives of $\mu_{i}T$ with respect to $\mu_{i}$s and $\alpha_{ij}$s are all zeros, it is sufficient to check the second derivatives of $L_{i,n}(\mu,\alpha|\beta)$ to examine the Hessian of $L(\mu,\alpha|\beta)$. Note that

[TABLE]

As the Hessian of $L_{i,n}(\mu,\alpha|\beta)$ is negative semidefinite, the Hessian of $L(\mu,\alpha|\beta)$ is and hence $L(\mu,\alpha|\beta)$ is concave. Due to the conditional concavity, a numerical optimizer can find

[TABLE]

for each $\beta$ (if it exists). Therefore, once we determine $L^{*}(\beta)$ for enough $\beta$s and if we can calculate the maximum of $L^{*}(\beta)$, it is possible to determine the global maximum of $L(\mu,\alpha,\beta)$. Given that this strategy takes a long time, it is only used in this study when examining the global maximum of the Hawkes log-likelihood. For the empirical studies, a typical quasi-Newton method is applied under the assumption that the numerical optimizer works well.

2.4 Simulation study

First, by comparing the theoretically calculated moment values derived in Subsection 2.2 with the sample mean obtained through simulation, we verify whether the simulation process is performed well. The simulation method is based on the thinning algorithm, which can generally be applied to nonhomogeneous Poisson processes (Lewis and Shedler, 1979).

Example 1**.**

Under the three kernel model, with presumed parameters

[TABLE]

[TABLE]

and $t=1,000$, by Propositions 1 and 7, we have

[TABLE]

Our simulation with $10,000$ generated paths shows that the sample means

[TABLE]

which confirms that the simulation process works well.

Next, we examine the estimation performance of the model. The concavity of the log-likelihood depends on the sample size and the branching ratio. If the number of observations is insufficient and the branching ratio or $\alpha$ is close to zero, then the log-likelihood function is probably not concave. The log-likelihood function is likely concave with sufficient sample size and significant $\alpha$, Note that the Hawkes process with insignificant $\alpha$ compared with $\beta$ is very close to the homogeneous Poisson; it is evident that the model with insignificant $\alpha$ is not well estimated. The simulation study shows that the more significant $\alpha$, the lower the sample size that is needed and if $\alpha$ is close to zero, the estimation could be incorrect.

Figure 1 illustrates examples of non-concave log-likelihood functions in a univariate Hawkes model. The left represents the small sample size case with parameters $\alpha=0.15,\beta=1.0$ and 10 observations, while the right represents the case of insignificant $\alpha=0.001$ with $\beta=1.0$ and sample size 500. Using a numerical optimizer, we determine the $L^{*}(\beta)$ for each $\beta$ on a sufficiently dense set. As non-concavities are observed in the figure, the numerical optimizer may fail to find the global maximum.

A small $\alpha$ value can cause an inaccurate estimation for $\beta$. In a multivariate Hawkes model, if $\beta$s are considered insignificantly different from each other, it would be preferable to assume that $\beta$s are equal in the model for parsimony. For example, consider the following model:

[TABLE]

where $\alpha_{12}=0.001$ is close to zero, and $\beta_{ij}$s differ insignificantly. With a simulated path with 5,000 observations, the estimates of $\beta$s are

[TABLE]

where $\hat{\beta}_{12}$ is inaccurate. With a constraint of all equal $\beta_{ij}$, the estimate of $\beta_{11}=\beta_{12}=\beta_{21}=\beta_{22}$ is 2.3048, which is much more reliable.

We conduct additional experiments to test the relationship between the sample size and the uniqueness of the local maximum. For computational complexity, the experiment is based on a one-dimensional Hawkes model. However, it is important to demonstrate that the sample size is crucial.

A path is simulated for each of the 0.1, 0.5, and 0.9 branching ratios. A numerical optimizer finds $L^{*}(\beta)$ for all $\beta$s over a sufficiently dense set. Examining computed $L^{*}(\beta)$ makes it possible to check the uniqueness of local maximum. We assume that the numerical optimizer would ultimately succeed in finding the global maximum of $L$ and count it a success, if there is a unique local maximum in $L^{*}$. Otherwise, we cannot be sure of its success and count this as a failure. Figure 2 shows success rates after repeating 100 experiments. As mentioned earlier, the lowest success rate is observed when the branching ratio is low. However, with more than 150 samples, the rate is greater than 90%. We use more than 10,000 observations in our empirical study. Furthermore, as the branching ratios are generally significant, we assume that the numerical results are reliable.

3 Empirical result

3.1 Basic result

We conduct empirical studies using high-frequency stock price data under the US equity and futures market from 2016222From 2016, the data provided by the New York Stock Exchange is recorded in nanoseconds. to 2019. The data comprise of the timestamps of changes in the national best bid and offer (NBBO) prices, derived from the trade and quote (TAQ) data released by the New York Stock Exchange.

Before conducting the main analysis, we visualize $L^{*}(\beta_{1},\beta_{2})$ under the two-kernel model for IBM on January 3, 2017 as shown in Figure 3. The global maximum is located at $\beta_{1}=880$ and $\beta_{2}=12$. The other estimates at the maximum are $\hat{\mu}_{1}=\hat{\mu}_{2}=0.089$, $\hat{\alpha}_{1s}=258.8$, $\hat{\alpha}_{1c}=102.6$, $\hat{\alpha}_{2s}=0.926$ and $\hat{\alpha}_{2c}=1.314$. The surface is not entirely concave (see the region where $\beta_{1}$ is close to zero in the figure). However, the function is concave and regular around the global maximum.

Now, we discuss the model-selection issue. The MLEs are performed on IBM stock’s mid-price process with one-, two-, three-, and four-kernel Hawkes models. In Figure 4, based on IBM data on Jan 03, 2018 we present a Q-Q (quantile-quantile) plot of the residuals versus the unit exponential distribution. The set of residuals is defined by

[TABLE]

where $\hat{\lambda}_{i}$ denote inferred intensities. The points are closer to the straight line, as the number of kernels increases up to three. Meanwhile, the three-kernel and the four-kernel models differ insignificantly.

We obtain similar results based on the Akaike information criterion (AIC). Figure 5 plots the daily AICs for various kernel models based on IBM data for Jan 2018. The solid black line, the blue dotted line and the red dash-dot line represent the one-kernel, two-kernel, and three-kernel models, respectively. The one-kernel model has a significant information loss compared to the others. The three-kernel model is slightly better than the two-kernel model. Although not depicted, the four-kernel model has very similar AICs to those of the three-kernel model. We observe similar patterns in the other cases. The model extensions are useful for up to two or three kernels but not for higher.

The estimates with a one second time unit under the two-kernel model are reported in Table LABEL:Table:4kernel, with standard errors in parentheses. The estimates of $\alpha$s of the first kernel are several hundred; $\hat{\beta}$s range from approximately 1,000 to 2,000, $\hat{\alpha}$s of the second kernel are less than 10, and $\hat{\beta}$s range from approximately 30 to 80.

Figure 6 and 7 illustrate the dynamics of the daily estimates under the two-kernel and three-kernel models, respectively. We use the NBBO mid-price processes of IBM from 2016 to 2018. We term each kernel of the two-kernel model as ultra-high-frequency (UHF) and very-high-frequency (VHF) kernels. Figure 6(a) reveals $\hat{\alpha}_{1s}$, $\hat{\alpha}_{1c}$ and $\hat{\beta}_{1}$ for the UHF and $\hat{\alpha}_{1s}$s range from approximately 100 to 200 whereas $\hat{\alpha}_{1c}$s range from approximately 200 to 1,000.

This implies that, if the mid-price changes, the intensities increase instantly at several hundred. This can lead to numerous additional activities at frequencies of hundreds per second. Moreover, as shown in the figure, $\hat{\beta}_{1}$ is quite large and the effect dissipates quickly. This fast responsiveness is possibly due to automated high-frequency trading and quotes. Meanwhile, both $\hat{\alpha}_{2s}$ and $\hat{\alpha}_{2c}$ for the VHF kernel ranged from approximately 0 to 15, which is much less than $\hat{\alpha}$s in the UHF kernel. Furthermore, $\hat{\beta}_{2}$ is distributed between 0 and 80. This implies a longer persistence than in the UHF kernel.

For the three-kernel model, we term each kernel to be ultra-high-frequency (UHF), very-high-frequency (VHF), and high-frequency (HF). The $\hat{\alpha}$s in the UHF kernel are hundreds, $\hat{\alpha}$s in the VHF kernel are up to 60 and $\hat{\alpha}$s in the HF kernel are less than 2.5. Table 1 presents the summary statistics of the parameter estimates for various models.

3.2 Responsiveness

Under the Hawkes model, the intensities are instantly excited when an event occurs and decay to the level they would have been at if the excitement had not occurred. The increased intensity may or may not cause additional event during this period. If the excitement is large, and the decay is slow, then the occurrence of an event is highly possible. Based on this argument, we compute the conditional expectation of the arrival time of an event caused by a single excitement. This analysis helps us understand the responsive speed of each kernel.

If the excitation parameter is $\alpha$ and the exponential decaying parameter is $\beta$, then the remaining amount of excitement in the intensity owing to $\alpha$ when time $s$ has elapsed since the event occurred is $\alpha\exp(-\beta s).$ Based on the survival analysis with a hazard function $\alpha\exp(-\beta s)$, the conditional probability for the inter-arrival time, $\tau$, of the next event owing to excitement $\alpha$ is

[TABLE]

Further, the conditional expectation of the next event owing to excitement $\alpha$ is

[TABLE]

which can be computed numerically. Using the above formula and the median estimates of $\alpha$s and $\beta$s for IBM in 2018, Table 2 presents the estimated conditional expectation of the arrival time for each kernel.

In the three-kernel model, the conditional expectation of the arrival time originating from the UHF kernel ranges from approximately 50 to 100$\mathrm{\SIUnitSymbolMicro s}$. The response time is very fast, and the corresponding activities may include automated actions such as flickering quotes and fleeting orders (Hasbrouck and Saar, 2009). The conditional expectation of arrival time by the VHF kernel ranges from approximately 1-2$\mathrm{ms}$, while that for the HF kernel, ranges from approximately 60-70 $\mathrm{ms}$. The expected arrival time by the base intensity is approximately 100-200$\mathrm{m}\mathrm{s}$.

Next, we examine the proportions of the causes of events among the kernels and the base intensity. According to the immigrant offspring argument, an event can be an external immigrant represented by the base intensity, or the descendant of an internally generated offspring described by the kernels. Based on the inferred intensities derived from the maximum likelihood estimates, we compute the probability that the parent of an event belongs to the source, that is, kernels or the base intensity. By comparing the estimated kernel components of the intensities at the time of an event, we calculate the empirical probability of the causes as the proportion of the intensity component. More precisely, in the three-kernel model, assuming that an upward move occurs at time $t$, the following

[TABLE]

represent the estimated probabilities that the cause of event belongs to base intensity, UHF, VHF, and UF kernels, respectively.

The annually aggregated empirical probabilities are presented in Table 3 from 2016 to 2019 for IBM and AAPL. Approximately 40% of the events are caused by the UHF kernel which is the largest proportion. The other kernels have similar percentages, ranging from approximately 15-25%. The percentages of the base intensity range from approximately 30-40% and 15-25% for the two-kernel and three-kernel models, respectively.

4 Conclusion

The multi-kernel Hawkes model is applied to the high-frequency price process and the basic moment properties are derived. Furthermore, the suitability of the MLE is examined in terms of conditional concavity. The empirical results reveal that high-frequency mid-price dynamics derived directly from raw data follow multi-kernel Hawkes processes. The diagnostics test show that two- or three-kernel model are appropriate for ultra-high-frequency modeling.

In the two-kernel model, the kernel is composed of the ultra-high-frequency (50-100$\mathrm{\SIUnitSymbolMicro s}$ in expected response time) and very-high-frequency (2-3$\mathrm{ms}$). In the three-kernel model, the kernel is composed of the ultra-high-frequency (50-100$\mathrm{\SIUnitSymbolMicro s}$), very-high-frequency (1-1.5$\mathrm{ms}$), and high-frequency (60-70$\mathrm{ms}$). The percentage estimated by the ultra-high-frequency kernel is approximately 40%, which is the largest proportion in both the two and three-kernel models. This percentage is approximately 20% for the very-high-frequency and high-frequency kernels. The estimated percentages of external events are 15-25% and 30-40% for the three-kernel and two-kernel models, respectively.

Acknowledgements

This work has supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT)(No. NRF-2021R1C1C1007692).

Appendix A Proofs

Proof for Proposition 1.

Eq. (4) implies

[TABLE]

Under the steady state condition, by taking expectation to the both sides of the above equation, we have

[TABLE]

and using $\mathbb{E}[\bm{\lambda}_{k}(t)]=\mathbb{E}[\bm{\lambda}_{k}(0)]$,

[TABLE]

and substituting the above to

[TABLE]

which is derived by taking the expectation on Eq. (2), we have the desired result. ∎

Proof for Propositoin 3.

By Lemma 2,

[TABLE]

Since

[TABLE]

we have

[TABLE]

∎

Proof for Propositon 5.

For simplicity, consider only the particular solution for $\frac{\mathrm{d}\mathbb{E}[\bm{\lambda}_{t}\bm{N}_{t}^{\top}]}{\mathrm{d}t}$. Assume that

[TABLE]

Since

[TABLE]

we have

[TABLE]

and

[TABLE]

Note that

[TABLE]

and

[TABLE]

and

[TABLE]

and

[TABLE]

Substituting the above to Eq. (13), we obtain the desired result. ∎

Proof for Propositon 7.

Using Lemma 2 with $\bm{a}=\bm{b}=0,\bm{f}_{x}=\bm{f}_{y}=\bm{\mathrm{I}}$ and Proposition 5, we obtain the result. ∎

Appendix B Estimation results