Perfect sampling of a single-server queue with periodic Poisson arrivals

In this paper we present algorithms for the perfect sampling of single-server time-varying queues with periodic Poisson arrivals under the first come first served (FCFS) discipline. The service durations have periodically time-dependent exponential (Mt/Mt/1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm M _t/\mathrm M _t/1$$\end{document}) or homogeneous general (Mt/G/1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm M _t/\mathrm G /1$$\end{document}) distributions. Assuming a cycle length of 1, we construct discrete dominating processes at the integer instants n∈{0,±1,…}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \in \{0, \pm 1, \ldots \}$$\end{document}. Perfect sampling of the Mt/Mt/1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm M _t/\mathrm M _t/1$$\end{document} queue is obtained using dominated CFTP (Kendall and Møller 2000) when the system is relatively lightly loaded or with the regenerative method (Sigman 2012) in the general case. For the Mt/G/1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm M _t/\mathrm G /1$$\end{document} queue, perfect sampling is achieved with dominated CFTP.


Introduction
Time-varying queueing models are more realistic than time-homogeneous queues, but they are not usually mathematically tractable [17, p. 697]. As noted by [14], computational methods and approximation techniques involved in time-varying queueing problems have long been regarded as challenging.
In this paper we consider cases in which the time-dependent stochastic processes follow periodic patterns. Hasofer [10] showed that the Laplace-Stieltjes Transform (LST) of virtual waiting time in the M t /G/1 queue is asymptotically periodic in time.
Harrison and Lemoine [9] proved that the virtual waiting time at any given time has its own limiting distribution, and there is one such distribution for each point within the period of the queue. Asmussen and Thorisson [4] extended the context to more general cases, where the inter-arrival times and service durations both depend on the arrival instant within the time period. They proved that with more conditions (such as Harris ergodicity [1, p. 202] of the phase parameter which the inter-arrival time and service duration depend on), the virtual waiting time and queue length also have time-dependent limiting distributions in periodic patterns.
Due to the complexity of the time-varying systems, only asymptotic solutions have been developed, and this has happened gradually over recent decades. By assuming some state at time 0, [22] and [13] found the transient distributions of the number of customers in the system (Q t ) in M t /M t /1 and M t /M t /c queues, respectively, using generating functions and Volterra integral equations. Zeifman et al. [21] approximated the limiting mean value (E(Q t )) of the M t /M t /1 queue with the transient distribution of the truncated time-varying birth and death processes by restricting their difference to some controllable extent. The asymptotic periodic solutions for M t /M t /1 and M t /M t /c systems were achieved by Margolius [14], where distributions and moments are given in terms of integral equations.
If we could draw samples directly from the steady-state distribution, statistical inference would become straightforward. Perfect sampling is an approach to draw a sample directly from the steady-state distribution without explicitly solving for it. The first well-known perfect sampling algorithm is commonly referred to as coupling from the past (CFTP), introduced by [15]. "Dominated CFTP" [11] is an important extension of CFTP. It enables the coupling of Markov chains with unbounded state spaces by reducing the number of past chains that need to be simulated.
Recently, [18] applied dominated CFTP to achieve perfect sampling of an M/G/c queue with a "super stable" (i.e., ρ < 1/c) condition. The backward simulation of an M/G/1 queue was implemented by running the coupled M/G/1 Processor-sharing (PS) model, which is time reversible. Connor and Kendall [7] showed how to generalize the dominated CFTP idea used by [18] to relax the super stable restriction to ρ < 1. They further proposed a sandwiching dominated CFTP algorithm for the perfect sampling of the M/G/c queue. It significantly reduces the expected runtime. However, in the time-varying circumstance (M t /G/1), it is hard to construct a dominating process which empties, or whose upper and lower envelopes coalesce.
Other methods are also available for perfect sampling. In [2, p. 420] the perfect sampling of regenerative processes was described and then applied to the M/G/c queue by [19]. A special case is perfect sampling of the GI/G/1 queue [2, p. 437], assuming that the Exponential Change of Measure (ECM) [1, p. 352] exists for the underlying random walk.
In practice, dominating processes are key elements for perfect sampling. They are processes defined on the same probability space, which bound the states of the target process, reducing the range of unknown values to a bounded set. In this paper, dominating processes are constructed by modifying the arrival instants or potential departure instants on each periodic cycle of the time-varying queues. In the M t /M t /1 queue, we can simulate the steady-state draw of the dominating process using the ECM method mentioned above for the GI/G/1 queue. For the M t /G/1 queue, we estimate the upper bound of the dominating process based on a coupled homogeneous queue.
In Sect. 2, we present our assumptions and notation. Section 3 presents perfect sampling of the M t /M t /1 queue in both the relatively lightly loaded and more general cases. In Sect. 4, we use dominated CFTP to achieve perfect sampling for the M t /G/1 queue. Section 5 concludes the paper.

