Performance of information criteria used for model selection of Hawkes process models of financial data

TL;DR

This study evaluates the effectiveness of AIC, BIC, and HQ information criteria in selecting the correct order of Hawkes process models for high-frequency financial data, considering sample size effects.

Contribution

It provides a systematic comparison of IC performance in model selection for Hawkes processes using simulated data with realistic parameters.

Findings

BIC outperforms AIC and HQ in correct model selection.

Sample size significantly influences IC accuracy.

Results align with theoretical predictions and empirical data.

Abstract

We test three common information criteria (IC) for selecting the order of a Hawkes process with an intensity kernel that can be expressed as a mixture of exponential terms. These processes find application in high-frequency financial data modelling. The information criteria are Akaike's information criterion (AIC), the Bayesian information criterion (BIC) and the Hannan-Quinn criterion (HQ). Since we work with simulated data, we are able to measure the performance of model selection by the success rate of the IC in selecting the model that was used to generate the data. In particular, we are interested in the relation between correct model selection and underlying sample size. The analysis includes realistic sample sizes and parameter sets from recent literature where parameters were estimated using empirical financial intra-day data. We compare our results to theoretical predictions and similar empirical findings on the asymptotic distribution of model selection for consistent and inconsistent IC.

Tables

Table 1. Table 1: Simulation parameters (Parameter Set 1)

	$\mu$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\beta }_1$	${\beta }_2$	${\beta }_3$
P=1	0.5	9	–	–	10	–	–
P=2	0.5	0.00066	100	–	0.001	300	–
P=3	0.5	0.00033	3.3	100	0.001	10	300

Table 2. Table 2: Parameter set taken from [ 26 ] (Parameter Set 2)

$\mu$	${\alpha }_1$	${\alpha }_2$	${\beta }_1$	${\beta }_2$
0.05	0.01761905	0.28	0.04761905	0.6666667

Table 3. Table 3: Absolute RMSE values for MLE of the exponential Hawkes models of order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} using Parameter Set 1 with varying time horizons T ∈ { 500 , 1000 , 2000 , 5000 } 𝑇 500 1000 2000 5000 T\in\{500,1000,2000,5000\} . The order of the model which was used for simulation coincides with the model used for fitting. Thus, the true parameter values are known and the RMSE is expected to decrease as the MLE improves.

		$\mu$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\beta }_1$	${\beta }_2$	${\beta }_3$	Average sample size
P=1	T=500	0.039664	0.42256	–	–	0.44731	–	–	2483
	T=1000	0.02763	0.32473	–	–	0.32928	-	–	5019
	T=2000	0.018738	0.21756	–	–	0.22197	-	–	9977
	T=5000	0.011276	0.13846	–	–	0.14544	-	–	24962
P=2	T=500	0.071796	6.0214	12.399	–	37.322	44.92	–	470
	T=1000	0.061258	0.00085434	7.9865	–	0.0023803	18.956	–	1121
	T=2000	0.049989	0.00020347	4.7732	–	0.00045077	11.415	–	2977
	T=5000	0.042938	0.00010551	2.3918	–	0.00018339	5.9634	–	12883
P=3	T=500	0.07713	0.3232	1.3408	9.8602	1.5085	9.7124	30.678	929
	T=1000	0.061804	0.00036118	0.29469	6.2401	0.0021189	0.76784	18.334	2207
	T=2000	0.051527	0.0001279	0.19655	3.8663	0.00058096	0.49288	11.77	5840
	T=5000	0.045562	0.000055755	0.10766	1.8741	0.00019317	0.27921	5.6723	25017

Table 4. Table 4: Relative RMSE values for MLE of the exponential Hawkes models of order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} using Parameter Set 1 with varying time horizons T ∈ { 500 , 1000 , 2000 , 5000 } 𝑇 500 1000 2000 5000 T\in\{500,1000,2000,5000\} . The order of the model which was used for simulation coincides with the model used for fitting. Thus, the true parameter values are known and the RMSE is expected to decrease as the MLE improves. The values are given in percent.

		$\mu$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\beta }_1$	${\beta }_2$	${\beta }_3$	Average sample size
P=1	T=500	7.9328	4.6951	–	–	4.4731	–	–	2483
	T=1000	5.526	3.6082	–	–	3.2928	–	–	5019
	T=2000	3.7475	2.4173	–	–	2.2197	–	–	9977
	T=5000	2.2551	1.5384	–	–	1.4544	–	–	24962
P=2	T=500	14.359	912330	12.399	–	3732200	14.973	–	470
	T=1000	12.252	129.45	7.9865	–	238.03	6.3186	–	1121
	T=2000	9.9978	30.829	4.7732	–	45.077	3.805	–	2977
	T=5000	8.5877	15.986	2.3918	–	18.339	1.9878	–	12883
P=3	T=500	15.426	97939	40.63	9.8602	150850	97.124	10.226	929
	T=1000	12.361	109.45	8.9301	6.2401	211.89	7.6784	6.1113	2207
	T=2000	10.305	38.758	5.956	3.8663	58.096	4.9288	3.9233	5840
	T=5000	9.1124	16.895	3.2624	1.8741	19.317	2.7921	1.8908	25017

Table 5. Table 5: Model selection for simulated data of an exponential Hawkes model of order P=1 using Parameter Set 1 with varying time horizons T ∈ { 500 , 1000 , 2000 , 5000 } 𝑇 500 1000 2000 5000 T\in\{500,1000,2000,5000\} . The numbers indicate how often the model order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} is selected among the 1000 1000 1000 samples and are given in percent. Bold numbers show which model was selected most often.

	Time horizon	P=1	P=2	P=3	Average sample size
AIC	T=500	92.8	6.9	0.3	2483
	T=1000	91.6	7.9	0.5	5019
	T=2000	92.1	7.6	0.3	9977
	T=5000	93.7	6.1	0.2	24962
