Exact Simulation of Hawkes Process with Exponentially Decaying Intensity

We introduce a numerically efficient simulation algorithm for Hawkes process with exponentially decaying intensity, a special case of general Hawkes process that is most widely implemented in practice. This computational method is able to exactly generate the point process and intensity process, by sampling interarrival-times directly via the underlying analytic distribution functions without numerical inverse, and hence avoids simulating intensity paths and introducing discretisation bias. Moreover, it is flexible to generate points with either stationary or non-stationary intensity, starting from any arbitrary time with any arbitrary initial intensity. It is also straightforward to implement, and can easily extend to multi-dimensional versions, for further applications in modelling contagion risk or clustering arrival of events in finance, insurance, economics and many other fields. Simulation algorithms for one dimension and multi-dimension are represented, with numerical examples of univariate and bivariate processes provided as illustrations.


Introduction
Throughout the recent subprime mortgage crisis and current European debt crisis, financial deterioration and losses can easily spread through the business network and financial market, and the risk of contagion has created enormous instability and uncertainty substantially threatening the financial systems, both regionally and globally. Financial models capable to capture the contagion or clustering effect are essentially needed. Selfexcited point processes, in particular, Hawkes processes early introduced by Hawkes (1971), nowadays become very popular and viable mathematical tools for modelling contagion risk and clustering arrival of events in finance, insurance and economics, see examples of applications in Errais et al. (2010), Embrechts et al. (2011), Dassios and Zhao (2012), and more recently, Aït-Sahalia and Hurd (2012) and Aït-Sahalia et al. (2013). They also have various applications in many other fields such as seismology and neurophysiology, see Ogata (1988) and Chornoboy et al. (1988).
For statistical analysis or practical implementations, simulation is a crucial issue as the procedure of the parameter estimation or calibration, in particular, for higher dimensional processes, heavily relies on it. Risk management and asset pricing for complex path-dependent financial derivatives, such as collateralized debt obligations (CDOs), require a high level of computational efficiency of simulation, as the analytic formulas are rather limited and often difficult to obtain. For instance, the model parameters for pricing CDOs can be calibrated to the market data via numerically solving optimisation problem using simulation, see Giesecke and Kim (2007).
There are two major simulation approaches in the literature: so-called intensity-based and cluster-based, since a Hawkes process can either be defined via a conditional stochastic intensity representation (e.g. Definition 2.1), or, via a Poisson cluster representation (e.g. Definition 2.2). As recently surveyed by Liniger (2009), the current standard algorithms for simulating Hawkes processes mostly remain the conventional thinning method via the acceptance/rejection procedure early introduced by Lewis and Shedler (1979) 1 that is generally used for simulating any point process with stochastic intensity. Møller and Rasmussen (2005) provided algorithms named perfect simulation 2 for simulating (marked) Hawkes processes (on R-time) using a cluster-based approach with a branching structure. Their algorithms start from time minus infinity and require certain stationarity condition, also see Møller and Rasmussen (2006). Giesecke and Kim (2007) introduced an algorithm named exact simulation 3 for Hawkes process (on R + -time) which is capable to sample interarrival-times directly via their underlying analytic distribution functions and hence avoid generating intensity paths and introducing discretisation bias for associated estimators, also see Giesecke et al. (2011). However, their method requires numerical evaluation of the inverse of these distribution functions via Brent's method which involves intensive computations.
In this paper, we focus on simulation algorithms rather than associated applications.
1 Also see Ogata's modified thinning algorithm by Ogata (1981) and Algorithm 7.5.IV in Daley and Vere-Jones (2002). 2 Here 'perfect' refers to the fact that the simulation on a finite time interval takes the past into account (without edge effects), also see Brix and Kendall (2002). 3 Here the 'exact' simulation means a method of drawing an unbiased associated estimator throughout the entire simulation process.
We have tailored a numerically efficient algorithm specifically for a Hawkes process (on R + -time) with exponentially decaying intensity, a special case of general Hawkes process that is most widely implemented in practice. Especially for this special case, our approach has many advantages, in particular, 1. it is able to exactly generate Hawkes processes by sampling interarrival-times directly via the underlying analytic distribution functions without numerical inverse, as each of these interarrival-times can be simulated exactly by splitting into two independent simple random variables; 2. it is straightforward to implement, and easily to generalise into higher dimensions; 3. it is flexible to simulate Hawkes processes starting from any time, conditional on any arbitrary fixed initial intensities, or, with any arbitrary distributions for the initial intensities; 4. it does not require stationarity condition for intensity dynamics, and hence Hawkes processes with unbounded intensities in a finite time can also be simulated; 5. it also can be easily adjusted to simulate Hawkes processes with stationary intensities.
The paper is organised as follows. Section 2 briefly reviews the Hawkes process with exponentially decaying intensity via definitions of intensity-based and cluster-based representations respectively, and some basic distributional properties are given. Section 3 provides the numerical algorithm of exact simulation for a Hawkes process in one dimension. For verification, numerical examples and comparison with exact theoretical results of means and variances are also given in Section 4. Section 5 represents how straightforward our method can be extended to a multi-dimensional version, and a numerical example of bivariate case is provided as an illustration. Hawkes (1971) introduced a general type of Hawkes process, where the current intensity of the arrival of points is determined by points in the past with associated weights. In this paper, we focus on a simple and slightly modified version, when the response function is decaying exponentially, and sizes of self-excited jumps in the intensity process are allowed to be either fixed or random. The randomness in jump sizes could provide higher flexibility for measuring the self-exciting effect, and this is a minor generalisation of the Markovian self-exciting process represented by Oakes (1975).

