DP2014-37

Perfect Simu lation for Models ofIndustry Dynamics
Takashi KAMIHIGASHIJohn STACHURSKI

November 21, 2014

Perfect Simulation for Models of Industry
Dynamics∗
Takashi Kamihigashia,b and John Stachurskic
November 19, 2014a Research Institute for Economics and Business Administration, Kobe University
b IPAG Business School, Paris, Francec Research School of Economics, Australian National University

Abstract

In this paper we introduce a technique for perfect simulation from the stationary distribution of a standard model of industry dynamics. The method can beadapted to other, possibly non-monotone, regenerative processes found in industrial organization and other fields of economics. The algorithm we propose is aversion of coupling from the past. It is straightforward to implement and exploits
the regenerative property of the process in order to achieve rapid coupling.
Keywords: Regeneration, simulation, coupling from the past, perfect sampling

1

Introduction

The dynamics of entry and exit by firms play an essential role in economic theory as
well as in real life. Productive new entrants replace unproductive incumbents, rejuvenating overall economic activity. There is a large and growing literature on this economic mechanism (see, e.g., Hopenhayn and Rogerson (1993), Cooley and Quadrini∗ We are grateful for financial support from ARC Discovery Award DP120100321 and the Japan Soci-ety for the Promotion of Science. Email: tkamihig@rieb.kobe-u.ac.jp, john.stachurski@anu.edu.au.

1

(2001) or Melitz (2003)), and much of this literature builds upon the model of entry andexit studied by Hopenhayn (1992). The stationary distribution of entry-exit models of
the type studied by Hopenhayn represents a cross-sectional distribution of firms thatis both consistent with the definition of equilibrium at any point in time and also invariant over time. For typical parameter values the stationary distribution is uniquely
defined but has no analytical solution.1
Simulation is a useful option for computing the cross-sectional distribution for a given
exit policy, since it is not difficult to write down an ergodic Markov process such that its
stationary distribution coincides with the cross-sectional stationary distribution. Thispermits approximate sampling from the stationary distribution simply by running theprocess from an arbitrary initial value until it is judged to have “nearly” converged.
Simulating until the distribution of the state is approximately stationary is referred
to as “burn-in.” Unfortunately the length of burn-in required is often the subject ofguesswork and heuristics.2 Moreover, regardless of how much burn-in is performed,
the resulting sample is never exactly stationary, and the size of the error is once again
unknown.3
In this paper we show that these problems can be overcome. By using a variation of
coupling from the past (CFTP) technique originally due to Propp and Wilson (1996),
we show that it is possible to perform perfect sampling—that is, to sample exactly fromthe stationary distribution of this class of models—for any specified exit threshold. In
particular, we develop an efficient algorithm that generates exact, IID draws from thestationary distribution. For each random seed, the algorithm terminates as soon as anexact draw has been generated, and it is guaranteed to terminate in finite time withprobability one. Hence there is no need for the heuristics used to judge burn-in time.
Moreover, by repeating the algorithm with independent seeds it becomes possible togenerate multiple independent draws from the stationary distribution.
Our work draws on a large literature on CFTP that mainly exists outside of economics,
1 The difficulty of obtaining analytical solutions for this model is related to the existence of a positive

threshold for productivity at which firms choose to exit. The threshold introduces a nonlinearity thatessentially rules out analytical methods.
2 While some methods for computing error bounds exist, they are rarely used for two reasons. First,
they are nontrivial to compute. Second, these bounds are often highly pessimistic, since any such bounds
must address the worst case scenario admitted by the model specification.
3 A related issue is that, for a given method, the size of the error is likely to vary with the parameters,
since the parameters change the structure of the problem. If the burn-in is not varied accordingly, this is
likely to cause bias.