BIC	T=500	99.8	0.2	0	2483
	T=1000	100	0	0	5019
	T=2000	100	0	0	9977
	T=5000	100	0	0	24962
HQ	T=500	98.9	1.1	0	2483
	T=1000	98.6	1.2	0.2	5019
	T=2000	99.2	0.8	0	9977
	T=5000	99.7	0.3	0	24962

Table 6. Table 6: Model selection for simulated data of an exponential Hawkes model of order P=2 using Parameter Set 1 with varying time horizons T ∈ { 500 , 1000 , 2000 , 5000 } 𝑇 500 1000 2000 5000 T\in\{500,1000,2000,5000\} . The numbers indicate how often the model order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} is selected among the 1000 1000 1000 samples and are given in percent. Bold numbers show which model was selected most often.

	Time horizon	P=1	P=2	P=3	Average sample size
AIC	T=500	48.2	50.3	1.5	470
	T=1000	0.2	99	0.8	1121
	T=2000	0	96.9	3.1	2977
	T=5000	0	93.7	6.3	12883
BIC	T=500	94.5	5.4	0.1	470
	T=1000	10.3	89.7	0	1121
	T=2000	0	100	0	2977
	T=5000	0	100	0	12883
HQ	T=500	76.9	22.9	0.2	470
	T=1000	2.1	97.8	0.1	1121
	T=2000	0	99.4	0.6	2977
	T=5000	0	99.6	0.4	12883

Table 7. Table 7: Model selection for simulated data of an exponential Hawkes model of order P=3 using Parameter Set 1 with varying time horizons T ∈ { 500 , 1000 , 2000 , 5000 } 𝑇 500 1000 2000 5000 T\in\{500,1000,2000,5000\} . The numbers indicate how often the model order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} is selected among the 1000 1000 1000 samples and are given in percent. Bold numbers show which model was selected most often.

	Time horizon	P=1	P=2	P=3	Average sample size
AIC	T=500	0	53.7	46.3	929
	T=1000	0	0.2	99.8	2207
	T=2000	0	0	100	5840
	T=5000	0	0	100	25017
BIC	T=500	0	96.5	3.5	929
	T=1000	0	25	75	2207
	T=2000	0	0	100	5840
	T=5000	0	0	100	25017
HQ	T=500	0	81.6	18.4	929
	T=1000	0	4.7	95.3	2207
	T=2000	0	0	100	5840
	T=5000	0	0	100	25017

Table 8. Table 8: Absolute RMSE values for MLE of the exponential Hawkes models of order P = 2 𝑃 2 P=2 using Parameter Set 2 with varying time horizons T ∈ { 600 , 900 , 1800 , 3600 , 7200 , 21600 } 𝑇 600 900 1800 3600 7200 21600 T\in\{600,900,1800,3600,7200,21600\} . The order of the model which was used for simulation coincides with the model used for fitting. Thus, the true parameter values are known and the RMSE is expected to decrease as the MLE improves.

	$\mu$	${\alpha }_1$	${\alpha }_2$	${\beta }_1$	${\beta }_2$	Average sample size
T=600	0.030176	0.076208	12489000	0.17692	27063000000	135
T=900	0.025411	0.052971	1704000	0.11054	15861000000	205
T=1800	0.016031	0.023916	0.078183	0.049322	4.0458	417
T=3600	0.010572	0.0095225	0.039458	0.020001	0.15516	853
T=7200	0.0070848	0.0055058	0.025844	0.011919	0.088505	1708
T=21600	0.0039548	0.0030036	0.014737	0.006448	0.051022	5144

Table 9. Table 9: Relative RMSE values for MLE of the exponential Hawkes models of order P = 2 𝑃 2 P=2 using Parameter Set 2 with varying time horizons T ∈ { 600 , 900 , 1800 , 3600 , 7200 , 21600 } 𝑇 600 900 1800 3600 7200 21600 T\in\{600,900,1800,3600,7200,21600\} . The order of the model which was used for simulation coincides with the model used for fitting. Thus, the true parameter values are known and the RMSE is expected to decrease as the MLE improves. The values are given in percent.

	$\mu$	${\alpha }_1$	${\alpha }_2$	${\beta }_1$	${\beta }_2$	Average sample size
T=600	60.353	432.53	4460300000	371.53	4059400000000	135
T=900	50.822	300.64	608590000	232.14	2379200000000	205
T=1800	32.061	135.74	27.923	103.58	606.87	417
T=3600	21.144	54.047	14.092	42.003	23.275	853
T=7200	14.17	31.249	9.2299	25.03	13.276	1708
T=21600	7.9096	17.047	5.263	13.541	7.6533	5144

Table 10. Table 10: Model selection for simulated data of an exponential Hawkes model of order P=2 using Parameter Set 2 with varying time horizons T ∈ { 600 , 900 , 1800 , 3600 , 7200 , 21600 } 𝑇 600 900 1800 3600 7200 21600 T\in\{600,900,1800,3600,7200,21600\} . The numbers indicate how often the model order P ∈ { 1 , 2 , 3 } 𝑃 1 2 3 P\in\{1,2,3\} is selected among the 1000 1000 1000 samples and are given in percent. Bold numbers show which model was selected most often.

	Time horizon	P=1	P=2	P=3	Average sample size
AICc/AIC	T=600	51.4	48	0.6	135
	T=900	35.1	63.8	1.1	205
	T=1800	6.7	90.2	3.1	417
	T=3600	0	96.9	3.1	853
	T=7200	0	94.7	5.3	1708
	T=21600	0	94.1	5.9	5144
AIC	T=600	49.3	50.1	0.6	135
	T=900	34.8	64.1	1.1	205
	T=1800	6.7	90.2	3.1	417
	T=3600	0	96.9	3.1	853
	T=7200	0	94.7	5.3	1708
	T=21600	0	94.1	5.9	5144
