Perfect Sampling of Hawkes Processes and Queues with Hawkes Arrivals

In this paper we develop the first perfect sampling algorithm for queues with Hawkes input, i.e. single-server queues with Hawkes arrivals and i.i.d. service times of general distribution. In addition to the stability condition, we also assume the excitation function of the Hawkes process has a light tail and the service time has finite moment generating function in the neighborhood of the origin. In this procedure, we also propose a new perfect sampling algorithm for Hawkes process with improved computational efficiency compared to the existing algorithm. Theoretical analysis and numerical tests on the algorithms' correctness and efficiency are also included.


Introduction
Many stochastic systems have arrival processes that exhibit clustering or self-exciting behavior, i.e. an arrival will increase the possibility of new arrivals.As a natural extension of the classic Poisson process, Hawkes processes are used widely to model arrivals with selfexcitement.Examples include order flows in stock market ( [1], [11]), risk events in financial systems ( [4], [15]) and social network events ( [21], [24]).
To study the impact of self-excitement on the performance of stochastic systems, several papers have analyzed queueing dynamics with customer arrivals following a Hawkes, or other type of self-exiting process.[16] studies the heavy-traffic limit of infinite-server queue with Hawkes arrivals.[12] and [19] provide analytic solution to the transient and steadystate moments for different infinite-server systems with Hawkes arrivals.In particular, [12] studies the systems with Markovian Hawkes arrivals and phase-type/deterministic service times, while [19] studies the cases with non-Markovian Hawkes arrivals and exponential service times.In addition, [18] studies an infinite-server queue with shot-noise arrivals.[13] proposes a so-called Queue-Hawkes model that combining Hawkes process with an infiniteserver queue to capture ephemeral self-exciting behaviors.To the best of our knowledge, analytic results on queueing processes with Hawkes arrivals are only available for infiniteserver systems in the literature.Due to the dependence between customer arrivals and sojourn times, it is difficult to obtain analytic results even for the most simple single-server queue with Hawkes arrivals (see, for instance, the discussion on page 941 of [19]).
In this paper, we apply simulation techniques to numerically compute the steady state of queues with Hawkes input.In detail, we develop the a perfect sampling algorithm that can generate i.i.d.samples exactly from the steady-state distribution of Hawkes/GI/1 queues.Our algorithm is applicable to a variety of queueing models with Hawkes arrivals.In detail, we assume the arrival process is a linear Hawkes process, which covers both Markovian and non-Markovian Hawkes process as studied in [19] and [13], and the service times are i.i.d.following a general continuous distribution.
Our algorithm is closely related to the literature on perfect sampling of queueing models, for instance [5], [6], [7], [9], [14], and [23], to name but a few.Most of the existing works have been focused on queues with arrivals modeled as Poisson or renewal processes.In those cases, the steady-state waiting time can be related to the maximum or the running maximum of a (possibly multi-dimensional) random walk, see the discussion on page 378 of [6].The running maximum of random walk is defined as the maximum from a positive time n to infinity, for any n.In our case, the steady-state waiting time of Hawkes/GI/1 queue can still be represented as max n≥0 R(n) for certain stochastic process R(n).However, R(n) is defined by a stationary version of the Hawkes arrivals and its increments have sequential dependence.Therefore, existing perfect sampling algorithms for queueing models can not be directly applied to queues with Hawkes arrivals.
There are two key steps in our algorithm.First, to simulate the process R(n) involves generating a stationary sample path of the Hawkes arrivals, namely, perfect sampling of Hawkes process, and therefore is far from trivial.The only existing perfect sampling algorithm for Hawkes processes in the literature is developed in [22].However, the complexity of the algorithm, in terms of the expected total number of random seeds generated, is infinite (see Propostion 2).Using importance sampling techniques, we proposed a new, and probably the first perfect sampling algorithm for Hawkes processes that has finite expected termination time.In particular, the complexity of our new algorithm is finite and has an explicit expression in terms of model and algorithm parameters.We believe this new algorithm can be applied to, in addition to the single-server queue studied in the current paper, other stochastic models that involve Hawkes processes as we mentioned previously.
Once the stationary Hawkes process and R(n) are simulated, the next key step of our algorithm is to find out max n≥0 R(n).The main idea is to construct a random walk R(m) coupled with the Hawkes/GI/1 queue such that its running maximum max m≥n R(m) dominates the process R(n) in a proper way.Then, we apply the techniques dealing with the running maximum of random walks, as developed in [6], to simulate max n≥0 R(n) jointly with max m≥n R(m), which completes our algorithm as max n≥0 R(n) equals in distribution to the steady-state waiting time.
The rest in the paper is organized as follows.We first introduce the definition of Hawkes process, the Hawkes/GI/1 queue model and the technical assumptions in Section 2. In Section 3, we introduce our perfect sampling algorithms for Hawkes processes and Hawkes/GI/1 queues, along with the main results of the paper.In Section 4, we implement the algorithms and report the numerical experiment results.Finally, Section 5 concludes the paper with a brief discussion on future research directions.Proofs of some technical results are included in the Appendices.

