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.