BIC	T=600	86.5	13.5	0	135
	T=900	79.8	20.2	0	205
	T=1800	42.7	57.2	0.1	417
	T=3600	5	95	0	853
	T=7200	0.1	99.9	0	1708
	T=21600	0	100	0	5144
HQ	T=600	69	30.8	0.2	135
	T=900	55.6	44.4	0	205
	T=1800	17.5	81.8	0.7	417
	T=3600	1	98.7	0.3	853
	T=7200	0	98.9	1.1	1708
	T=21600	0	99.2	0.8	5144

Equations

$$\lambda (t|{\mathcal{F}}_t):=\underset{\Delta t\to 0}{\lim }\mathbb{P}(N(t+\Delta t)-N(t)=1|{\mathcal{F}}_t)/\Delta t,$$

$$\lambda (t)=\mu +{\int }_0^{t}g(t-\tau )\mathrm{d}N(\tau ),$$

$$g(t)=\sum _{m=1}^P{\alpha }_m{\mathrm{e}}^{-{\beta }_{m}t}.$$

$$\lambda (t)=\mu +\sum _{m=1}^P{\alpha }_m\sum _{i=1}^k{\mathrm{e}}^{-{\beta }_m(t-t_i)}$$

$$\Lambda :=\mathbb{E}[\lambda (t)]=\frac{\mu }{1-{\int }_0^{\infty }g(\nu )\mathrm{d}\nu },$$

$$\phi (t)=\mu +\sum _{m=1}^P{\int }_0^t{\alpha }_m{\mathrm{e}}^{-{\beta }_m(t-u)}\phi (u)\mathrm{d}u.$$

$$\tilde{\phi }(s)$$

$$={\int }_0^{\infty }{\mathrm{e}}^{-st}\mu \mathrm{d}t+\sum _{m=1}^P{\alpha }_m{\int }_{t=0}^{\infty }{\mathrm{e}}^{-st}{\int }_{u=0}^t{\mathrm{e}}^{-{\beta }_m(t-u)}\phi (u)\mathrm{d}u\mathrm{d}t$$

$$=\frac{\mu }{s}+\sum _{m=1}^P{\alpha }_m{\int }_{u=0}^{\infty }{\mathrm{e}}^{-su}\phi (u){\int }_{t=u}^{\infty }{\mathrm{e}}^{-(s+{\beta }_m)(t-u)}\mathrm{d}t\mathrm{d}u$$

$$=\frac{\mu }{s}+\sum _{m=1}^P\frac{{\alpha }_m}{s+{\beta }_m}{\int }_{u=0}^{\infty }{\mathrm{e}}^{-su}\phi (u)\mathrm{d}u=\frac{\mu }{s}+\left(\sum _{m=1}^P\frac{{\alpha }_m}{s+{\beta }_m}\right)\tilde{\phi }(s),$$

$$\tilde{\phi }(s)=\frac{\frac{\mu }{s}}{1-{\sum }_{m=1}^P\frac{{\alpha }_m}{s+{\beta }_m}}.$$

$$\tilde{\phi }(s)=\frac{\mu }{s}\frac{{\prod }_{m=1}^P(s+{\beta }_m)}{{\prod }_{m=1}^P(s+{\beta }_m)-{\sum }_{m=1}^P{\alpha }_m{\prod }_{k\ne m}(s+{\beta }_k)}.$$

$$\tilde{\phi }(s)=\frac{\mu (s+{\beta }_1)}{s(s+{\beta }_1-{\alpha }_1)}=\frac{\mu }{{\beta }_1-{\alpha }_1}\left(\frac{{\beta }_1}{s}-\frac{{\alpha }_1}{s+{\beta }_1-{\alpha }_1}\right),$$

$$\phi (t)=\frac{\mu }{{\beta }_1-{\alpha }_1}\left({\beta }_1-{\alpha }_1{\mathrm{e}}^{-({\beta }_1-{\alpha }_1)t}\right),t>0.$$

$$\tilde{\phi }(s)=\frac{\mu (s+{\beta }_1)(s+{\beta }_2)}{s[(s+{\beta }_1)(s+{\beta }_2)-{\alpha }_1(s+{\beta }_2)-{\alpha }_2(s+{\beta }_1)]}$$

$$\tilde{\phi }(s)=\frac{P(s)}{Q(s)}=\sum _{i=1}^3\frac{P(s_i)}{Q^{\prime }(s_i)(s-s_i)}=\mu \left(\frac{A_1}{s}+\frac{A_2}{s-s_2}+\frac{A_3}{s-s_3}\right),$$

$$s_1=0,s_2=\frac{1}{2}(\gamma -\xi ),s_3=\frac{1}{2}(\gamma +\xi )$$

$$\text{with}\gamma ={\alpha }_1+{\alpha }_2-{\beta }_1-{\beta }_2\text{and}\xi =\sqrt{{\gamma }^2-4({\beta }_1{\beta }_2-{\alpha }_1{\beta }_2-{\alpha }_2{\beta }_1)}.$$

$$A_1(s-s_2)(s-s_3)+A_{2}s(s-s_3)+A_{3}s(s-s_2)\stackrel{!}{=}(s+{\beta }_1)(s+{\beta }_2)$$

$$A_1+A_2+A_3$$

$$-A_1(s_1+s_2)-A_2{s}_3-A_3{s}_2$$

$$A_1{s}_1{s}_2$$

$$A_1=\frac{{\beta }_1{\beta }_2}{s_1{s}_2}=\frac{{\beta }_1{\beta }_2}{({\gamma }^2-{\xi }^2)/4}=\frac{{\beta }_1{\beta }_2}{{\beta }_1{\beta }_2-{\alpha }_1{\beta }_2-{\alpha }_2{\beta }_1}$$

$$-A_1{s}_3+A_2(s_2-s_3)={\beta }_1+{\beta }_2+s_2.$$

$$A_2$$