2where the technique is often used for models with large but discrete state spaces.4
More recently, researchers have developed techniques for implementing CFTP methods in continuous state settings. Murdoch and Green (1998) showed that CFTP canin principle be used in continuous state settings when the underlying Markov processsatisfies Doeblin’s condition. This condition requires the existence of a nonnegativelower bound function that (a) integrates to a positive value, (b) depends only on thenext state, and (c) is pointwise dominated by the transition density function (which
depends on both the current state and the next). Theoretical work along the same lines
can be found in Foss and Tweedie (1998) and Athreya and Stenflo (2003).
Although these results are fundamental, they can be difficult to apply. Murdoch andGreen (1998) admit that their basic method, which is in principle applicable to ourmodel, may have “a limited range of application for two reasons.” First, the function
associated with Doeblin’s condition “may be too small for practical use” to generate
exact draws in a reasonable length of time. Second, their method requires the user
to draw from scalar multiples of the lower bound transition density and a residual
kernel. It can be nontrivial or even impossible to explicitly calculate and draw fromthese distributions. If approximations are required, this to some degree defeats the
purpose of CFTP.
For these reasons, CFTP methods tend to be popular only in specific settings, perhaps
the most notable of which is where the underlying Markov process is stochasticallymonotone. For such processes, efficient and straightforward CFTP methods are available. Corcoran and Tweedie (2001) developed general results on CFTP particularly
suitable for monotone Markov processes. An application to economics can be found in
Nishimura and Stachurski (2010), where monotonicity makes the algorithm straightforward to implement and analyze.
Here we develop a CFTP algorithm that is designed to produce exact draws from thestationary distribution of Hopenhayn’s entry-exit model, which is not monotone. We
do however exploit some monotonicity properties from Hopenhayn’s model in ouralgorithm. We show that the algorithm terminates successfully in finite time withprobability one by using both the monotonicity of productivity for incumbents andthe regenerative property introduced by new entrants. Our algorithm is distinct from
Murdoch and Green’s method discussed above (Murdoch and Green (1998)), in thatit does not use Doeblin’s condition, and does not require explicit knowledge of the
4 Applications range from statistical mechanics to page ranking and the design of peer-to-peer filesharing systems. See, for example, Propp and Wilson (1996), Kijima and Matsui (2006), Huber (2003)
and Levin et al. (2009).

3

transition density.5 As long as one can simulate the overall Markov process, one cansample exactly from the stationary distribution using the algorithm.6
Aside from related models of industry dynamics, our techiques can also potentially beapplied to other non-monotone regenerative models, such as those found in various
intertemporal decision problems. One example is the problem of optimal replacementof a part or machine, the performance of which degrades stochastically over time (see,e.g., Rust (1987)).

2Preliminaries

2.1

The Entry-Exit ModelIn this section we briefly review a benchmark model of firm dynamics due to Hopenhayn (1992). The model is set in a competitive industry where entry and exit is endogenously determined. In the model there is a large number of firms that produce
a homogeneous good. The firms face idiosyncratic productivity shocks that follow a
Markov process on S := [0, 1]. The conditional cumulative distribution function forthe shock process is denoted by F (φ0 | φ). Following Hopenhayn (1992), we impose thefollowing restrictions:
Assumption 2.1. F is decreasing in its second argument and, for any e > 0 and any
φ ∈ S, there exists an integer n such that F n (e | φ) > 0.7
We let P denote the stochastic kernel on [0, 1] corresponding to F. That is, P(φ, A) :=R
1 A (φ0 ) F(dφ0 | φ) for φ ∈ S and A ∈ B, where B represents the Borel sets on [0, 1]
and 1 A is the indicator function of A. Incumbent firms exit the industry whenevertheir current productivity falls below a reservation value xt . Letting Mt be the mass ofentrants at time t and ν be the Borel probability measure from which the productivityof entrants is drawn, the sequence of firm distributions {µt } on S satisfies µt+1 ( A) =R
P(φ, A)1{φ ≥ xt }µt (dφ) + Mt+1 ν( A) for all A ∈ B. At the stationary equilibrium,
both x and M are constant, and a stationary distribution µ is a Borel probability8 mea5 The assumptions used to show the probability one termination of the algorithm in fact imply Doe-blin’s condition for some n-step transition, but our proof of this property does not use the latter.
6 Of course in computer implementations exactness is modulo the errors associated with floating
point arithmetic and imperfect random number generators.
7 F n (· | φ ) is the conditional distribution for productivity after n periods, given current productivity
φ.
8 We focus only on normalized measures, since other cases are just scalar multiples.