Hawkes Process
Following [19], we provide two equivalent definitions for Hawkes process, namely the conditional intensity and cluster representation definitions.The conditional intensity definition clearly demonstrates self-excitement of Hawkes arrivals while the cluster representation provides an alternative probabilistic construction of Hawkes process, which will be used in our simulation algorithm.Definition 1. (Conditional Intensity) A Hawkes process is a counting process N (t) that satisfies as ∆t → 0, where F(t) is the associated filtration and λ is called the conditional intensity such that where t 1 , t 2 ,... are the arrival times, the constant λ 0 > 0 is called the background intensity and the function h : R + → R + is called excitation function.
According to (1), an arrival will increase the future intensity function and the possibility of new arrivals, so arrivals of a Hawkes process are self-exciting.Now we introduce another equivalent definition (or construction) of Hawkes process which represents it as a branching process with immigration.Definition 2. (Cluster Representation) Consider a (possibly infinite) T ≥ 0 and define a sequence of events {t n ≤ T } according to the following procedure: 1.A set of immigrant events {τ m ≤ T } arrive according to a Poisson process with rate λ 0 on [0, T ].

2.
For each immigrant event τ m , define a cluster C m , which is a set of events, as follows.
Each event ∈ C m is indexed by k ≥ 1 and is represented by a tuple (k, t k m , pa k m ), where t k m is the event's arrival time and pa k m ≥ 0 is the index of its parent event.Following this representation, we denote the immigrant event as (1, τ m , 0).Then, the counting process N (t) corresponding to the event sequence {t n } is equivalent to the Hawkes process defined by conditional intensity function (1).

The cluster
For each cluster C m , we define a non-decreasing function S m (t) = |{t n m : τ m ≤ t n m ≤ t + τ m }|, i.e. the number of events in C m that arrive on interval [τ m , τ m + t] for t ∈ R. By definition, for each m, S m (t) is an increasing function such that S m (t) ≡ 0 for t < 0, S m (0) = 1 and S m (∞) = |C m |.Then, the Hawkes counting process N (t) = ∞ m=1 S m (t−τ m ).We define for each cluster C m its cluster length L m as length of the time interval between the immigrant event and the last event in this cluster, i.e.
We call τ m and δ m max k t k m the arrival time and departure time of the cluster C m , respectively. Let Then, according to Definition 2, h 1 is the expected number of nextgeneration events generated by each single event.For any event t k m that is not an immigrant, the arrival time of its parent event is t pa k m m following our notation.We define the birth time Following the property of non-homogeneous Poisson process, conditional on

positive random variables that follow the probability density function f (t)
h(t)/h 1 for t ≥ 0. Given the clear meaning of h 1 and f (•) in the cluster representation of Hawkes process, in the rest of the paper, we shall denote by (λ 0 , h 1 , f (•)) as the parameters that decide the distribution of a Hawkes process.