$$=\frac{4{\beta }_1{\beta }_2-2(\xi -\gamma )({\beta }_1+{\beta }_2)+(\xi -\gamma {)}^2}{2\xi (\xi -\gamma )}=\frac{(\xi -\gamma -2{\beta }_2)(\xi -\gamma -2{\beta }_1)}{2\xi (\xi -\gamma )}$$

$$=\frac{(\xi -{\alpha }_1-{\alpha }_2+{\beta }_1-{\beta }_2)(\xi -{\alpha }_1-{\alpha }_2-{\beta }_1+{\beta }_2)}{2\xi (\xi -\gamma )}.$$

$$A_3=\frac{(\xi +\gamma +2{\beta }_1)(\xi +\gamma +2{\beta }_2)}{2\xi (\xi +\gamma )}=\frac{(\xi +{\alpha }_1+{\alpha }_2+{\beta }_1-{\beta }_2)(\xi +{\alpha }_1+{\alpha }_2-{\beta }_1+{\beta }_2)}{2\xi (\xi +\gamma )}.$$

$$\phi (t)=\mu \left(A_1+A_2{\mathrm{e}}^{s_{2}t}+A_3{\mathrm{e}}^{s_{3}t}\right).$$

$$\mathbb{E}[N(t)]={\int }_0^t\phi (\tau )\mathrm{d}\tau$$

$$\log \mathcal{L}(t_1,\dots ,t_n|\theta )=-{\int }_0^T\lambda (t|\theta )\mathrm{d}t+{\int }_0^T\log (\lambda (t|\theta ))\mathrm{d}N(t).$$

$$\log \mathcal{L}(t_1,\dots ,t_n|\theta )=-\mu T$$

$$+\sum _{t_k<T}\log \left(\mu +\sum _{m=1}^P{\alpha }_m\sum _{t_i<t_k}{\mathrm{e}}^{-{\beta }_m(t_k-t_i)}\right).$$

$$\log \mathcal{L}(t_1,\dots ,t_n|\theta )$$

$$+\sum _{t_k\le t_n}\log \left(\mu +\sum _{m=1}^P{\alpha }_m{A}_m(k)\right),$$

Full text

Hawkes process, self-exciting process, model selection, information criterion, AIC, BIC, HQ

I Introduction

Technological advancement made it possible to record detailed data of all trades on financial markets. This development called for suitable econometric models that incorporate the time structure of durations between trades. Previously, models were designed such that this information was lost due to aggregation of data to equidistant time grids. However, empirical studies of high-frequency trading data show that intra-day trades have a typical pattern: there is high trading activity at the beginning and end of the trading day whereas there is low trading activity during lunch hours in the middle of the trading day (see for example Bertram_2004 ). Engle and Russell were among the first to propose a point process approach to modelling durations between trades (Engle_Russell_1997 ; Engle_Russell_1998 , Engle_2000 ). The proposed autoregressive conditional duration (ACD) model is closely related to the popular GARCH model for volatility clustering. It is also known under the name multiplicative error model (see Hautsch_2012 for further details on this topic).

However, self-exciting point processes have gained vast popularity among econometricians and financial mathematicians. Especially Hawkes processes (Hawkes_1971b, ; Hawkes_1971a, ) offered an intuitive notion of endogenous and exogenous components contributing to (trade) event clustering, which is sometimes referred to as “market reflexivity” (Filimonov_Sornette_2012 , Hardiman_2013 ). Additionally, from a theoretical point of view, Daley_2003 draw the analogy between the role of Hawkes processes spectral approximations of point processes and the importance of autoregressive models for mean square continuous processes.

Hawkes processes were originally used for seismic data (Hawkes_1973 , Ogata_1988 ), but their characteristic property of self-excitation and event clustering are appealing properties for mimicking similar phenomena found as stylized facts in intra-day financial data. Bowsher_2007 was among the early works to establish the connection between Hawkes processes and financial modelling. As there is an intensity-based as well as cluster-based definition, there exist various simulation and estimation techniques which take advantage of either perspective on Hawkes processes. To mention a few, for simulation we have the thinning approach (Ogata_1981, ), the time-change approach based on the random time-change theorem (Meyer_1971, ) and applied specifically to Hawkes processes for instance in (Ozaki_1979, ), exact simulation (Dassios_Zhao_2013, ) and perfect simulation (Moller_Rasmussen_2005, ). Concerning estimation techniques, the standard maximum likelihood approach can for example be found in Ozaki_1979 . Beyond that, Hawkes_1973 used a spectral estimation approach, Rasmussen_2013 proposes a Bayesian estimation technique and an application of the expectation maximization (EM) algorithm can be found in Veen_Schoenberg_2008 .

These tools for handling Hawkes processes numerically paved the way for applications on various types of financial data such as mid-price changes, order books, extreme price movements (among others) gathered from liquid stocks, futures, indices or foreign exchange markets. For details, the review paper by Bacry_2015 gives a very good summary of recent literature on Hawkes models in finance.

Essentially, for parametric estimation, there are two kernels which are widely used in the literature to fit financial data: the exponential kernel and the power law kernel. Whereas the power law asymptotics are additionally backed up by results from non-parametric estimation literature as in Bacry_2012 , the exponential kernel case is analytically more tractable and is still applied in recent literature (Hardiman_2013 , Rambaldi_2015 , Lallouache_Challet_2016 ).

Today’s computing power not only allows accurate recording of high-frequency trades, but enables us to fit almost arbitrarily complex models to previously gathered data. Recent proposals to model such data include intensities of Hawkes processes that can be expressed as weighted sums of exponential and power law kernels. The natural question arises as to how many terms should be included in such a model to be best suited in describing the data. Information criteria (IC) offer quantitative methods to discriminate between (possibly numerous) models. There are two competing objectives when it comes to selecting an “optimal” model order: On the one hand we would like to capture and describe the observed phenomena within the data as accurately as possible but, on the other hand, it is important to keep the complexity of the model to a minimum. A complex model can lead to numerical instabilities and superfluous parameters that do not carry much descriptive power. Information criteria are quantitative tools to manage this trade-off situation. Our aim in this paper is to test how well this model selection method could work for a Hawkes process intensity of weighted sums of exponential terms using simulated data.