4

sure µ satisfyingZµ( A) =P(φ, A)1{φ ≥ x }µ(dφ) + Mν( A)( A ∈ B ).

(1)

It follows from (1) and µ(S) = P(φ, S) = ν(S) = 1 that M = M( x, µ) := µ{φ ∈ S :φ < x }. As a result, we can also write (1) asµ( A) =
whereZ

Q(φ, A)µ(dφ)

Q(φ, A) := P(φ, A)1{φ ≥ x } + ν( A)1{φ < x }.

(2)

(3)

Equation (2) states that µ is a stationary distribution for the stochastic kernel Q in theusual sense of time invariance. As shown by Hopenhayn (1992), the kernel Q has only
one stationary distribution. For the purposes of this paper we will treat x as given. For
typical parameter values the stationary distribution has no analytical solution.

2.2

Simulation

It is not difficult to produce an ergodic Markov process suitable for simulation suchthat its stationary distribution (i.e., time-invariant distribution) coincides with the crosssectional distribution µ in (2). In essence, we need a method for sampling from thestochastic kernel Q. The first step is to simulate from the conditional distribution
P(φ, ·) = F (· | φ). In particular, we seek a random variable U and a function g suchthat D( g(φ, U )) = F (· | φ) for all φ ∈ S. (Here D( X ) indicates the distribution of random variable X.) This can be achieved via the inverse transform method, where U is
uniform on [0, 1] and g(φ, u) = F −1 (u | φ).9 Now consider the process {Φt } defined byΦt+1 = g(Φt , Ut+1 )1{Φt ≥ x } + Zt+1 1{Φt < x }(4)where {(Ut , Zt )} is IID with D( Zt ) = ν and D(Ut ) = Uniform[0, 1]. In what follows
we call (4) the simulation model.
Lemma 2.1. The simulation model is a Markov process with stochastic kernel Q.
Proof. Fix A ∈ B. Applying the law of iterated expectations to (4), we haveP{Φt+1 ∈ A} = E [Et 1{Φt+1 ∈ A}]= E [Et 1{Φt+1 ∈ A}1{Φt ≥ x } + Et 1{Φt+1 ∈ A}1{Φt < x }]= E [ P(Φt , A)1{Φt ≥ x } + ν( A)1{Φt < x }].
9 Here F −1 (· | φ ) is the generalized inverse of F (·| φ ). That is, F −1 ( u | φ ) : = inf{ z : F ( z | φ ) ≥ u }.

5

Hence if µt := D(Φt ) for all t, then µt+1 ( A) = E Q(Φt , A) =words, {Φt } is Markov with stochastic kernel Q.

R

Q(φ, A)µt (dφ). In other

Under some additional assumptions one can also show that the simulation model isuniformly ergodic. Hence it can be used to estimate µ via ordinary forward simulation,
starting from an arbitrary initial condition. As discussed in the introduction, we adopt
a different approach, which permits exact sampling.

3

Perfect Simulation

Before turning to our algorithm it can be helpful to recall some features of CFTP in a
simpler setting. To this end we first review a version for finite state Markov chains andthen move on to industry dynamics.

3.1

Coupling From the Past