Assumptions and notation
To ensure consistency and clarity, we present our assumptions and notation which we will require in the subsequent analyses.
Let x + be the non-negative part of x; i.e., if x > 0, then x + equals x; otherwise it is 0. -All queueing systems involved in this paper are work conserving and non preemptive. Queueing disciplines are generally presumed to be First Come First Served (FCFS) unless described otherwise. -Arrivals to the queue form a time-varying Poisson process, with non-trivial periodic arrival rate λ(t) ≥ 0. Without loss of generality, assume the length of the period is 1. We also assume that "potential service events" form a time-varying Poisson process with periodic rate μ(t) ≥ 0, also with period 1. These are the departures that would occur if the system were to remain busy; there may be fewer actual departures if the queue empties. Both λ(t) > 0 and μ(t) > 0 except possibly at discrete points, so their integrals are strictly increasing. Time "0" is congruent with our time point of interest. For both λ(t) and μ(t), t is measured in clock time: λ(t) = λ(t + 1), and μ(t) = μ(t + 1), ∀t ∈ R. -Defineλ for t ∈ (0, 1]. These functions are strictly increasing on the defined interval, according to the definitions above of λ(t) and μ(t). Therefore their inverse functions exist, and are denoted by F −1 λ (x) and F −1 μ (x), x ∈ (0, 1], respectively. -To ensure the stability of the M t /M t /1 queue, the occupancy ρ must be less than unity; i.e., See [4,Theorem 4.5] for more details. For the M t /G/1 queue, we define μ = 1/E(B) where B is the homogeneous service duration, and for stability, -For the time-varying queues, let N A k be the number of arrivals on the interval The N A k 's constitute an i.i.d. sequence of random variables. -In the M t /M t /1 queue, let N D k be the number of potential departures on the interval The N D k 's also constitute an i.i.d. sequence of r.v.'s, and they are independent of the N A k 's. When the server is idle, potential departure events have no effect. -In the algorithms that follow we will couple homogeneous queues to the timevarying queues that we are studying. Denote by Q N t and Q H t the numbers of customers at time t in the time-varying and homogeneous queues, respectively. Similarly, let V N t and V H t be the unfinished workloads in these two systems, respectively. These processes are all right continuous in t.

Perfect sampling of the M t /M t /1 queue
In this section, we present perfect sampling for the M t /M t /1 queue using one of two methods, depending upon the occupancy level in the queue. For the lightly loaded case where the minimum service rate is greater than the maximum arrival rate, dominated CFTP works as a straightforward solution. In the general setting, where we only havē λ <μ, we achieve the perfect sampling using an ECM to sample from the GI/G/1 queue and the regenerative method to extend this to the time-varying queue.