Preliminaries
Our aim is to develop algorithms for practical implementation and we assume the Hawkes processes start from time zero. The definition of this univariate Markovian Hawkes process on R + -time via conditional stochastic intensity representation is given in Definition 2.1.

Definition 2.1 (Intensity-based) Hawkes process with exponentially decaying intensity is
is a history of the process N t with respect to which {λ t } t≥0 is adapted; • a ≥ 0 is the constant reversion level; • λ 0 > 0 is the initial intensity at time t = 0; • δ > 0 is the constant rate of exponential decay; • {Y k } k=1,2,... are sizes of self-excited jumps, a sequence of i.i.d. positive random variables with distribution function G(y), y > 0; • {T k } k=1,2,... and {Y k } k=1,2,... are assumed to be independent of each other.
Jump-sizes {Y k } k=1,2,... in the intensity measure the levels of compact of contagion and can also be fixed constants; δ captures the persistence of contagion with a common rate of exponentially decaying after each jump. A sample figure of joint process (N t , λ t ) is represented in Figure 1.
On the other hand, this Hawkes point process can also be alternatively constructed as a marked Poisson cluster process on R + -time with the clusters following a recursive branching structure, see Hawkes and Oakes (1974), Daley and Vere-Jones (2002) and Rasmussen (2011). The definition via marked Poisson cluster representation is given by Definition 2.2, and offers a nice economic interpretation: the immigrants and offsprings could be considered as primary shocks and associated aftershocks respectively. Note that, as this point process is defined on R + , there are no immigrants before time 0 and hence no edge effects as concerned by Møller and Rasmussen (2005). with the following branching structure: • The immigrant and its mark (T m , Y m ) is said to be of generation 0.
• Recursively, given generations 0, 1, ..., n in C m , each (T j , Y j ) ∈ C m of generation n generates a Poisson process of offspring of generation n + 1 on (T j , ∞) with intensity Y j e −δ(t−T j ) , t > T j , where mark Y j ∼ G, independent of generations 0, 1, ..., n.
(e) C consists of the union of all clusters, i.e. C = S m=1,2,... C m .
Definition 2.1 and Definition 2.2 are equivalent. The stationarity condition for this , although this is not required in the simulation algorithm we developed later. We provide the expectation and variance of λ t and expectation of N t conditional on λ 0 in Proposition 2.3, which will be used later in Section 4 for numerically validating our algorithm. Fundamental distributional properties including Proposition 2.3, Laplace transform of λ t , and probability generating function of N t and the size of clusters can be easily derived by using formulas in Dassios and Zhao (2011) for a more generalised point process (so-called dynamic contagion process), also see further extensions in Dassios and Zhao (2013a) and Dassios and Zhao (2013b).

Proposition 2.3
The expectation and variance of λ t conditional on λ 0 are given by (y); and the expectation of N t is given by