To begin, consider an irreducible aperiodic Markov chain { Xt } on finite state space S
with unique stationary distribution µ. Suppose that the dynamics of the model are
given by Xt+1 = h( Xt , Wt+1 ) where {Wt } is IID. If we start at any X0 and simulate{ Xt } by drawing W1 , W2 , . . . and successively applying the maps h(·, W1 ), h(·, W2 ), . . .,then the distribution of Xt will converge to (but in general never reach) the stationary
distribution µ. A second immediate observation is that if { Xt } and { Xt0 } are two chains
simulated using the same sequence of shocks {Wt }, and if Xτ = Xτ0 for some τ, thenXt = Xt0 for all t ≥ τ. We say that { Xt } and { Xt0 } couple at time τ.
Now consider an IID sequence of shocks {Wt }0t=−∞ indexed on {. . . , −2, −1, 0}, andlet { Xt }0t=−∞ be the process generated from these shocks, starting in the infinite pastand terminating at t = 0. Intuitively, since the burn-in is infinite, the distribution of X0
will be exactly µ. On the other hand, because the chain { Xt }0t=−∞ is infinite and hencecannot be simulated, it appears that X0 is not observable. In fact this is not the case.
To understand how X0 can be observed, suppose that we fix T ∈ N, take the finite
subset of shocks W−T +1 , W−T +2 , . . . , W0 and then, for each point x in the state space S,
construct one version of the process { Xt ( x )}0t=−T that starts at this point x at time − T
and runs up until time zero. All of these processes { Xt ( x )}0t=−T are updated using thesame shocks W−T +1 , WT +2 , . . . , W0 . We will call them tracking processes. Because there
are only finitely many tracking processes and each is of finite length, it is possible to
6

10
8
6
4
2
0
10
8
6
4
2
0

−T

0

−T

0

Figure 1: Coalescence of the tracking processes
simulate them. Moreover, following the tracking processes can allow us to pin down
the value of X0 .
To understand the idea, consider first the two different visualizations of the tracking
processes shown in the top panel of figures 1 and 2 respectively. Here the state spaceis {1, . . . , 10}, and the two figures show different realizations from the same model.
Notice that in the simulations, some of the processes couple and then run together up
until time zero. In figure 1, all of the tracking processes have coupled by time zero,and only one final value is observed. If this is the case we will say that the tracking
processes coalesce. In figure 2 multiple paths are still present at time zero, and theprocesses fail to coalesce.
Now consider again the process { Xt }0t=−∞ generated from the entire sequence of shocks{Wt }0t=−∞ . Since this process must pass through one point x in the state space at − T,and since it receives the same shocks W−T +1 , W−T +2 , . . . , W0 as the tracking processes
from that time forwards, over the period − T, − T + 1, . . . , 0 it must follow the same
path as the tracking process { Xt ( x )}0t=−T that started at x. Of course we do not know
which of the tracking processes it is following. However, if the tracking processes coalesce, then there is only one final observation. This observation must be a realization of
the time zero value X0 of the process { Xt }0t=−∞ that started in the infinite past, and istherefore a draw from the stationary distribution µ. Such an outcome is illustrated in
the bottom panel of figure 1. For comparison, an unsuccessful outcome is illustrated
in the bottom panel of figure 2. Here there are three final values, and we do not know
which is the time zero value of the process { Xt }0t=−∞ . In this case we can take one step
7

10
8
6
4
2
0
10
8
6
4
2
0

−T

0

−T

0

Figure 2: The tracking processes failing to coalesce
further back in time, drawing the additional shock W−T while preserving the existing
shocks W−T +1 , W−T +2 , . . . , W0 , recalculate the tracking processes, test for coalescence,and so on. This procedure will eventually terminate with an exact draw from µ.10

3.2

CFTP with the Simulation Model

While the technique presented in the previous section works when there are a finite
number of states, it is clearly not possible to directly implement the same idea when
the state space is infinite. However, variations on the idea can potentially be foundby exploiting the structure of a given model. In particular, for the simulation model
defined in (4), we now show how the monotonicity of g and the renewal component
can be combined to pin down the final value X0 (which is Φ0 in this case).
To begin, let {(Ut , Zt )}t∈Z be an infinite sequence of IID shocks indexed on Z and witheach pair (Ut , Zt ) having the product distribution Uniform[0, 1] × ν. To simplify notation we will let gt := g(·, Ut ), so that, for example, gt · · · g1 φ := gt ◦ gt−1 ◦ · · · ◦ g1 (φ)
is exogenous productivity at t, given time zero productivity φ ∈ S. To further simplify
notation, let
ht (φ) := g(φ, Ut )1{φ ≥ x } + Zt 1{φ < x },so that (4) becomes Φt+1 = ht+1 Φt := ht+1 (Φt ).
10 Using finiteness of the state space, aperiodicity and irreducibility, it can be shown that termination