Perfect sampling of the M t /M t /1 queue with inf μ(t) > sup λ(s)
Let μ l = inf μ(t) and λ u = sup λ(t). Assume A stable M/M/1 queue can be generated with arrival and service rates of λ u and μ l , respectively, since μ l > λ u . Based on its homogeneous arrival and potential departure events, the time-varying inputs of the coupled M t /M t /1 queue are simulated as follows: -The time-varying arrival events are filtered from the homogeneous arrivals using the thinning method [17, p. 697, Method 1], thinning arrivals to rate λ(t) ≤ λ u .
-Time-varying potential departure events are reproduced based on the homogeneous ones by supplementing the homogenous events with events from a time-varying Poisson process at rate μ(t) − μ l . These extra events are generated with the interevent time method [17, p. 702, Method 3]. Thus, the composite potential departure process is a superposition of the homogeneous and the time-varying parts, due to the aggregation property of independent Poisson processes. Under the coupling scheme described above, it is easy to see that the coupled M/M/1 queue dominates the time-varying one in Q t (the number of customers in the system), because the arrivals in the time-varying queue are a subset, and the potential departures a superset, of those in the coupled M/M/1 queue.
Conceptually, we start the dominating homogeneous queue and the coupled timevarying queue infinitely long ago. At time 0, both of them are in steady state. In a past time τ ∈ R − , if Q H τ = 0, then the coupled time-varying queue must also be empty at this time. By running it forward with the time-varying events generated as above, we get a steady-state draw of the time-varying queue at time 0.
In practice, only a finite number of values are needed. The algorithm is described as follows: 1. We simulate the M/M/1 dominating queue's stationary value at time 0. In station- Filter the arrival events generated by step 2 above by thinning as follows: an arrival event at time ζ is retained with probability λ(ζ )/λ u , otherwise it is deleted from the collection of arrival events. 4. Supplement the potential departures according to a Poisson process with rate μ(t)− μ l . Based upon our assumptions in Sect. 2, since 1 0 μ(t)dt =μ < ∞, we can proceed as follows: Let t 0 = τ , and repeat the two sub-steps below for where the subscript t n indicates that F t n (x) depends on t n . (2) Assign t n+1 = t n + X t n . If N > 0, then append {t 1 , . . . , t N } to the collection of potential departure events, and sort these in ascending order, yielding the aggregate collection of potential departure events for the time-varying queue. 5. Starting from the empty state at time τ − with the set of arrival and potential departure instants generated in the previous two steps, run the time-varying system forward until time 0 and output the state Q N 0 as a steady-state draw at an integral time for the M t /M t /1 queue.
Append t to Arrivals 10: else 11: Append t to In this section, we use the regenerative method to perform perfect sampling of the M t /M t /1 queue with the general stationary condition, i.e.,λ <μ.
We start with a general description of the regenerative method to obtain a stationary draw from our distribution of interest, which we adapt from [19] in what follows. Let X n , n = 0, 1, . . . denote the number of customers in a stable queue just before the (n + 1) st arrival, with X 0 = 0. Then {X n } n≥0 is a positive recurrent non-delayed discrete-time regenerative process with X n = 0 as the regenerative setting. Assume its cycle length is T ∈ N, with E(T ) < ∞, where the cycle length is defined as Explicitly, a generic cycle with length T can be defined as It is easy to simulate i.i.d. cycles and the sequentially generated ones are denoted by with corresponding cycle lengths T ( j) . Denote by T e a random variable which has the equilibrium distribution of the cycle length. Suppose we can sample T e , and let J = min{ j ≥ 1 : T ( j) ≥ T e }, then we have a steady-state draw of {X n } n≥0 as (Note that, if T (J ) = T e , then X = 0.) The interested reader will find a proof of the correctness of the regenerative method to generate a stationary draw in [19]. More details on regenerative methods can be found in [2, p. 420] and [3].
The key to our algorithm is to construct a dominating process which can be simulated in steady-state. If we start with a realization of the time-varying system, and concentrate all arrivals to the end of the interval (k − 1, k], k ∈ N, and all potential departures to the beginning of it, the modified process would dominate the original one in the sense that at whole integer times it will have at least as many customers waiting. Intuitively, since more potential departure events might be lost to an empty queue due to the postponing of the arrival events, there would tend to be more customers remaining in the system after the arrivals than at the same instant in the original process. This idea is explicitly stated by the following proposition.
Proof At the non-integer points, we define It is obvious that since no matter what the system's state is, the N A k arrivals provide a lower bound. It is clear that when k = k 0 the inequality is true. Assume that when k = m, m ≥ k 0 , m ∈ Z, the inequality holds, then when k = m + 1 we have one of the following situations: , the time-varying queue keeps busy on this interval), so that (1) If L t > 0, ∀t ∈ (m, m + 1), then and it follows that (2) Otherwise ∃t ∈ (m, m + 1) such that L t = 0, then it must be the case that Thus in both cases, the inductive step is established, and the result follows.
It is clear that {L k } k≥0 is a non-delayed regenerative process with L k = 0 as the regenerative setting. So its cycle length can be defined as Then we have

