I Introduction
In 1971, Hawkes (Hawkes (1971a, b)) introduced a class of selfexciting processes to model contagious processes. In their simpler version, these are point processes with the following conditional intensity
where , is a Hawkes selfexciting counting process, represents the history of the process, is a branching ratio that must be smaller than 1 for stability, and is a suitable kernel (
must be a probability density function for a positive random variable).
In 1988, Ogata (Ogata (1988)) proposed the use of a powerlaw kernel for selfexciting processes of Hawkes type in order to reproduce the empirical Omori law for earthquakes. Ogata’s models are also known as Epidemic Type Aftershock Sequence models or ETAS models. Within this framework, it is natural to replace Ogata’s powerlaw kernel with a MittagLeffler kernel and this will be the main contribution of this chapter. We will first introduce the MittagLeffler distribution for positive random variables, then we will define a “fractional” version of Hawkes processes. Spectral properties and intensity expectation will be discussed using the fact that the Laplace transform of the oneparameter MittagLeffler function of argument with is known analytically. Finally, we will present a simple algorithm based on the thinning method by Ogata Ogata (1981) that simulates the conditional intensity process.
Ii MittagLeffler distributed positive random variables
Consider the oneparameter MittagLeffler function
(II.1) 
with . If computed on for , the MittagLeffler function has the meaning of survival function for a positive random variable
with infinite mean. This function interpolates between a stretched exponential for small times and a powerlaw with index
for large times. Its signchanged first derivative(II.2) 
is the probability density function of , where is the twoparameter MittagLeffler function defined as
(II.3) 
Notice that the oneparameter MittagLeffler function coincides with the two parameter one with and . For a suitable function defined for positive , let us define its Laplace transform as
The functions and have explicit Laplace transforms. The survival function has the following Laplace transform
(II.4) 
and the probability density function has the following Laplace transform
(II.5) 
Moreover, they have an explicit representation as an infinite (actually continuous) sum of exponential functions Gorenflo (1997);
(II.6) 
with
(II.7) 
leading to
(II.8) 
These functions play an important role in fractional calculus. For instance is the solution of the following anomalous relaxation problem
(II.9) 
where is the Caputo derivative defined as
(II.10) 
with initial condition .
Iii The fractional Hawkes processes
It becomes natural to use as kernel for a version of Hawkes processes that we can call fractional Hawkes processes. The conditional intensity of the process is given by
(III.1) 
where and is a Hawkes selfexciting point process, leading to
(III.2) 
with the branching ratio for stability. Hainaut Hainaut (2019) gives a different definition of fractional Hawkes process. Let us use his notation in the following remark. In his paper, he considers the timechanged intensity process where, in his case, the conditional intensity of the selfexciting process is the solution of a meanreverting stochastic differential equation
where , and are suitable parameters and the driving process is given by
where is a counting process and
s are independent and identically distributed marks with finite positive mean and finite variance. The timechange
is the inverse of a stable subordinator. Our definition (III.1) is much simpler and it is directly connected with ETAS processes given that the kernel has powerlaw tail with index . Mainardi (2000). In particular, given the explicit Laplace transform of and its representation in terms of infinite sum of exponentials, we can derive some explicit formulas.iii.1 Spectral properties
From equation (11) in Hawkes (1971a), we get the following equation for the covariance density for .
(III.3) 
where represents the asymptotic stationary value of the conditional intensity as derived in equation (3.9) below. If we now take the Laplace transform of (III.3), we get
(III.4) 
Now, using (II.8) and setting
we can write
so that the last summand in (III.4) becomes
(III.5)  
In principle, a numerical approximation of the integral in equation (III.5) plugged into (III.4) and coupled with equation (II.5) can lead to an explicit approximate expression for the Laplace transform . This will be the subject of a further paper. An alternative approach is given in Hawkes (1971b) where the Bartlett spectrum is defined for real as
(III.6) 
where, because if events cannot occur multiply, the complete covariance density contains a delta function
Then it is shown in Hawkes (1971b) p. 441 that
(III.7) 
where, in our case, we would have
The proof of this result depended on the assumption that the exciting kernel decays exponentially asymptotically. However, this is not true for the MittagLeffler distribution, which decays as a power law. Bacry and Muzy Bacry (2016) prove a more general result using Laplace transforms in the complex plane, more easily digested from section 2.3.1 in Bacry (2015). Then the Laplace transform of the covariance density is given by
(III.8) 
where
is the Laplace transform of the exciting kernel. In our case for , we have
and
Equations (III.7) and (III.8) look much the same, apart from a change of notation and a factor . The difference, however, is that in (III.7) we are dealing with real while, in (III.8), is a general complex variable and we can choose its domain to obtain wellbehaved functions.
iii.2 Intensity expectation
Let us consider the expectation for both a stationary and nonstationary fractional Hawkes process. In the stationary case (process from ), from equation (III.1), we get , leading to
(III.9) 
On the contrary in the nonstationary case (process from ), we can modify equation (III.1) as follows
(III.10) 
so that the timedependent expectation obeys the equation
(III.11) 
Taking Laplace transforms, we get
so that
Using equation (II.5) yields
(III.12) 
Equation (III.12) can be inverted numerically (or analytically for ) to give as well as the expected number of events from to as
Also, based on a continuous version of HardyLittlewood Tauberian theorem Feller (1971) we get
as given by equation (III.9). This result is exemplified in Fig. 1 for , and . In that case, we get, for
From Fig. 1, one can see that goes up very fast at first and, then, slowly converges to its asymptotic value . This is presumably due to the long tail of the MittagLeffler kernel.
Iv Simulation
In order to simulate the intensity process introduced in equation (III.1), we use the thinning algorithm introduced by Ogata Ogata (1981) (see also Zhuang and Touati 2015 report Zhuang (2015)). The function ml.m described in Garrappa (2015) is needed to compute the MittagLeffler functions described above and can be retrieved from the Matlab file exchange. The algorithm is as follows