occurs in finite time with probability one.

8

Now fix T ≥ 1. For each φ ∈ S, there is a corresponding “tracking process” that starts
at time − T with value φ, and then updates with maps h−T +1 , h−T +2 , . . . , h0 , obtaining
the value h0 · · · h−T +1 φ at time zero. We say that the tracking processes coalesce if, forsome T ∈ N, the set of final values
h 0 · · · h − T +1 ( S ) : = { h 0 · · · h − T +1 φ : φ ∈ S }(5)

is a singleton. What we will now show is that under mild conditions coalescence occurs with probability one, and, moreover, that it is not necessary to keep track of thefull continuum of tracking processes in order to find the value of the singleton. In particular, we show that, conditional on a certain event described below, the set of final
values h0 · · · h−T +1 (S) has only finitely many possibilities. Hence coalescence occurs
whenever these finite possibilities take the same value. All of these finite possibilities
are computable. To begin describing them, let T > 1 be given, letΣ T := {k ∈ N : 1 ≤ k < T and g−T +k · · · g−T +2 · g−T +1 1 < x },and let σT := min Σ T . Intuitively, σT is the number of periods that an incumbent firm
survives, given that it starts at time − T with maximal productivity 1 and faces theshock sequence {Ut }t∈Z . Clearly σT is only defined when Σ T is nonempty. However,the probability that Σ T is nonempty converges to one as T → ∞ by assumption 2.1.
Moreover, it is remarkable that if Σ T is nonempty, then the set h0 · · · h−T +1 (S), whichcontains the final values of the tracking processes started at − T, can have only finitely
many values:
Lemma 3.1. If Σ T is nonempty, then h0 · · · h−T +1 (S) ⊂ Λ T , where the latter is the finite set
given by
Λ T := { h0 · · · h−T +k+1 Z−T +k : k = 1, . . . , σT + 1}.
(6)
The proof of lemma 3.1 is given in section 4. The intuition for the result can be obtainedby considering figure 3. In the figure, T = 10. Tracking processes are plotted for 50
different initial values of φ ∈ S. (Ideally, tracking processes would be plotted fromevery φ ∈ S, but this is clearly impossible.) For this particular realization of shocks,
the set Σ T is nonempty because the process starting from 1 at time −10 falls below x
at t = −2 (and hence σT = 10 − 2 = 8). As is clear from the figure, the fact that theprocess starting from 1 at time −10 falls below x at t = −2 implies that all tracking
processes fall below x at least once between −10 and −2 (recall that the productivity ofincumbents is monotone). Moreover, if any collection of tracking processes fall below
x at some point in time t, they subsequently couple, taking the common value Zt+1 at
9

1

x

0 10

9

8

7

6

5

4

3

2

1

0

Figure 3: Tracking processes with T = 10 and σT = 8
t + 1 and being equal from then on. As a result, by t = −1, there are at most σT + 1 = 9
distinct tracking processes. Their time zero values are included in the set Λ T definedin lemma 3.1. In particular, Λ T is the time zero values of the processes that start below
x at dates −10, −9, . . . , −2.
To see the importance of lemma 3.1, let {Φt }t∈Z be a stationary, doubly-indexed process on the same probability space as {(Ut , Zt )}t∈Z that obeys Φt+1 = ht+1 Φt for allt ∈ Z. The common marginal distribution of Φt is µ. Since Φ−T lies somewhere in
S, we know that Φ0 = h0 · · · h−T +1 Φ−T ∈ h0 · · · h−T +1 (S). Moreover, if the set Σ T isnonempty, then lemma 3.1 yields the inclusion h0 · · · h−T +1 (S) ⊂ Λ T , and Φ0 lies in
the finite observable set Λ T . In particular, if Λ T is a singleton, then the value of Φ0 isrevealed as the value of that singleton.
Figures 4 and 5 show simulations with successful and unsuccessful coalescence respectively. In each figure, the top panel shows only the tracking processes. (As with
figure 3, the full continuum of tracking processes cannot be plotted, so we show only50.) The bottom panel shows the tracking processes and the path of {Φt }. In reality,
the path of {Φt } is not observed. However, in figure 4, there is only one final, coalesced value, and Φ0 must take this value. Hence Φ0 is observed. On the other hand,in figure 5, Φ0 is equal to one of two final values, and we have no way of identifying
which one it is.
Now let us consider how to use our results to sample from µ by generating observations of Φ0 . In order to avoid conditioning on coalescence by a certain point in time,
we wish to reveal the value of Φ0 for every random seed. This can be done by fixing
10