Stationary Hawkes Process
For the Hawkes process with parameters (λ 0 , h 1 , f (•)) to be stable in long term, intuitively, each cluster should contain a finite number of events on average.Therefore, we shall impose the following stability condition on the Hawkes process throughout the paper, which is also a common assumptions in the literature ( [10], [17]): Under this condition, the Hawkes process has a unique stationary distribution ( [10]).Actually, we can directly construct a stationary Hawkes process using the cluster representation as follows.Note that the arrival process of the immigrant, or equivalently, the clusters, is a homogeneous Poisson process and can be extended to time interval (−∞, ∞).For this two-ended Poisson process, we index the sequence of immigrant arrival times by {±1, ±2, ...} such that τ −1 ≤ 0 < τ 1 and generate the clusters {C ±m:m=1,2,... } independently for each m following the procedure in Definition 2.Then, the events that arrive after time 0 form a stationary sample path of a Hawkes process on [0, ∞), namely is a stationary Hawkes process.

Hawkes/GI/1 Queue
We consider a single-server queue where customers arrive according to a stationary Hawkes process with parameters (λ 0 , h 1 , f (•)).We denote by {U n } the corresponding sequence of inter-arrival times.Upon arrival, customers are served FIFO (first-in-first-out) and their service times V n are i.i.d.positive random variables of general distribution with probability density function g(•).
Denote by W (t) the virtual waiting time process (also known as the workload process) of this single-server queue.Mathematically, W (t) can be defined as a reflected process as follows.Let R(t) In the setting of a single-server queue, the function L(t) is nothing but the server's idle time by time t.A single-server queue is said to be stable if the distribution of W (t) converges as t → ∞.The following proposition states that the Hawkes/GI/1 queue is stable if and only if its service rate is higher than the stationary arrival rate of customers.Proposition 1.The Hawkes/GI/1 queue is stable if and only if Proof of Proposition 1.The stability condition directly follows Propositions 1.1 and 1.2 on page 267 of [2].For a single-server queue with stationary inter-arrival time sequence {U n } and service time sequence For the Hawkes/GI/1 queue, each arrival of the Hawkes process is associated with a service time.Therefore, we shall also include service-time information in the cluster representation of Hawkes process.In detail, we denote m is the service time of the k-th customer in cluster C m .Our goal is to simulate the steady-state virtual waiting time for the Hawkes/GI/1 queue, so we shall impose the stability condition throughout the paper.Besides, in order to carry out importance sampling procedures in the simulation algorithm, we also need some technical assumptions on the Hawkes/GI/1 model.These assumptions are satisfied by a large class of queueing models, such as Hawkes/GI/1 queue with arrivals that are Markovian Hawkes, or non-Markovian Hawkes having excitation function with finite support, and service times that are exponential or phase-time.Now we close this section by summarizing our assumptions on the Hawkes/GI/1 queue.Assumption 1.The stability condition (4) holds.Assumption 2. The birth time of the Hawkes process b and the service time of customers V are of continuous distribution.Besides, there exists θ 0 > 0 such that i.e. the random variables b and V have finite moment generating function in a neighborhood of the origin.Assumption 3. The distributions of birth time and service time can be simulated exactly, and moreover, we can simulate from exponential tiltings (i.e., the natural exponential family) associated with these distributions.

Simulation Algorithms
Our goal is to simulate the steady-state virtual waiting time W ∞ of Hawkes/GI/1 queue.We now construct an expression of W ∞ in the same spirit of [20].First, we extend the stationary Hawkes arrival process N (t) backward in time to (−∞], and for t < 0, define N (t) as the number of arrivals on [t, 0].The i.i.d.service time sequence {V n } can be natural extended to n ≤ −1.Following the same argument as in Proposition 1 of [5], we can construct the stationary distribution W ∞ as According to (5), we can simulate W ∞ in two steps 1. Generate the sample path of R(t).

Simulate the maximum max t≥0 R(t).
To simulate the process R(t) in Step 1 is essentially to simulate a stationary sample path of the Hawkes process backward in time.We shall explain how to do this in Section 3.1, in which we develop a novel and efficient perfect sampling algorithm for Hawkes process.Then we explain Step 2 in Section 3.2 and this completes our algorithm.