Set the initial time , a counter and a final time .

Compute , for some small value of .

Generate a positive exponentially distributed random variable
with the meaning of a waiting time, with rate . 
Set .

Generate a uniform random variate between and .

If , set and update time to , else, just set .

Return to step 2 until exceeds .

Return the set of event times (or epochs)
An implementation of this algorithm in Matlab is presented below.
T=5;
t=0;
n=0;
alpha=0.9; % For stability
mu=1;
epsilon=1e10;
beta=0.7;
SimPoints=[];
while t<T
t
M=mu+sum(alpha*ml((t+epsilonSimPoints).beta,beta,beta,1). …
*(t+epsilonSimPoints).(beta1));
E=exprnd(1/M,1,1);
t=t+E;
U=rand;
if (U<(((mu+sum(alpha*ml((tSimPoints).beta,beta,beta,1). …
*(tSimPoints).(beta1))))/M))
n=n+1;
SimPoints = [SimPoints, t];
end
end
index=find(SimPoints<10);
SimPoints=SimPoints(index);
Two examples of intensity process simulated up to are presented in Fig. 2 and Fig. 3 for and , respectively.
V Outlook
In this paper, we defined a “fractional” Hawkes process and we studied its spectral properties and the expectation of its intensity. Thanks to explicit expressions for Laplace transforms, we could derive some analytical expressions that are only asymptotically available for powerlaw kernels of Pareto type as originally suggested by Ogata. We also presented an explicit simulation of the intensity process based on the socalled thinning method.
Further work is needed to better characterize our process. In particular, we did not deal with parameter estimation, the multivariate version of the process, and we did not use the model to fit earthquake data (or any other data for what matters including financial data). We do hope that all this and more can become the subject of an extensive future paper on this process.
References
 Bacry (2015) E. Bacry, I. Mastromatteo, and J.F. Muzy, Hawkes processes in finance, Market Microstructure and Liquidity, 1, 2015. DOI: http://dx.doi.org/10.1142/S2382626615500057, pp59.
 Bacry (2016) E. Bacry, and J.F. Muzy, First and second order statistics characterization of Hawkes processes and nonparametric estimation, IEEE Transactions On Information Theory, 62(4), 2184–2201, 2016.

Feller (1971)
W. Feller, An introduction to probability theory and its applications. Vol. II. Second edition, John Wiley & Sons, 1971.
 Garrappa (2015) R. Garrappa, Numerical evaluation of two and three parameter MittagLeffler functions, SIAM Journal of Numerical Analysis, 53(3), 1350–1369, 2015.
 Gorenflo (1997) R. Gorenflo and F. Mainardi, Fractional Calculus: Integral and Differential Equations of Fractional Order in Fractals and Fractional Calculus in Continuum Mechanics, A. Carpinteri and F. Mainardi (eds.) Springer 223–276, 1997.
 Hainaut (2019) D. Hainaut, Fractional Hawkes processes, UCLouvain preprint, 2019.
 Hawkes (1971a) A.G. Hawkes, Spectra of some selfexciting and mutually exciting point processes, Biometrika 58, 83–90, 1971.
 Hawkes (1971b) A.G. Hawkes, Spectra of some mutually exciting point processes, J. Royal Statistical Society B 33(3), 438–443, 1971.
 Mainardi (2000) F. Mainardi, M. Raberto, R. Gorenflo and E. Scalas, Fractional calculus and continuoustime finance II: the waitingtime distribution, Physica A, 287, 468–481, 2000.
 Ogata (1981) Y. Ogata, On Lewis’ simulation method for point processes, IEEE Transactions on Information Theory, 27, 23–31, 1981.
 Ogata (1988) Y. Ogata, Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes, 83(401), 9–27, 1988.
 Zhuang (2015) J. Zhuang and S. Touati, Stochastic simulation of earthquake catalogs, Community Online Resource for Statistical Seismicity Analysis, 2015. Available at http://www.corssa.org.