1

0 60

50

40

30

20

10

0

50

40

30

20

10

0

1

0 60

Figure 4: Successful coalescence from T = 50

1

0 60

50

40

30

20

10

0

50

40

30

20

10

0

1

0 60

Figure 5: Failure of coalescence11the seed, which determines the values of the shock processes, and then taking T larger
and larger until coalescence occurs. Algorithm 1 gives details.
Algorithm 1: Generates an exact draw from µ
fix T to be an integer greater than 1;
draw (U0 , Z0 ), . . . , (U−T +1 , Z−T +1 ) independently from their distributions;
repeatcompute the set Σ T ;if Σ T is nonempty thencompute the set Λ T ;if Λ T is a singleton thenset X to be the value of that singleton ;break ;
end
end
draw (U−T , Z−T ) and append to list (U0 , Z0 ), . . . , (U−T +1 , Z−T +1 ) ;set T = T + 1 ;end
return X ;We still need to show that (a) the distribution of X is µ, and (b) that the algorithm willterminate with probability one. We now show that these claims are true under thefollowing condition, which is satisfied by many standard distributions.
Assumption 3.1. If G ⊂ S is a nonempty open set, then ν( G ) > 0.
Let T ∗ := min{ T ∈ N : Σ T is nonempty and Λ T is a singleton }. In other words, T ∗ isthe length of the simulation run required for algorithm 1 to terminate. Evidently T ∗ isa random variable.
Proposition 3.1. If assumption 3.1 holds, then
1. There exists a γ ∈ (0, 1) such that P{ T ∗ > t} ≤ γt .
2. The random variable X returned by algorithm 1 has distribution µ.
Two comments are in order. First, item 1 tells us in particular that P{ T ∗ < ∞} = 1. Inother words, the algorithm terminates in finite time with probability one. Thus proposition 3.1 not only gives probability one occurrence but also provides a geometric rate.
12

3.0
2.5
2.0
1.5
1.0
0.5
0.00.0

0.2

0.4

0.6

0.8

1.0

Figure 6: Stationary density
Second, in view of item 2, algorithm 1 gives us a way to generate an exact sample fromµ. If we then repeat the process with independent shocks we will generate multiple
independent draws. The proof of proposition 3.1 is given in section 4.11
The web page https://github.com/jstac/hh sampling contains a simple implementation of algorithm 1. We tested the code by following Hopenhayn and Rogerson (1993)
in taking the distribution ν for new entrants to be uniform, and the process for incumbents to be Φt+1 = gt+1 Φt = a + ρΦt + et+1 where {et } is IID with distribution
N (0, σ2 ). To bound the process we added reflecting barriers at 0 and 1. The parameters
were set to a = 0.36, ρ = 0.4 and σ = 0.1, while x was set to 0.49, so that approximately
40% of incumbents exit within 5 years (Hopenhayn, 1992, p. 1127). For these parameters, running the program on a standard workstation without parallelization produces
about 36,000 independent draws from µ per second.12
Figure 6 shows the density computed from 36,000 observations combined with a standard nonparametric kernel density estimator (using a Gaussian kernel). Figure 7 shows
a 95% confidence set for the cumulative distribution function corresponding to µ, basedon the same observations and calculated using the Kolmogorov distribution of the supnorm deviation between true and empirical cdfs. The Kolmogorov result is applicablehere because the draws are exact and IID. The true distribution function lies entirely
between the two bands with 95% probability.
11 A final point of intuition before the proofs is that on first reading it might appear that X is a draw