The paper is organized as follows: Section II is devoted to the exponential Hawkes-P model. After a short definition and discussion of the average intensity, we move on to the simulation procedure and the parameter estimation method via maximum likelihood. In Section III, we give a short introduction to information criteria and discuss the consistency property. Finally, we describe the setup of the Monte-Carlo experiment and give the numerical results in Section IV.

II A Hawkes model with exponential kernels

For a self-exciting point process $(N(t))_{t\geq 0}$, the conditional intensity function is formally defined by

[TABLE]

where $\mathcal{F}_{t}$ represents the known history up to time $t$. We assume the conditional intensity function to be of the form (conditioning on history removed for the sake of simplicity)

[TABLE]

where we have for the response function $g(\tau)\geq 0$ $\forall\tau\in\mathbb{R}^{+}$ and $\mu>0$ is the baseline intensity. The term containing the response function can be identified with the self-excitation property and is therefore referred to as the endogenous part of the intensity whereas the baseline intensity is the exogenous part. The above intensity function defines a Hawkes process with finite past, as we assume the counting process $(N(t))_{t\geq 0}$ to start at [math]. Note that we deviate from the original definition, where usually the integral in (2) is evaluated over $(-\infty,t]$.

In particular, we are interested in the case when the response function can be written as a weighted sum of exponentials:

[TABLE]

Then, the intensity function is given by

[TABLE]

with $\mu,\alpha_{m},\beta_{m}>0$ and $\{t_{1},\ldots,t_{k}\}$ are the jump times of $N(t)$ up to time $t$. In short, we will call this process exponential Hawkes-$P$ process, where $P$ is the order of the process.

We will consider this class of Hawkes processes as a possible parametric model for durations between trades.

Average intensity: stationary vs. non-stationary case

In Hawkes_1971a the average intensity for a stationary Hawkes process with infinite past has been calculated to be

[TABLE]

where $n:=\int_{0}^{\infty}g(\nu)\,\mathrm{d}\nu$ is called the branching ratio. This result follows essentially by taking the expectation on both sides of (2). In particular, for the exponential kernel we have $n=\sum_{m=1}^{P}\frac{\alpha_{m}}{\beta_{m}}$ and the stationarity condition is $n<1$. The special case of $n=1$ also allows stationary processes which are treated in Bremaud_Massoulie_2001 .

However, for the Hawkes process with finite past associated with the intensity in (4), we apply a different approach using Laplace transforms: Let $\varphi(t):=\mathbb{E}[\lambda(t)]$ now be the average intensity function of a non-stationary Hawkes process. Then, taking expectations on both sides of (2) yields

[TABLE]

The Laplace transform of $\varphi$ is given by

[TABLE]

where in (7) we we are able to apply Fubini’s theorem since the integrand is positive. Finally, we have an algebraic equation which can be solved for $\tilde{\varphi}$:

[TABLE]

For $P>1$ we could write alternatively:

[TABLE]

This gives an analytic expression for the Laplace transform of the intensity function. From Equation (8) we can see that it is reasonable to demand the usual stationarity condition $\sum_{m=1}^{P}\frac{\alpha_{m}}{\beta_{m}}<1$ in order to ensure that the right hand side term is well defined.

In general, the evaluation of the average intensity function can be done by (numerical) Laplace inversion. However, for lower model orders (up to $P=4$) it is possible to invert the Laplace transform analytically. We will show this for first and second order in the following examples.

Example 1** (Formula for the average intensity in the case $P=1$).**

For $P=1$ the expression in (8) simplifies to

[TABLE]

where we used a partial fractions decomposition in the last step. This allows us to analytically invert the Laplace transform:

[TABLE]

Example 2** (Formula for the average intensity in the case $P=2$).**

For $P=2$ we have

[TABLE]

Starting from order $P=2$, the explicit formulas can be quite complicated.

Let $R$ and $Q$ denote the polynomial in the numerator and the denominator of the right hand side expression in (12) respectively. Then, assuming $Q$ has only real valued roots of single multiplicity denoted by $s_{1},s_{2},s_{3}$, the partial fractions decomposition is given by

[TABLE]

where

[TABLE]

The partial fractions decomposition implies that

[TABLE]

and comparing coefficients of $s^{2},s$ and $1$ on both sides of the equation yields

[TABLE]

Then we get

[TABLE]

by solving (19) for $A_{1}$ and inserting (14).

Now multiply (17) by $s_{2}$ and add (18) to get

[TABLE]

Solving for $A_{2}$ we get

[TABLE]

Multiplying (17) by $s_{3}$, adding (18) and following similar steps as for $A_{2}$ yield

[TABLE]

The Laplace inversion gives

[TABLE]

Note that with the condition $\sum_{m=1}^{P}\frac{\alpha_{m}}{\beta_{m}}<1$ it follows that the roots $s_{2}$ and $s_{3}$ are real and negative.

From both examples, we can see that for large times $t$ the exponential terms in Equations (11) and (24) become negligible and the remaining expressions coincide with the intensity function of the stationary case. In a small Monte-Carlo (MC) experiment, we simulated 1000 paths of a Hawkes process with 1000 events (see also empirAgg2.m). The parameters are $\mu=0.5$, $\alpha_{1}=3.1$, $\alpha_{2}=5.9$, $\beta_{1}=9.9$ and $\beta_{2}=10$. Figure 1 shows a plot of the empirically observed average number of events against the theoretically expected number of events. Plotting such figures might be useful for validation of a simulations algorithm. Recall the relation between average intensity function $\varphi$ and expected number of events of a point process $(N(t))_{t\geq 0}$:

[TABLE]