5)
, and N A k is independent of these r.v.'s as defined in Sect. 2. The limiting random variable L ∞ is defined by lim k→∞ L k . A segment of a sample path of the dominating process is shown in Fig. 1. It resembles the Late Arrival System discrete queue of [6], but the differences preclude us from using the LAS model directly. Instead, we exploit the form of (3.5) directly.
The differences Z k constitute an i.i.d. sequence, which we generically denote by Z = N A − N D . In light of (3.5) we find So the limiting random variable L ∞ , defined by lim k→∞ L k , satisfies To perform the Exponential Change of Measure (ECM) [2, p. 129], solve for γ > 0, where The changed measure is given by where P stands for the original measure. Equation (3.7) has the equivalent form λe γ +μe −γ −λ −μ = 0.

Consequently
due to the convexity of g(θ ) as shown above. This implies that {S k } k≥0 becomes a random walk with positive drift under the measure of P γ . So Z (under the measure P γ ) can be simulated by generating N A * and N D * from their respective distributions, and then taking the difference. Now using P γ , define a strictly increasing process with ladder heights S τ (n) , n = 0, 1, . . ., where Then L is a stationary draw of L ∞ . Thus the stationary draw of L ∞ given by where N A ∼ Poi(λ), and N A is independent of L .

Algorithm for perfect sampling of the M t /M t /1 queue
Based on the constructed dominating process, whose stationary state can be simulated, the perfect sampling of the M t /M t /1 queue is available using the regenerative method.
1. Simulate a random variable (denoted by T e ) from the equilibrium distribution of the cycle length (3.2) of the regenerative process of {L k } k≥0 , which dominates the queue length process in the time-varying queue at the integer time points. We obtain T e as follows. At time 0, sample a stationary draw of the dominating process, denoted by L 0 , using the method presented in the previous subsection. Continue simulating this process forward until it becomes 0. According to equations (3.5) and (3.6), 3. Use the order statistics method of simulating the time-varying Poisson process (see [17, p. 700, Method 2]) to construct time-varying events (arrival and potential departure instants) according to N A k and N D k (1 ≤ k ≤ T (J ) ) generated in cycle C (J ) . The coupling scheme specified in Proposition 1 implies that the time-varying instant can be generated by shifting the coupled homogeneous one on each interval (k − 1, k]. Let t N and t H be the instants in the time-varying and homogeneous systems, respectively, on this interval. For each such instant of a homogeneous arrival or service event on the interval, we know that independent of all other events, it will be uniformly distributed on (k − 1, k]. The corresponding timevarying instant can be computed as (3.8) where F −1 corresponds to the inverse of functions F λ (t) or F μ (t) defined in equation (2.1). From time 0, where the system is empty, simulate forward with these inputs to restore the time-varying queue. Output Q N T e as the stationary draw of the number of customers in the M t /M t /1 queue at integer time points.

Remark 1 (1) According to the regenerative method, L (J )
T e is a steady-state draw from the dominating process. Since it is coupled with the time-varying queue, at this time point (T e ), the corresponding sample of Q N T e is also stationary. (2) At time 0, even if the stationary draw of L 0 equals zero, we still continue simulating forward. (3) Since {L k } k≥0 is the dominating process, we can also directly output Q N 0 = 0 if L 0 = 0. But in this case, the condition of stopping the generic cycle simulation becomes (4) Unfortunately the regenerative method has infinite expected runtime [5,20]. Thus some runs may take so long that a practitioner would abort them, introducing a bias into the simulation.