from the distribution ν of new entrants, rather than the stationary distribution µ, since coalescence to asingle value always occurs after a re-entry event. However, this is not the case because the algorithm
still runs on from this event until we reach time zero, at which point the final value is recorded. See
figure 4.
12 Our workstation has a 2.67GHz Intel CPU and 4 gigabytes of RAM.

13

1.0
0.8
0.6
0.4
0.2
0.00.0

0.2

0.4

0.6

0.8

1.0

Figure 7: 95% confidence set for the stationary distribution

4

ProofsIn the following proofs T is fixed, and we write σ for σT to simplify notation.
Proof of lemma 3.1. Let Σ T be nonempty. As a first step, we show that if φ ≥ x, thenthere exists a j ∈ {1, . . . , σ} such that h−T + j · · · h−T +1 φ < x. To see that this is so, fix
φ ≥ x and suppose that the statement fails. In other words, h−T + j · · · h−T +1 φ ≥ x forj ∈ {1, . . . , σ}. We know that if y ≥ x, then hi y = gi y. It follows that h−T +σ · · · h−T +1 φ =g−T +σ · · · g−T +1 φ. But thenx ≤ h−T +σ · · · h−T +1 φ = g−T +σ · · · g−T +1 φ ≤ g−T +σ · · · g−T +1 1 < x,where the second inequality is due to monotonicity of gi , and then third is by thedefinition of σ. Contradiction.
To complete the proof, pick any φ ∈ S. Our claim is that h0 · · · h−T +1 φ ∈ Λ T . Suppose
first that φ < x. In this case we have h0 · · · h−T +1 φ = h0 · · · h−T +2 Z−T +1 , which is anelement of Λ T . Next, suppose that φ ≥ x. In light of the preceding argument, there
exists a j ∈ {0, . . . , σ } with h−T + j · · · h−T +1 φ < x, and hence
h−T + j+1 · h−T + j · · · h−T +1 φ = Z−T + j+1 ,from which we obtainh 0 · · · h − T +1 φ = h 0 · · · h − T + j +2 · h − T + j +1 · h − T + j · · · h − T +1 φ

= h0 · · · h−T + j+2 Z−T + j+1 .
14

Since j ∈ {0, . . . , σ}, the right-hand side is an element of Λ T . This completes the proof.

Proof of proposition 3.1. We start with the first claim, which is existence of a γ ∈ (0, 1)
such that P{ T ∗ > t} ≤ γt .
Let n be an integer such that F n ( x | 1) > 0, existence of which is due to assumption 2.1.
Fixing j ∈ N, let
Ej := { g−( j−1)n−1 · · · g− jn 1 < x } ∩ { Z−( j−1)n−1 < x, . . . , Z− jn < x }.
The events { g−( j−1)n−1 · · · g− jn 1 < x } and { Z−( j−1)n−1 < x, . . . , Z− jn < x } are independent because the first event depends only on U−( j−1)n−1 , . . . , U− jn and the second
depends only on Z−( j−1)n−1 , . . . , Z− jn . As a result,δ := P( Ej ) = F n ( x | 1)ν([0, x ))n .
The constant δ is strictly positive as a result of assumption 3.1. We claim that if theevent Ej occurs, then Σ jn+1 is nonempty and Λ jn+1 is a singleton. To simplify notation,
we treat only the case of j = 1.
So suppose that E1 occurs. Clearly Σn+1 contains n, and hence is nonempty. To see thatΛn+1 is a singleton, observe that since σ = σn+1 is the smallest element of Σn+1 , wemust have σ ≤ n. As a consequence,Λn+1 = { h0 · · · h−n+k Z−n−1+k : k = 1, . . . , σ + 1}