For small times we can observe the transient exponential behavior which vanishes for large times. In particular, the slope of the two theoretical functions are approximately equal for large times and indicate that the intensity function of the non-stationary case converges to the stationary case. Also, we can verify the edge effect when simulating a Hawkes process with finite past, which will be briefly discussed in the next section.

II.1 Simulation

As seen in the previous section, simulating a Hawkes process with finite past in order to approximate a Hawkes process with infinite past will cause the simulated process to be non-stationary at the beginning of the simulation time. This phenomenon is also known as edge effect as offspring of events that might have occurred in the past are omitted. For further details on this see Moller_Rasmussen_2005 ; Moller_Rasmussen_2006 .

However, similar to Dassios_Zhao_2013 , we explicitly want to work with a Hawkes process with finite past. Therefore, we view the edge effect as an inherent property of the model rather than an artifact of the simulation. Besides, the exact simulation algorithm in Dassios_Zhao_2013 , though applicable to multidimensional exponential models, does not directly apply to our proposed model due to the lack of identification of the exogenous and endogenous part of the intensity. This leaves us with the popular thinning algorithm going back to Lewis_Shedler_1979 and Ogata_1981 . We used an implementation of the thinning algorithm to simulate the process on a time interval $[0,T]$ (see hawkesThinning.m) and compare models up to order $3$. We first generate sample data that serve as a technical example for the estimation and model selection methods. The parameter settings are given in Table 1.

Connection to empirical findings in financial literature

In order to enhance the practical relevance of our experiments and results we would also like to use parameter settings which allow intensities which can also be observed in empirical studies. Concerning the exponential Hawkes-P model, Hardiman_2013 found that the use of the single exponential intensity function might give misleading results, which is also confirmed by Rambaldi_2015 . However, this does not necessarily hold for exponential Hawkes processes of higher order: Lallouache_Challet_2016 found that Hawkes models with exponential intensity kernels of order $P$ greater than one perform better than the single exponential model and comparably well to power law models when applied to FX data. This is why we include a parameter set that was estimated in this paper for our MC experiment (see Table 2).

II.2 Maximum likelihood estimation (MLE) and goodness of fit

Although we are primarily interested in the performance of model selection, we must make sure that the MLE gives reasonable results. This is because we expect a close connection between the quality of the MLE and the subsequent model selection result. A poor MLE due to numerical problems or lack of data is likely to compromise the model selection. For example, a correctly selected model order can be meaningless if the estimated model itself fails to describe and predict key features or quantities of the data we are interested in. In the following subsections we briefly present the fitting procedure as well as the root mean squared error as our chosen measure for goodness of fit.

II.2.1 Fitting via MLE

The fitting algorithm follows the theory in Ozaki_1979 which is a standard maximum likelihood procedure. For a self-exciting point process with intensity $\lambda$ the log-likelihood for data $0<t_{1}<\ldots<t_{n}<T$ is given by

[TABLE]

Let $\theta=(\mu,\alpha_{1},\ldots,\alpha_{P},\beta_{1}.\ldots,\beta_{P})$ be the vector of parameters for the Hawkes-$P$ model. Inserting Eq (4) into Eq (26) gives

[TABLE]

Moreover, Ozaki_1979 shows that the log-likelihood can be calculated recursively, which reduces the computational burden from $\mathcal{O}(n^{2})$ to $\mathcal{O}(n)$: Assume that $T=t_{n}$, i.e. the last event is the last time point of observation. Then

[TABLE]

To obtain the MLE of the parameters we maximize the log-likelihood function with respect to the parameters subject to the stationarity condition:

[TABLE]

We assume the $\beta$ parameters to be ordered to avoid identification problems. The maximization (or rather the minimization of the negative log-likelihood) is typically done numerically as the estimators are not available in closed form. We used the standard MATLABTM function fmincon for constrained problems. The optimization routine can be found in the supplementary files fitting.m, conditions.m and LogLik_iter.m.

II.2.2 Goodness of fit

Important asymptotic properties of the MLE for Hawkes processes have been studied and proven by Ogata_1978 (see also in the appendix in Rambaldi_2015 for a brief summary). In particular, we may assume the MLE to be consistent, i.e. with sample size tending to infinity the MLE converge to the true values of the parameters. In order to verify these results with our MC experiment, we use the RMSE (root mean squared error) as a measure for the goodness of fit: Let $\theta$ be a generic model parameter to be estimated and $\hat{\theta}$ the corresponding estimator. We are given $N=1000$ samples and have the parameter estimates $\hat{\theta}^{(k)}$, $k=1,\ldots,N$. For each sample we calculate the (absolute) root mean squared error to be

[TABLE]

and the relative root mean squared error

[TABLE]

It is easy to calculate the above quantities as the true model values are known in our mock data setting.

III Information criteria and model selection

In the following sections we will define the IC we are interested in and briefly describe relevant theoretical concepts. Based on that, we present the results of a simple MC experiment to assess model selection using IC. We are aware that certain conditions of our experiment are not given in reality and therefore also discuss the limitations to the conclusions we may draw from the numerical results.

Definitions and theoretical properties

This section follows introductory work which can be found in Claeskens_Hjort_2008 and references therein.

Definition 1**.**

For a given model fitted to data via MLE let $\mathcal{L}$ be the maximal log-likelihood value, $k$ the number of parameters and $n$ be the sample size of the data set. Then we define:

Akaike’s information criterion (AIC) (see Akaike_1973 )

[TABLE] 2. 2.

Bayesian information criterion (BIC) (see Schwarz_1978 )

[TABLE] 3. 3.

Hannan and Quinn information criterion (HQ) (see Hannan_1979 and Hannan_1980 )

[TABLE]

The exponential Hawkes-P process from Section II is an example of a nested series of models with $1+2P$ parameters and are therefore quite suitable for calculating information criteria. The formulas for the IC were implemented in the function IC.m. The above information criteria are of the form

[TABLE]