Perfect Sampling of Hawkes Process
To the best of our knowledge, the only perfect sampling algorithm for Hawkes process in the literature is given by [22], which is based on the cluster representation of Hawkes process.However, this algorithm is not very efficient, and it on average needs to generate an infinite number of random numbers before termination (as we shall prove in Proposition 2).In this part, we propose a novel perfect sampling algorithm for Hawkes process, whose complexity (in terms of the expected variables generated) is finite and has an explicit expression in model and algorithm parameters.Our algorithm is also based on the cluster representation, but exploits importance sampling and acceptance-rejection techniques to largely improve the simulation efficiency.

Existing Framework
As both our algorithm and [22] are based on cluster representation of Hawkes processes, we first briefly review the algorithm in [22] to introduce the outline of our algorithm.Then, we point out the major bottlenecks in their algorithm, and in Section 3.1.2,we provide our solutions to these bottlenecks and introduce the new algorithm.
Recall that the arrival process of the immigrant events is a homogeneous Poisson and can be extended into a double-ended counting process {τ m } ∞ m=−∞ .Let M (t) = max{m : τ m ≤ t} be the index of the last immigrant that arrives by time t.The counting process N (•) corresponding to a stationary Hawkes process on [0, ∞) can be decomposed as, for any t > 0, Here L m is the cluster length as defined in (2).The second equality just says that if a cluster C m that arrives before time 0 has cluster length L m < −τ m (minus its arrival time), all of its events will arrived before 0, and therefore, it will have no impact on dynamic of the stationary Hawkes process after time 0. In the above decomposition, N 1 (•) is the set of events from clusters arriving after time 0, which is basically a Poisson compound of i.i.d.clusters.To simulate N 1 (•), we just need to simulate the arrivals of clusters according to a Poisson process of rate λ 0 , and then, for each arrival, simulate a cluster independently according to Step 3 as described in Definition 2. Therefore, the key step in perfect sampling of Hawkes process is to simulate N 0 (•), or equivalently, to simulate all the clusters that have arrived before time 0 and last after time 0. By slightly notation abusing, we shall also use N 0 to denote the set of clusters that have arrived before time 0 and last after time 0, i.e., be the probability that a cluster last for more than t units of time.Then, by Poisson thinning theorem, the clusters in N 0 arrive on time horizon (−∞, 0] according a non-homogeneous Poisson process with intensity function γ(t) = λ 0 p(−t), for t ≤ 0. Based on this observation, [22] proposes the following procedure to simulate N 0 : Step 1: Sample a non-homogeneous Poisson process with intensity function γ(t) on (−∞, 0] and obtain the arrivals {τ −1 , ..., τ −K }.
Step 2: For each k ∈ {1, 2, ..., K}, repeatedly simulate a cluster sample C until its cluster length The above algorithm has two major bottlenecks in the design that significantly affect its computational efficiency.First, in Step 1, the function p(•) does not have an explicit expression, so [22] uses a while loop to approximate p(•) by iteration.Second, in Step 2, a naive acceptance-rejection procedure is used to obtain a cluster with cluster length > −τ k .As a consequence, the total number of cluster simulation rounds, before an acceptance occurs, could be very large when −τ k is large.The following proposition shows that in fact, Step 2 needs to spend, on average, infinite rounds of cluster simulation before termination.Proposition 2. Let N C be the total number of clusters generated in Step 2 by naive acceptancerejection.Then, Proof of Proposition 2. For any fixed t > 0, the expected number of clusters generated to obtain one sample with cluster length L > t equals to P (L > t) −1 = p(t) −1 .Therefore,

New Algorithm with Improved Efficiency
The bottleneck in Step 2 is indeed a rare event simulation problem, namely, for a given t > 0, to simulate C m conditional on the event {L m > t} which could have very small probability for large t.In our new algorithm, we use importance sampling to design a more efficient procedure to sample from the conditional distribution of C m given L m > t.Besides, we also combine importance sampling with a more sophisticated acceptance-rejection procedure to avoid evaluating p(t) in Step 1, thus further improve the computational efficiency.
We shall apply exponential tilting to do the importance sampling in Step 2. However, exponential tilting with respect to the cluster length L m is not easy.Instead, we shall do exponential tilting with respect to the total birth time where b k m is the birth time of event k in cluster C m as defined in (3), and b 1 m = 0 for the immigrant event.There are two reasons for us to use B m to do the exponential tilting.First, as the total birth time is always larger than cluster length, for any cluster that arrives at time −t to last after time 0, it must have total birth time > t.Therefore, to simulate N 0 , it is sufficient to find the set of clusters {C m : m ≤ −1, B m > −τ m }.Second, to apply importance sampling to a cluster C m using exponential tilting with respect to total birth time B m is much easier.In Proposition 3 below, we summarize the properties of the random variable B m that are useful for our simulation algorithm, and explain explicitly how to do exponential tilting with respect to B m .Proposition 3. Consider a cluster C m with parameter (h 1 , f (•)).Let L m and B m be its cluster length and total birth time, respectively.The following statements are true: (1) B m ≥ L m .
(3) Let P be the probability distribution of a cluster C m .Let Q be the importance distribution of the cluster under exponential tilting by parameter 0 < η < θ 1 with respect to the total birth time B m , i.e.

