# A Fractional Hawkes process

We modify ETAS models by replacing the Pareto-like kernel proposed by Ogata with a Mittag-Leffler type kernel. Provided that the kernel decays as a power law with exponent β + 1 ∈ (1,2], this replacement has the advantage that the Laplace transform of the Mittag-Leffler function is known explicitly, leading to simpler calculation of relevant quantities.

• 6 publications
• 1 publication
• 1 publication
09/22/2020

### Deep Neural Tangent Kernel and Laplace Kernel Have the Same RKHS

We prove that the reproducing kernel Hilbert spaces (RKHS) of a deep neu...
04/30/2019

### First digit law from Laplace transform

The occurrence of digits 1 through 9 as the leftmost nonzero digit of nu...
10/04/2021

### Generalized Kernel Thinning

The kernel thinning (KT) algorithm of Dwivedi and Mackey (2021) compress...
06/08/2016

### Rotation Invariant Angular Descriptor Via A Bandlimited Gaussian-like Kernel

We present a new smooth, Gaussian-like kernel that allows the kernel den...
09/08/2020

### Controlled Drift Estimation in the Mixed Fractional Ornestein-Uhlenbeck Process

This paper is devoted to the determination of the asymptotical optimal i...
10/29/2010

### Fractionally Predictive Spiking Neurons

Recent experimental work has suggested that the neural firing rate can b...
02/03/2020

### Stochastic geometry to generalize the Mondrian Process

The Mondrian process is a stochastic process that produces a recursive p...

## I Introduction

In 1971, Hawkes (Hawkes (1971a, b)) introduced a class of self-exciting processes to model contagious processes. In their simpler version, these are point processes with the following conditional intensity

 λ(t|Ht)=limh→0E(N(t,t+h)|Ht)h=λ+α∫t−∞f(t−u)dN(u),

where , is a Hawkes self-exciting 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 power-law kernel for self-exciting 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 power-law kernel with a Mittag-Leffler kernel and this will be the main contribution of this chapter. We will first introduce the Mittag-Leffler 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 one-parameter Mittag-Leffler 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 Mittag-Leffler distributed positive random variables

Consider the one-parameter Mittag-Leffler function

 Eβ(z):=∞∑n=0znΓ(nβ+1), (II.1)

with . If computed on for , the Mittag-Leffler 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 power-law with index

for large times. Its sign-changed first derivative

 fβ(t):=−dEβ(−tβ)dt=tβ−1Eβ,β(−tβ) (II.2)

is the probability density function of , where is the two-parameter Mittag-Leffler function defined as

 Eγ,δ(z):=∞∑n=0znΓ(nγ+δ). (II.3)

Notice that the one-parameter Mittag-Leffler function coincides with the two parameter one with and . For a suitable function defined for positive , let us define its Laplace transform as

 ~f(s)=L(f(t),s)=∫∞0f(t)e−stdt.

The functions and have explicit Laplace transforms. The survival function has the following Laplace transform

 L(Eβ(−tβ);s)=sβ−11+sβ, (II.4)

and the probability density function has the following Laplace transform

 L(fβ(t);s)=11+sβ. (II.5)

Moreover, they have an explicit representation as an infinite (actually continuous) sum of exponential functions Gorenflo (1997);

 Eβ(−tβ)=∫∞0e−θtKβ(θ)dθ, (II.6)

with

 Kβ(θ)=1πθβ−1sin(βπ)θ2β+2θβcos(βπ)+1, (II.7)

 fβ(t)=∫∞0θe−θtKβ(θ)dθ. (II.8)

These functions play an important role in fractional calculus. For instance is the solution of the following anomalous relaxation problem

 dβg(t)dtβ=−g(t), (II.9)