⊂ {h0 · · · h−n+k Z−n−1+k : k = 1, . . . , n + 1}.
We claim that on the set E1 we haveh0 · · · h−n+k Z−n+k−1 = Z0

for anyk ∈ {1, . . . , n + 1}.

(7)

To prove that (7) holds, observe that on E1 the values Z−1 , . . . , Z−n are all less than x.
As a result, we haveh0 Z−1 = Z0h0 h−1 Z−2 = h0 Z−1 = Z0h0 h−1 h−2 Z−3 = h0 h−1 Z−2 = h0 Z−1 = Z0
and so on. Together, these equalities give (7). As a consequence, we conclude that Σn+1
is nonempty and Λn+1 is a singleton whenever E1 occurs, and, more generally, Σnj+1
15

is nonempty and Λnj+1 is a singleton whenever Ej occurs. The events E1 , E2 , . . . are
independent and have positive probability δ. Using the definition of T ∗ , we then haveP{ T ∗ > nj} = P{ T ∗ ≥ nj + 1} ≤ P ∪ij=1 Eic = (1 − δ) j
for all k ∈ N. Setting γ := (1 − δ)1/n gives the first result stated in proposition 3.1
The second claim in proposition 3.1 is that the value X returned by algorithm 1 has
distribution µ. It is clear from the algorithm that X is the value contained in Λ T when
the latter is a singleton. In the discussion preceding the proposition it was show thatwhen Λ T is a singleton the value of the singleton is Φ0 . Hence X = Φ0 . Moreover,{Φt } is by assumption a strictly stationary process, and hence the distribution of Φ0 isexactly µ.

References
Athreya, Krishna B and Örjan Stenflo (2003) “Perfect sampling for Doeblin chains,”
Sankhyā: The Indian Journal of Statistics, pp. 763–777.
Cooley, Thomas F and Vincenzo Quadrini (2001) “Financial markets and firm dynamics,” American Economic Review, pp. 1286–1310.
Corcoran, JN and RL Tweedie (2001) “Perfect sampling of ergodic Harris chains,” The
Annals of Applied Probability, Vol. 11, pp. 438–451.
Foss, Sergey G and RL Tweedie (1998) “Perfect simulation and backward coupling,”
Stochastic Models, Vol. 14, pp. 187–203.
Hopenhayn, Hugo A (1992) “Entry, exit, and firm dynamics in long run equilibrium,”
Econometrica, pp. 1127–1150.
Hopenhayn, Hugo and Richard Rogerson (1993) “Job turnover and policy evaluation:
A general equilibrium analysis,” Journal of Political Economy, pp. 915–938.
Huber, Mark (2003) “A bounding chain for Swendsen-Wang,” Random Structures &Algorithms, Vol. 22, pp. 43–59.
Kijima, Shuji and Tomomi Matsui (2006) “Polynomial time perfect sampling algorithm
for two-rowed contingency tables,” Random Structures & Algorithms, Vol. 29, pp. 243–
256.
16

Levin, David A, Yuval Peres, and Elizabeth L Wilmer (2009) “Markov chains and mixing times. With a chapter by James G. Propp and David B. Wilson,” American Mathematical Society, Providence, RI.
Melitz, Marc J (2003) “The impact of trade on intra-industry reallocations and aggregate industry productivity,” Econometrica, Vol. 71, pp. 1695–1725.
Murdoch, Duncan J and Peter J Green (1998) “Exact sampling from a continuous state
space,” Scandinavian Journal of Statistics, Vol. 25, pp. 483–502.
Nishimura, Kazuo and John Stachurski (2010) “Perfect simulation of stationary equilibria,” Journal of Economic Dynamics and Control, Vol. 34, pp. 577–584.
Propp, James Gary and David Bruce Wilson (1996) “Exact sampling with coupled
Markov chains and applications to statistical mechanics,” Random Structures and Algorithms, Vol. 9, pp. 223–252.
Rust, John (1987) “Optimal replacement of GMC bus engines: An empirical model ofHarold Zurcher,” Econometrica, pp. 999–1033.

17