An example
Let These parameters are the same as those used by [14], and haveλ = 1 andμ = 4. The regenerative method described in the previous subsection has been applied for the simulation. On a unit cycle, we chose 100 points at equal spacing from 0 to 1 and generated 10,000 samples for each point. Since only Q N 0 is generated in each trial, to get the samples at different points we changed the phases of the sinusoid functions in each run. For every point we repeated the algorithm 10,000 times. Although it is time consuming (around 100 times more than simulating the successive 99 points by continuing running the system after sampling the first point with the regenerative method), the simulation shows that our method works quite well. In this example, we did not need to abort any long runs. Since the occupancy (0.25) is quite small, they are very rare.
The idle probability and expected number of customers at time t ∈ (0, 1) of the time-varying queue are illustrated in Fig. 2. The gray areas indicate pointwise 95 % confidence intervals. The time average of Q N t is around 0.383. The simulated values match very well with the analytical results derived by [14,Section 3], which involve solving a Volterra equation of the second kind numerically.
For a more efficient simulation, we could repeat the algorithm 10,000 times for just one phase, and continue simulating the process through a whole cycle using standard forward simulation methods. In this case, samples at different time points in the simulation would be correlated and could be studied as draws from their joint distribution. In this section, we present a perfect simulation algorithm with a service time distribution that does not vary with time. Service durations (denoted by B) are drawn from some general distribution G(·). We require that E(B 2 ) < ∞ in order to ensure that the algorithm for the backward simulation of the coupled M/G/1 queue has finite expected runtime. Since the service requirement of a customer can be considered to be known at an arrival instant, we can take the perspective of analyzing the unfinished workload to explore this time-varying system. As a result, this case is easier to handle than the M t /M t /1 queue. As was the case before, the first step is to find a dominating process. In the next subsection we construct a process which dominates the time-varying queue in the unfinished workload. Using dominated CFTP [11] we achieve perfect sampling of the M t /G/1 queue.

The dominating process
Proof We use mathematical induction to prove this proposition. Clearly V N k ≤ V H k +1 for k = t 0 , since both are 0. For larger k, let ξ k be the additional workload that arrives during the interval (k − 1, k] (the same for both queues). Let η N k be the amount of work done on these new customers during the interval (k − 1, k] in the time-varying queue, and let η H k be the counterpart in the homogeneous queue. On the interval (k −1, k], unless the server finishes the remaining workload (carried from previous intervals) within that cycle, it cannot direct any capacity to serve the new arrivals on this interval. Therefore we obtain Note that η N k ∈ [0, 1), and that η N k = 0 if V N k−1 > 1, as there remains earlier work in the system to be done.
(by the inductive hypothesis) ≤ 1 which again agrees with our result.

Backward simulation of the M/G/1 queue
The previous subsection established that the upper bound of the dominating process can be estimated with the unfinished workload of the coupled homogeneous M/G/1 queue. To perform the dominated CFTP algorithm, we simulate the M/G/1 queue backwards using the time reversibility of its Processor-sharing (PS) variant. (Under the PS discipline, the customers share the server, i.e., when n customers are present, the server devotes 1/n of its capacity to each. The customers attain service at rate 1/n. The first such customer leaves the system once the attained service reaches the minimum of the residual service times, unless another customer arrives first in which case the rates are reduced to 1/(n + 1) [1, p. 63].) Ross [16, p. 280] showed that the M/G/1 PS model is time reversible and at stationarity the number of customers (Q) in the system has the geometric distribution where ρ is the occupancy; the completed (or unfinished) workload of the M/G/1 PS queue is where the Y i 's are i.i.d. and they follow the equilibrium distribution of the service duration.
Since the sample paths of the unfinished workload for the M/G/1 FCFS queue are identical to those of the M/G/1 PS model, we can use the latter to achieve the backward simulation of the M/G/1 FCFS queue. Sigman [18,Algorithm1.1,Step1] described this algorithm, which we restate as follows for the sake of completeness: In Algorithm 2, the vector Instants is used to store the departure instants (when running the PS model forward) and arrival instants (after being reversed), and Services the associated service requirements. They are initially empty. When looking backwards, ζ is the stationary age of the busy period of the M/G/1 FCFS queue.
It is well known (for example, [12, pp. 213-214]) that so we need the existence of E(B 2 ) in order for this algorithm to work in finite time.