where is the Caputo derivative defined as

 dβg(t)dtβ=1Γ(1−β)ddt∫t0g(τ)(t−τ)βdt−t−βΓ(1−β)g(0+), (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

 λ(t|Ht)=limh→0E(N(t,t+h)|Ht)h=λ+α∫t−∞fβ(t−u)dN(u), (III.1)

where and is a Hawkes self-exciting point process, leading to

 λ(t|Ht)=λ+α∑ti

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 time-changed intensity process where, in his case, the conditional intensity of the self-exciting process is the solution of a mean-reverting stochastic differential equation

 dλt=κ(θ−λt)dt+ηdPt,

where , and are suitable parameters and the driving process is given by

 Pt:=Nt∑k=1ξk,

where is a counting process and

s are independent and identically distributed marks with finite positive mean and finite variance. The time-change

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 power-law 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 .

 μ(τ)=α[Λfβ(τ)+∫τ0fβ(τ−v)μ(v)dv+∫∞0fβ(τ+v)μ(v)dv], (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

 ~μ(s)=α[Λ~fβ(s)+~fβ(s)~μ(s)+L(∫∞0fβ(τ+v)μ(v)dv;s)]. (III.4)

Now, using (II.8) and setting

 h(θ)=θKβ(θ)=1πθβsin(βπ)θ2β+2θβcos(βπ)+1,

we can write

 fβ(τ)=∫∞0h(θ)e−θτdθ,

so that the last summand in (III.4) becomes

 L(∫∞0fβ(τ+v)μ(v)dv;s) = ∫∞0e−sτ[∫∞0[∫∞0h(θ)e−θ(τ+v)dθ]μ(v)dv]dτ (III.5) = ∫∞0h(θ)1θ+s~μ(θ)dθ.

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

 f(ω)=12π∫e−iωτμ(c)(τ)dτ, (III.6)

where, because if events cannot occur multiply, the complete covariance density contains a delta function

 μ(c)(τ)=Λδ(t)+μ(t).

Then it is shown in Hawkes (1971b) p. 441 that

 f(ω)=Λ2π(1−G(ω))(1−G(−ω)), (III.7)

where, in our case, we would have

 G(ω)=∫∞0e−iωταfβ(τ)dτ=α1+(iω)β.

The proof of this result depended on the assumption that the exciting kernel decays exponentially asymptotically. However, this is not true for the Mittag-Leffler 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

 ~μ(c)(s)=Λ(1−~Φ(s))(1−~Φ(−s)), (III.8)

where

 ~Φ(s)=∫∞0e−sτΦ(τ)dτ

is the Laplace transform of the exciting kernel. In our case for , we have

 Φ(τ)=αfβ(τ)

and

 ~Φ(s)=α1+sβ.

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 well-behaved functions.

### iii.2 Intensity expectation

Let us consider the expectation for both a stationary and non-stationary fractional Hawkes process. In the stationary case (process from ), from equation (III.1), we get , leading to

 Λ=λ1−α. (III.9)

On the contrary in the non-stationary case (process from ), we can modify equation (III.1) as follows

 λ(t|Ht)=λ+α∫t0fβ(t−u)dN(u), (III.10)

so that the time-dependent expectation obeys the equation

 Λ(t)=λ+α∫t0fβ(t−u)Λ(u)du. (III.11)

Taking Laplace transforms, we get

 ~Λ(s)=λs+α~fβ(s)~Λ(s),

so that

 ~Λ(s)=λs11−α~fβ(s).

Using equation (II.5) yields

 ~Λ(s)=λs1+sβ(1−α)+sβ. (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

 E[N(t)]=∫t0Λ(τ)dτ.

Also, based on a continuous version of Hardy-Littlewood Tauberian theorem Feller (1971) we get

 limt→∞Λ(t)=λ1−α

as given by equation (III.9). This result is exemplified in Fig. 1 for , and . In that case, we get, for

 Λ(t)=2−et/4erfc(√t/2).

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 Mittag-Leffler 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 Mittag-Leffler functions described above and can be retrieved from the Matlab file exchange. The algorithm is as follows

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

2. Compute , for some small value of .

3. Generate a positive exponentially distributed random variable

with the meaning of a waiting time, with rate .

4. Set .

5. Generate a uniform random variate between and .

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

8. 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=1e-10;
beta=0.7;

SimPoints=[];

while t<T
t
M=mu+sum(alpha*ml(-(t+epsilon-SimPoints).beta,beta,beta,1). …
*(t+epsilon-SimPoints).(beta-1));
E=exprnd(1/M,1,1);
t=t+E;
U=rand;
if (U<(((mu+sum(alpha*ml(-(t-SimPoints).beta,beta,beta,1). …
*(t-SimPoints).(beta-1))))/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 power-law kernels of Pareto type as originally suggested by Ogata. We also presented an explicit simulation of the intensity process based on the so-called 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 non-parametric 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 Mittag-Leffler 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 self-exciting 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 continuous-time finance II: the waiting-time 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.