where $M_{k}$ is a model associated with parameter number $k$ and $c(k,n)$ is a suitably chosen penalty term that accounts for the complexity of the model, i.e. the number of parameters. Within a given set of models to choose from, the “best” model is the one which minimizes the IC value. In other words, the selected model should give the best fit to the data, i.e. have a large log-likelihood value, while being as parsimonious as possible, i.e. use few parameters. Therefore, formula (35) represents the trade-off situation we have discussed previously.

Remark 1*.*

The AIC was derived from estimating the Kullback-Leibler distance between the “true” model distribution and the estimated one. Hurvich_Tsai_1989 proposed a correction of the AIC for small samples:

[TABLE]

We shall follow the recommendation in Burnham_Anderson_2004 and use the AICc whenever $n<40k_{\text{max}}$ as a rule of thumb, where $k_{\text{max}}$ is the maximal number of parameters used among the candidate models. 2. 2.

The BIC was first derived in a Bayesian estimation approach, but is also valid in the frequentist context and there is an alternative derivation of the BIC from the frequentist perspective (see Burnham_Anderson_2004 for details) 3. 3.

The HQ is designed to have the slowest growing penalty term that still renders the IC to be strongly consistent (see later for a more precise definition). The proof makes use of the law of iterated logarithm. Besides, the HQ was originally defined more generally as

[TABLE]

but $c$ was chosen to be $1$ in a subsequent example. Claeskens_Hjort_2008 point out that the choice of $c$ is not clear and renders the information criterion less relevant for practitioners.

Similar to the consistency property of the MLE, it is a desirable property to have the IC selecting the correct model order with high probability when the underlying sample size increases. To be more precise:

Definition 2**.**

Let $n$ be the underlying sample size, $\mathcal{J}$ be the set of models among all competing models that minimize the Kullback-Leibler distance to the true model and let $\mathcal{J}_{0}\subset\mathcal{J}$ be the subset of models with minimal (parameter) dimension. Then, an IC is said to be consistent if there is a $j_{0}\in\mathcal{J}_{0}$ such that

[TABLE]

i.e. the probability that the IC will choose a model with smallest dimension minimizing the Kullback-Leibler distance converges to 1.

An IC is strongly consistent if the assertion in (38) holds almost surely:

[TABLE]

Remark 2*.*

The above definition follows the notation in (Claeskens_Hjort_2008, , p. 101), but the original proof of sufficient conditions for consistency and strong consistency are shown in Sin_White_1996 , (where consistency actually goes under the name of weak consistency).

As a matter of fact, the AIC fails to be consistent as the penalty term does not depend on the sample size. The asymptotic distribution of the associated model selection was analyzed for autoregressive models for example in Shibata_1976 . BIC and HQ on the other hand are found to be strongly consistent. As a consequence, their asymptotic distribution of model selection is bound to converge to a delta on the most parsimonious Kullback-Leibler minimizing model. The respective convergence rates for AIC and BIC were analyzed in Zhang_1993 for another regression model.

Consistency from a practical perspective

From the previous section one might conclude that the non-consistent AIC would be inferior to the consistent BIC and HQ. However, the situation is more complicated: We have to keep in mind that consistency is an asymptotic property. This means that in theory the consistent IC will eventually outperform the AIC for almost all cases if the sample size is sufficiently large. Unfortunately, practitioners just have a limited amount of data available and it is very difficult to judge whether the sample size belongs to the asymptotic region. Indeed, empirical studies suggest for various statistical models that the AIC outperform the BIC in small sample cases111By “small samples” we refer to the situation that the sample size is not sufficiently large enough for the asymptotic consistency results to hold, but large enough such that effects similar to the paradox discussed by Freedman_1983 can be safely excluded.: As an example among regression models, Hurvich_Tsai_1989 ; Hurvich_Tsai_1990 compared different IC on simulated data especially to promote the (still inconsistent) AICc as a modification of the AIC for smaller samples. More recently, Javed_Mantalos_2013 applied IC (AIC, BIC, HQ, AICc) in a MC simulation of (nonlinear) GARCH models. Their results suggest that the AIC outperforms the BIC and HQ for higher-order GARCH processes.

As a consequence of the above discussion, we can make the idea and objective of our MC experiment more precise: First, we need to point out that the numerical results of the simplistic setting of our MC experiment do not directly translate to how empirical data should be handled. IC are one of many tools for model-selection and cross-validation. We do not expect to find a “best” IC, but rather want to verify the theoretical properties of the different IC for Hawkes processes. In particular, due to the fact that most theoretical results have been derived for regression models only, our work may help to shed light on asymptotic regions and convergence rates of consistent IC and the asymptotic distribution of selected orders of the AIC for this model class. The verification of theoretical properties will be the main aim for the MC simulation using Parameter Set 1 whereas for our empirical Parameter Set 2 we following the advice given in Burnham_Anderson_2004 and use realistic sample sizes.

IV Numerical results

Using the thinning algorithm described in Section II.1 we simulated four different data sets containing 1000 samples. Three of them correspond to each row of Parameter Set 1 in Table 1 and one data set consists of samples of an exponential Hawkes-2 model with parameter values from Parameter Set 2 shown in Table 2. Especially for Parameter Set 2 the time horizon $T$ can be assumed to be given in seconds. It ranges from $10\text{\,}\mathrm{min}$ to $6\text{\,}\mathrm{h}$ to reflect typical intra-day financial data sets.

In order to check how well the estimation method works for our parameter sets, we first assume that the correct model order $P$ is known and run a MLE of the parameters of the true model underlying each data set. Subsequently, we are able to calculate the RMSE as a measure of distance between the true and the estimated parameter values. The absolute and relative RMSE values for Parameter Set 1 can be found in Table 3 and Table 4 respectively. For Parameter Set 2, see Tables 8 and 9. We observe that the RMSE decreases with increasing sample size. This is to be expected as the MLE is known to be consistent.