Algorithm for perfect sampling of the M t /G/1 queue
The perfect sampling of the M t /G/1 queue is performed using dominated CFTP. Imagine that an M t /G/1 queue and the coupled dominating process {V H k + 1} k∈Z , 4. Starting from −ζ 1 with unfinished workload V N −ζ 1 = V H −ζ 1 + 1, run the timevarying queue forward with the collection of time-varying arrival instants and the service durations from the previous steps, and output V N 0 as a stationary draw of the unfinished workload at time 0 of the M t /G/1 queue. Proposition 3 By following the algorithm for simulating the M t /G/1 queue specified above, the output of V N 0 is a stationary draw of the unfinished workload at time 0. Proof Assume an M t /G/1 queue and an M/G/1 queue were started infinitely long ago and coupled in the way described in Proposition 2. So V N k ≤ V H k + 1, ∀k ∈ Z. The algorithm as constructed above ensures that there exists idle time on interval ( −ζ 1 , −ζ ) in the coupled time-varying queue. When mapping the homogeneous arrivals to time-varying ones, the maximum work originally finished after −ζ that could be rearranged prior to −ζ is V H −ζ . However, since V H t becomes zero on the interval [ −ζ , −ζ + 1), it follows that V H −ζ < 1. According to the dominance, it follows that V N −ζ 1 ≤ V H −ζ 1 + 1. Thus based on the homogeneous system the extra work introduced from the customers arriving prior to −ζ 1 is no greater than 1. Consequently, the 2 units of idle time on interval ( −ζ 1 , −ζ ) are guaranteed by the algorithm to be sufficient to ensure that there is idle time for the time-varying queue on this interval.
Finally, this idle time existence ensures the coalescence of all possible chains of the time-varying queue. At time −ζ 1 , all possible chains of the M t /G/1 queue have unfinished workload ranging from 0 to V H −ζ 1 + 1. They are driven by the same arrival instants and associated homogeneous service durations, and it is easy to see that this coupling is monotone. So when the upper bound chain (corresponding to V N −ζ 1 = V H −ζ 1 +1) becomes zero, so do all the coupled chains, ensuring coalescence. The usual CFTP argument then establishes our result.

Examples
Here we illustrate sampling the stationary unfinished workload in the M t /G/1 queue with Pareto and Erlang distributions of the service durations. In both cases, the arrival rates are λ(t) = 3 + 3 sin(2π t).
The average unfinished workload and 95 % confidence intervals are illustrated in Fig. 3. It is evident from the figure that they behave in periodic patterns with the same periodic lengths as the arrival processes.

Conclusions
In this paper, we have presented algorithms for perfect sampling of the M t /M t /1 and M t /G/1 queues, where the time-varying ingredients have periodic patterns. The stationary queue length and virtual waiting time (unfinished workload) have timedependent distributions.
For the M t /M t /1 queue with relatively light load, i.e. sup λ(s) < inf μ(t), we implemented perfect sampling with the dominated CFTP by constructing a homogeneous M/M/1 queue with arrival and service rates of sup λ(s) and inf μ(t), respectively. In the general case whereλ <μ, we constructed a discrete dominating process, and applied the regenerative method to get the time-dependent steady-state draws of the time-varying queue.
Perfect sampling of the M t /G/1 queue was achieved using dominated CFTP. The dominating process was computed based on the unfinished workload of the coupled homogeneous queue (M/G/1), which can be simulated backwards using its Processorsharing variant. Unlike most CFTP implementations, our algorithm allows direct calculation of the (random) starting time in the past, rather than the usual trial and error approach. This is in common with the dominated CFTP algorithm for the M t /M t /1 queue, and also with the dominated CFTP algorithm by [18]. As with other CFTP algorithms, it yields a draw from the stationary distribution at time 0. Other time points in the cycle are simulated simply by adjusting the origin of the clock.