dQ(C
Then, sampling a cluster from the importance distribution Q is equivalent to sampling a cluster with parameter (h Given Proposition 3, we are ready to give the whole procedure of our simulation algorithm for N 0 as described in Algorithm 1.The proof of Proposition 3 is given in Appendix A.  5.Return N 0 . Algorithm 1 contains two importance-sampling steps.In Step 2, instead of simulating the arrivals of clusters following the non-homogeneous Poisson with intensity function γ(t), it simulates a non-homogeneous Poisson with larger intensity function γ(t) ≥ γ(t) (by Markov's Inequality).In Step 4, it applies importance sampling to generate the conditional distribution of clusters.In the end, it utilizes one step of acceptance-rejection to transform the two importance sampling probability laws jointly into the target distribution.
The first main result of the paper is stated as Theorem 1 below, in which we provide a theoretical guarantee that the output of Algorithm 1 follows exactly the distribution of N 0 , and an explicit expression of algorithm complexity in terms of model and algorithm parameters.
Theorem 1.The list of clusters generated by Algorithm 1 exactly follows the distribution of N 0 .In particular, (1) The arrival times of the clusters follow a non-homogeneous Poisson process with intensity γ(t) = λ 0 p(−t) for t ∈ (−∞, 0]. (2) For each cluster C m in the list, given its arrival time τ m , it follows the conditional distribution of a cluster given that the cluster length > −τ m .
Besides, the expected total number of random variables generated by Algorithm 1 before termination is Proof of Theorem 1.To prove Statement (1), by Poisson thinning theorem, it suffices to show that, for each m, the acceptance probability of the cluster C m equals to According to Proposition 3, the importance distribution Q and the target distribution P satisfies .
Therefore, the probability for cluster C m to be accepted in Step 3 is Therefore, we obtain Statement (1).
From the above calculation, we can also see that, given m and τ m , and any event A ∈ σ(C m ), the joint probability Therefore, the accepted sample of C m indeed follows the conditional distribution of C m given {L m > −τ m }, and we obtain Statement (2).
To check (8), we first note that the expected total number of random variables generated by Algorithm 1 is equal to the expected total number of clusters multiplied by the average number of random variables generated in one cluster.In Step 2, the number of clusters generated is a Poisson random variable with mean ∞ 0 γ(−t)dt = λ 0 exp(ψ B (η))/η.The average number of events in each cluster is , where the last equality follows from (7).Besides, for each cluster, the algorithm also need to simulate an extra random number in the acceptance-rejection step.Therefore, the expected total number of random variables is which closes the proof.
Remark 1.The complexity result (8) not only guarantees that Algorithm 1 terminates in finite time in expectation, it also provides some guidance to the optimal choice of η that reduces the computational cost.
Reversing the Time So far we have focused on simulating a stationary Hawkes process forward in time.But to simulate the steady-state waiting time W ∞ following (5), we need to generate the sample path of stationary Hawkes process backward in time on the interval (−∞, 0].We now explain briefly how to do this using Algorithm 1.
Recall that τ m and δ m are the arrival time and departure time of cluster C m as defined in Section 2.1.By definition, the cluster length L m = δ m − τ m .Since L m are i.i.d.distributed and independent of τ m , by Poisson thinning theorem, {δ m } also follows a homogeneous Poisson process of rate λ 0 just as {τ k }.Therefore, to simulate the Hawkes process backward in time, we can first apply Algorithm 1 to simulate N 0 , i.e. list of clusters that depart after time 0 but arrive before time 0, and then simulate those clusters that depart before time 0 according to a Poisson process with rate λ 0 .In detail, we can use the following procedure to simulate a stationary Hawkes process backward in time on [−t, 0] for any t > 0: 1.Call Algorithm 1 to simulate N 0 .
3. For each δ m , simulate a cluster of events following Step 3 in Definition 2. Adjust the event times accordingly such that the arrival time of the last event equals to δ m .
Once a stationary Hawkes process can be simulated backward in time, we can simulate the process R(t) in ( 5) for any t > 0. The next step is to find out max t≥0 R(t).

Perfect Sampling of Single-Server Queue with Hawkes Arrivals
In this part, we explain the second key step in our algorithm, namely, to simulate max t≥0 R(t).
According to (5), once max t≥0 R(t) is simulated, we just return its value as an exact sample of W ∞ .
Recall that {U n : n ≤ −1} is the sequence of inter-arrival times of the Hawkes process.Then, for any t ≥ 0, because the left side is the arrival time of the first customer after time −t.Besides, the equality holds when t is the arrival time of an event.Therefore, With slightly notation abusing, let's denote R(n) For a GI/GI/1 queue, the corresponding {R(n) : n ≥ 1} is basically a random walk with negative drift.For example, in [14], the perfect sampling algorithm for GI/GI/1 is presented as a direct application of the simulation algorithm for the maximum of random walks with negative drift.In our case, however, due to the self-exciting behavior of arrivals, U n has sequential dependence and as a result, the distribution of R(n) is more complicated than a random walk.We shall deal with this issue by constructing an auxiliary random walk R(m) coupled with the Hawkes process, such that R(m) has negative drift and its running maximum max m≥n R(n) "dominates" the process R(n) in a proper way.Then, we can bound and learn the exact value of max n≥0 R(n) by simulating the running maximum max m≥n R(m), applying the techniques developed in [6].

The Auxiliary Random Walk
We first explain our construction of the random walk R(m).For a stationary Hawkes process backward in time, let's index its clusters in the order of their departure times.(In our previous notation, we index clusters in the order of their arrival times.)In detail, we shall denote the m-th cluster that depart before time 0 as cluster C −m and its departure time as δ −m < 0. Similarly, clusters that depart after time 0 are also indexed by positive integers in the order of their departure times.We now denote by A −k − k i=1 U −i as the arrival time of the k-th customer before time 0. By definition, δ −m is the arrival time of the last customer in cluster C −m .As a result, for each m ≥ 1, there must exist k m ≥ 0 such that We define the auxiliary random walk R(m) as the total service requests in cluster C −1 , ..., C −m minus the units of time that elapsed, i.e.

R(m)
Besides, for each m we denote by J(m) as the total service requests of customers who arrive before time δ −m and belong to clusters that depart after time δ −m , i.e.
The next proposition shows that the auxiliary random walk R(m) has negative drift, and its increment "dominates" the increment of R(k) in a certain sense.Its proof is given in Appendix B. (2) For all m 2 ≥ m 1 ≥ 1, and any

Perfect Sampling Algorithm for W ∞
The following corollary is a direct consequence of Proposition 4, and provides a way to find max k≥0 R(k) by simulating the running maximum of the auxiliary random walk R(m).
Corollary 1. Suppose there exist m 2 ≥ m 1 ≥ 1 such that the following statements are true: Then, we can conclude Following Statement (2) in Proposition 4, we have Corollary 1 implies that we can stop simulating the stationary Hawkes process (backward in time) and return W ∞ = max 0≤k≤km 2 R(k) when a pair of random times m 1 and m 2 satisfying Conditions (a) and (b) are detected.For given m 1 , and the Hawkes process up to cluster C −m , it is straightforward to check whether m satisfies Condition (a).Besides, since R(m) is a random walk with strictly negative drift, the smallest m 2 that satisfying Condition (a) is finite.Condition (b) involves checking whether the running maximum of R(m) exceeds a given level.Based on the above observations, we provide the outline of our perfect sampling algorithm for W ∞ d = max k≥0 R(k) as describe in Algorithm 2 In Step 4 of Algorithm 2, we directly apply the importance sampling technique in [6] to simulate the Bernoulli random variable jointly with the auxiliary random walk and to check Condition (b).To do this, we first need to verify that all assumptions in [6] are satisfied.Proposition 5.The following statements are true under Assumptions 1 and 2: Then, there exists θ 2 > 0 such that for all 0 ≤ θ < θ 2 ,

Numerical Experiments
We implement Algorithm 1 and Algorithm 2 in Python to test the performance and correctness of our perfect sampling algorithms.As an example of algorithm application, we investigate the effect of self-excitement on the steady-state waiting time distribution numerically using our perfect sampling algorithms.
Algorithm 1 Performance Test We consider a Hawkes process with parameters λ 0 = 1, h 1 = 0.5, f (t) = 2 exp(−2t).Then, the stationary intensity rate of this Hawkes process is λ 0 /(1 − h 1 ) = 2.To test the correctness of our algorithm, we apply Algorithm 1 to simulate a stationary Hawkes process on time interval [0, 1].If our algorithm is correct, the average number of events generated on this time interval should equal to E[N (1)] = 2, regardless of the choice of algorithm parameter η.To illustrate algorithm efficiency, we also record the number of random variables generated in each simulation round.In Table 1, we report, for different values of η, the 95% confidence interval for E[N (1)], along with the average number of random variables generated, based on 10000 rounds of simulation.Algorithm 2 Performance Test We consider a Haweks/GI/1 queue in which customers arrive according to a Hawkes process with parameters λ 0 = 1, h 1 = 0.5, f (t) = 2 exp(−2t), and the service times are i.i.d.exponential with mean 1/3.Since there is no theoretical benchmark to compare with, we shall compare the empirical distribution of the samples generated by Algorithm 2 with that of the samples generated by simulating Hawkes/GI/1 for a long time such that the system is close to its steady sate.In figure 1, we plot the histograms of 10000 i.i.d.samples generated by Algorithm 2 and 10000 i.i.d.samples of W (100) obtained by simulating the single-server queue from empty state, i.e.W (0) = 0.When simulating W (100), to mitigate the transient bias caused by Hawkes arrivals, we feed in the single-sever queue with stationary sample paths of Hawkes process generated by Algorithm 1.
In order to compare the efficiency of our algorithm with naive simulation, in each simulation round of Algorithm 2, we also record T ps = −τ −m 2 +∆ which is minus of the departure time of the last cluster generated by the algorithm before termination.Basically, T ps is length of the sample path generated by Algorithm 2 and we call T ps the PS sample path length.In this light, E[T ps ] can be used as a measurement for the computational cost of Algorithm 2. To evaluate the efficiency of naive simulation, we estimate the mixing time   We conjecture that, this is because, the perfect sampling algorithm can terminate much earlier before the time point T when the system get close to the steady state on average, say T = 40, if the simulated sample path coalesces with the steady state, i.e. when a Bernoulli B = 0 is detected in Step 4 of Algorithm 2, very soon.We plot the empirical distribution of T ps in Figure 3(b) and find that more than 50% of the sample paths are shorter than 20.On the other hand, a significant proportion (over 10%) are longer than 50.Intuitively, these sample paths could contribute to the average transient bias of E[W (T )] for all T < 50 and thus affect the accuracy of naive simulation with fixed T < 50.
Impact of Self-excitement on Waiting Time Now we apply our simulation algorithms to estimate the impact of self-excitement on the steady-state distribution of waiting times.In particular, we investigate a set of Hawkes/GI/1 queues with equal stationary customer arrival rate and same service time distribution, but different levels of self-excitement in customer arrivals.To see the impact of self-excitement on the steady-state distribution of the virtual waiting time, we apply Algorithm 2 to estimate E[W ∞ ] and V ar(W ∞ ).
In detail, we consider 5 Hawkes/GI/1 queues indexed by i = 1, 2, .., 5.The parameter set of the Hawkes process in queue i is (λ i 0 , h i 1 , f i (•)).We set f i (t) = 2 exp(−2t) for all i and h i 1 ∈ {0.3, 0.4, 0.5, 0.6, 0.7}.For each i, λ i 0 satisfies λ i 0 /(1−h i 0 ) ≡ 2 so that all 5 Hawkes/GI/1 queues have equal stationary customer arrival rate.We also assume that the service times are exponential with rate µ = 3 in all queues.For each queue, we run 10000 rounds of simulation using Algorithm 2 and report the estimated E[W ∞ ] and V ar(W ∞ ) in Table 2. From the simulation results, we can see that self-excitement behavior of customer arrivals could increase not only the mean of waiting time but also its level of dispersion measured by the variance-to-mean ratio.

Conclusion
In this paper, we develop a simulation algorithm to generate i.i.d.samples exactly from the steady state of Hawkes/GI/1 queues.To the best of our knowledge, this is the first perfect sampling algorithm for queueing models with arrivals that have self-exciting behavior.As a key component fo the algorithm, we also develop a new perfect sampling algorithm for Hawkes process that is much more efficient compared to existing algorithm in the literature.Both algorithms utilize importance sampling techniques on the cluster representation of Hawkes process.This approach, we believe, can probably be extended to other classes of counting processes with cluster representation, such as multi-dimensional Hawkes process with mutual-excitement and Poisson cluster processes.Besides, as single-server queues are the basic building blocks of more complicated queueing models, our approach can probably lead to perfect sampling methods for other queueing models with self-exciting arrivals.

A Proof of Proposition 3
We prove the three statements of Proposition 3 one by one.
(1).To see B m ≥ L m , probably the most straightforward way is to represent the cluster as a tree.Let the immigrant event τ m be the root node and link each event t k m to its parent event by an edge of length b k m .Then, by their definitions, B m is equal to the total length of all edges in the tree while L m is equal to the length of the longest path(s) from the root node to a leaf node.Therefore, B m ≥ L m .
( It is known in literature ( [8]) that the c.g.f.ψ S (θ) of S m is well-defined in a neighborhood around 0. Under Assumption 2, the c.g.f. of each b k m is well-defined on [0, θ 0 ].Since B m is the compound sum of (S m − 1) i.i.d.birth times, ψ B (θ) is also well-defined in a neighborhood around 0.
To obtain (7), recall that τ m = t 1 m is the immigrant event and suppose its next-generation events are t 2 m , t 3 m ,..., t Λ+1 m .Then, Λ is a Poisson r.v. with mean h 1 .By the self-similarity of branching process, the total birth time B m equals to Λ k=1 (t k+1 m − t 1 m ) plus i.i.d.copies of total birth times B Let K (k) be i.i.d.copies of K for k ≥ 1 and Λ be the number of next-generation events generated by the immigrant.Then, by self-similarity of the branching process, where ψ V is the c.g.f. of the service time.Therefore, and as a result, ψ R,η (θ) = ψ V,η (θ) + h 1 exp(ψ K (η))(exp(ψ K,η (θ)) − 1) − log λ 0 + η λ 0 + η + θ .

Figure 1 :
Figure 1: Comparison of the estimated distributions of W ∞ from 10000 i.i.d.samples generated by perfect sampling Algorithm 2 and long-term simulation respectively.

Figure 2 :
Figure 2: (a) Comparison between the mixing time of Hawkes/GI/1 queue and the average sample path length generated by Algorithm 2. (b) Empirical distribution of sample path length generated by Algorithm 2.
(a) mixing time v.s.PS sample path length (b) PS sample path length distribution ). Recall that B m = |Cm| k=2 b k m and b k m are i.i.d. for given |C m |.Let S m = |C m |.
Update k ← k + 1 and add the newly generated events into set C m .Repeat the above iteration until no more events are generated.

Table 1 :
95% confidence interval for E[N (1)] and expected number of random variables generated in one simulation round, estimated from 10000 i.i.d.sample path of Hawkes process on [0, 1] simulated by Algorithm 1 with different η.

Table 2 :
Estimated mean, variance and VMR (variance-to-mean ratio) of steady-state virtual waiting time for Hawkes/GI/1 queues with different level of self-excitement.