Finally, we assume that the true model order is not known, but needs to be selected by the IC. Consequently, for each data set we have to fit all possible model orders $P=1,2,3$ and to calculate the associated IC values. In the following we discuss the results of the model selection.

We first consider Parameter Set 1. For simulated data with model order $P=1$ we can see in Table 4 that the relative RMSE is comparably low even for the smallest samples corresponding to the time horizon $T=500$. The model selection in Table 5 confirms that the smallest sample size might already be enough to guarantee high success rates (over $90\%$) of all IC. Nevertheless, already in the lowest order case, we can observe the different behavior of consistent and inconsistent IC. For BIC and HQ, the success rate improves with increasing average sample size. In particular, the relation seems to be monotone and, in the case of the BIC, the success rate reaches $100\%$ already for $T=1000$, The HQ performs slightly worse than BIC, but is still well over $90\%$ and very close to $100\%$ for $T=5000$. However, the AIC behaves in a more concerning manner. Even for large sample sizes the model selection using the AIC allows a comparably large probability ($>6\%$) to select a higher order than $P=1$. As the AIC is not a consistent IC, we cannot exclude the possibility that these results already approximate the asymptotic distribution of model selection of the AIC. As mentioned earlier this asymptotic distribution is typically different from the delta distribution with mass one on the true model order. Additionally, the numerical results show that increasing the average sample size does not necessarily increase the success rate of model selection. For instance, moving from $T=500$ to $T=1000$ we can observe a decrease in success rate in the AIC case.

In the case of model order $P=2$, there is the possibility of both over- and underestimation. We observe quite large RMSE for the parameters $\alpha_{1}$ and $\beta_{1}$, especially for smaller samples corresponding to $T=500$ and $T=1000$. This could be one of the factors affecting the model selection for $T=500$ in Table 6: there is a significant proportion of underestimation among all IC, most notably the high underestimation rate of almost $95\%$ of the BIC. The AIC seems to perform best in this setting for $T=500$ with success rates slightly above $50\%$, but also with $48\%$ underestimation. For larger samples, the BIC and HQ select the correct model order with very high probability (around $90\%$ or even larger) and the BIC reaches $100\%$ success rate at $T=2000$. Again, we have the adverse effect that the success rates of the AIC decrease with growing average sample size. Even for the largest average sample size for $T=5000$ there is a relatively large probability of overestimation of over $6\%$.

For data simulated with $P=3$ we have a similar behavior as with $P=2$. Again, Table 4 reports large RMSE for the parameters $\alpha_{1}$ and $\beta_{1}$ in small sample cases. As $P=3$ is the highest selectable model order, this excludes cases of overspecification. This means that the we can observe the same pattern in model selection of the AIC as for the BIC and HQ: Starting at $T=500$, there are mostly cases of underestimates followed by improving success rates as the sample size increases. All IC reach $100\%$ success rate for $T=2000$. However, it is very likely that we would be able to observe the tendency of the AIC to overestimate if we included higher orders $P>3$ in the model selection set.

When working with Paramerter Set 2, we chose the time horizons $10\text{\,}\mathrm{min}$, $15\text{\,}\mathrm{min}$, $30\text{\,}\mathrm{min}$, $1\text{\,}\mathrm{h}$, $3\text{\,}\mathrm{h}$ and $6\text{\,}\mathrm{h}$. At first, there are large RMSE values for $T=600$ and $T=900$ (see Tables 8 and 9), which shows that the sample sizes are so small that we cannot ensure good estimates of the MLE method. Especially estimates of $\alpha_{2}$ and $\beta_{2}$ have large RMSE. This situation corresponds to the case $T=500$ in the setting of Parameter Set 1. When we compare with the corresponding model selection in Table 10, we observe the same phenomenon of underestimation is most severe for the BIC, less for the HQ and least for the AIC. As samples are quite small for these cases and may fulfill the rule of thumb discussed in point (i) in Remark 1, we included the combined model selection rule AICc/AIC in the table. It applies the AICc whenever $n<3\cdot 40=120$ and the AIC otherwise. The numerical results for the combined AICc/AIC selection rule are very similar to the standalone AIC and even slightly worse for T=600 and T=900.

When we move on to larger samples from $30\text{\,}\mathrm{min}$ to $1\text{\,}\mathrm{h}$, there is a noticeable change in the RMSE values. More precisely, the RMSE values decrease faster for the second exponential term, i.e. $\alpha_{2}$ and $\beta_{2}$, which leads to the first exponential term with $\alpha_{1}$ and $\beta_{1}$ to contribute more to the overall estimation error. There is a noticeable increase in the rate of correct model selection among all IC ranging over $90\%$ for $T=3600$.

Finally, for large samples with time horizons from $3\text{\,}\mathrm{h}$ up to $6\text{\,}\mathrm{h}$ represent data of half up to an entire trading day respectively. The relative RMSE of each parameter is less than $20\%$ and the rate of correct model selection for the consistent IC (BIC and HQ) is close to $100\%$. However, the success rate of the AIC decreases to about $94\%$ with a $6\%$ probability of overestimation.

V Summary and Outlook

Concerning the performance of model selection, the results of our MC experiment can be summarized as follows. In alignment with similar studies for regression models, we can observe that the inconsistent AIC outperforms the other two IC when the MLE is applied to smaller samples. In contrast, the consistent IC (BIC and HQ) perform excellently for sufficiently large samples and we can observe a monotonic improvement in their success rate when increasing the sample size. In spite of the concerns presented above in Remark 1 (iii), the numerical results show that the HQ should not be excluded as a well performing IC. More concerning, and not as commonly observed in previous studies, is a non-monotonic relation between sample size and success rate of model selection of the AIC.

Future research in this direction can be on the equally popular power law intensities and further model selection methods like the focused information criterion as well as model averaging.

Funding

E. Scalas and M. Trinh were supported by the Strategic Development Fund of the University of Sussex.

Supplemental material

Please contact either Enrico Scalas or Mailan Trinh.

[TABLE]

V.1 Appendices