Simulation Algorithm
The algorithm of exact simulation is given by Algorithm 3.1 for a univariate Hawkes process with exponentially decaying intensity (defined by Definition 2.1 or Definition 2.2) and random sizes of self-excited jumps in the intensity process. This algorithm is very easy to implement, because each of the random interarrival-times of jumps in the Hawkes process can be simulated exactly by decomposing it into two independent and simpler random variables without inverting the underlying cumulative distribution function. In addition, simulating intensity processes does not require stationarity conditions both for the cases of one-dimension in Algorithm 3.1 and multi-dimension in Algorithm 5.1.

Algorithm 3.1 (Univariate)
The simulation algorithm for one sample path of one-dimensional Hawkes process with exponentially decaying intensity ¦ (N t , λ t ) © t≥0 conditional on λ 0 and N 0 = 0, with jump-size distribution Y ∼ G andK jump-times {T 1 , T 2 , ..., TK}: 3. Record the (k + 1) th jump-time T k+1 in the intensity process λ t by 4. Record the change at the jump-time T k+1 in the intensity process λ t by where 5. Record the change at the jump-time T k+1 in the point process N t by Proof. Given the k th jump-time T k , the point process has the intensity process {λ t } T k ≤t<T k +S k+1 following the ODE with the initial condition λ t t=T k = λ T k . Obviously, (3) has the unique solution and the cumulative distribution function of the (k + 1) th interarrival-time S k+1 is given by By the inverse transformation method, we have However, we can avoid inverting this function F S k+1 (·) of (4) by decomposing S k+1 into two simpler and independent random variables S (1) • For the simulation of S (1) k+1 , since We can invert this simple function explicitly by Note that, S k+1 is a defective random variable as

and the condition for simulating a valid S
(1) • For the simulation of S Hence, for the simulation of S k+1 , we have where S (1) k+1 and S (2) k+1 are given by (5) and (6), respectively. Therefore, the (k + 1) th jump-time T k+1 in the Hawkes process is given by T k+1 = T k + S k+1 , and the change in λ t and N t at time T k+1 then can be easily derived as given by (1) and (2), respectively.
Algorithm 3.1 applies to any arbitrary distribution assumption G for self-excited jumpsizes {Y k } k=1,2,... , and of course also to the case when jump-sizes are fixed which is more widely used in the existing literature. By slightly adjusting the algorithm, it is straightforward to simulate process starting from any time with any arbitrary fixed value or distribution for λ 0 , and this flexibility is useful for practical implementations, such as applications in modelling credit risk and insurance risk, see Dassios and Zhao (2011) and Dassios and Zhao (2012). This algorithm can generate non-stationary process, and also is flexible to extend to simulate the process with stationary intensity by additionally assuming that λ 0 follows the stationary distribution 4 . For instance, if jump-sizes follow exponential distribution of parameter β, to simulate a stationary process, we only need set the initial intensity by

One Simulated Sample Path of Hawkes process
For instance, one simulated sample path of Hawkes process (N t , λ t ) with time T = 100 is provided in Figure 2. We can observe the clustering arrival of points of the Hawkes process by plotting the histogram of N t , and this contagion effect cannot be captured in classical Poisson models. This would be useful, for instance, for modelling the contagion risk of clustering defaults during economic crisis. For comparison, the theoretical expectations E[λ t |λ 0 ] and E[N t |λ 0 ] (given by Proposition 2.3) are also represented.

Comparison between Theoretical Formulas and Simulation Results
To analysis (in error percentage with respect to the exact theoretical result) are represented in Figure 3 and Table 1. Every point (marked by a star * ) in Figure 3 is calculated based on 100,000 simulated sample paths of (N t , λ t ).

Multi-dimensional Hawkes Process
By slightly modifying Algorithm 3.1, our method is straightforward to extend to multidimensional cases, for instance,D−dimendional point process with the underlying intensity process  2. Simulate the (k + 1) th interarrival-time W k+1 by , and each S [] k+1 can be simulated in the same way as S k+1 as given by Step 2 of Algorithm 3.1.

Numerical Example: A Bivariate Hawkes Process
Bivariate point processes are widely used in practice, for instance, for modelling the arrivals of (well-defined) events such as transactions, quote updates or price changes observable in the market, see Engle and Lunde (2003) and Bowsher (2007  Here, β [1,1] , β [2,2] provide measurement for self-contagion effect, and β [1,2] , β [2,1] provide measurement for cross-contagion effect. A simulated sample path with T = 100 is represented in Figure 4 by using Algorithm 5.1.