DP2014-09

Perfect Simu lation for Models
of Industry Dynamics
Takashi KAMIHIGASHIJohn STACHURSKI
March 13, 2014

Perfect Simulation for Models of Industry
Dynamics∗
Takashi Kamihigashia,b and John Stachurskic
March 7, 2014
a IPAG Business School, Paris, Franceb Research Institute for Economics and Business Administration, Kobe University
c 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

As is well-known, coupling from the past (CFTP) algorithms can generate perfect samples from otherwise intractable stationary distributions (see, e.g., Propp and Wilson∗ We

are grateful for financial support from ARC Discovery Outstanding Researcher Award
DP120100321 and the Japan Society for the Promotion of Science. Email: tkamihig@rieb.kobe-u.ac.jp,
john.stachurski@anu.edu.au.

1

(1996)). This makes them an attractive alternative to ordinary forward simulation,
where errors are typically difficult to assess. Parallelization is simple in most cases,
making it possible to rapidly generate large numbers of independent and identically
distributed draws. Because the resulting sample is IID and exact, it can be used toobtain unbiased estimates and confidence intervals for moments and distributions ofinterest.
The CFTP technique is often used for models with large but discrete state spaces,
where standard methods for computing stationary distributions are infeasible.1 Morerecently, researchers have developed techniques for implementing CFTP methods incontinuous state settings. Murdoch and Green (1998) showed that CFTP can in principle be used in continuous state settings when the underlying Markov process satisfies Doeblin’s condition. This condition requires the existence of a nonnegative lower
bound function that (a) integrates to a positive value, (b) depends only on the next
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 exact sampling.
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.
In this paper we study regenerative processes associated with industry dynamics that
1 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).

2

take the form

Φt+1 = g(Φt , Ut+1 )1{Φt ≥ x } + Zt+1 1{Φt < x }.

(1)

Here Φt is the state variable taking values in a closed interval, x is a point in the interiorof that interval, g is a given function, 1{·} is an indicator function and {Ut } and { Zt }are IID processes. The function g is monotone increasing in its first argument, and suchthat the state is driven below x with positive probability in finite steps. When Φt < x
the process regenerates, being drawn from the distribution of Zt+1 . If this distribution is stochastically dominated by that of g( x, Ut+1 ), then the process is stochasticallymonotone in the usual sense (see, e.g., Corcoran and Tweedie (2001)). As we discuss
below, this is not the case in many models of interest. The purpose of this paper is todevelop a CFTP algorithm that is designed to produce exact draws from the stationary
distribution even in the non-monotone case.
Regenerative models of the form (1) appear in many models of industry dynamics—
typical examples include Hopenhayn (1992), Melitz (2003) and Cooley and Quadrini(2001). Regeneration plays an essential role in the theory of industrial organization,
where Schumpeter’s notion of “creative destruction” summarizes the idea that new
and more productive firms replace older ones, rejuvenating overall economic activity(Schumpeter (1942)). In this context, the model in (1) has the following interpretation:
A large number of firms produce in a given industry. While incumbent, their productivity evolves according to Φt+1 = g(Φt , Ut+1 ), but when their productivity drops
below a threshold level x they exit, and are replaced by a new firm with productivity
drawn according to the distribution of Zt . The value x is obtained by the solution toa dynamic programming problem, taking in to account revenue, costs and so on (although in this paper we take it as given). The stationary distribution of the model canbe regarded as representing the time-invariant cross-sectional distribution of productivity in the industry.
The assumption that g is increasing in its first argument is based on the idea that if firm
A has higher productivity than firm B today, then firm A is expected to have higher
productivity than firm B tomorrow provided that both firms survive. On the other
hand, if firm B does not survive and the productivity of firm A is close to the thresholdlevel x, then a new entrant replacing firm B is likely to be more productive than firm
A. This characteristic of entry and exit makes the entire process non-monotone.
Other applications for the regenerative model (1) can be found in various intertemporaldecision problems. For example, consider the problem of optimal replacement of a partor machine, the performance of which degrades stochastically over time. Typically, the3

solution to the planning problem involves replacement when some measure of performance falls below a certain threshold (e.g., Rust (1987)). A similar idea is rebooting or
restarting operating systems, as discussed by Cotroneo et al. (2014). The performanceof these systems can also degrade over time as a result of memory leaks, software conflicts and so on. An essential feature of these models is again non-monotonicity: thevery purpose of replacement and rebooting is to break monotonicity and rejuvenate
performance.
While we do not assume that the entire process (1) is monotone, we do assume that g isincreasing in its first argument, as discussed above, and we heavily exploit this property in our algorithm. We show that the algorithm terminates successfully in finite timewith probability one by using both the monotonicity of productivity for incumbents
and the regenerative property introduced by new entrants. Our algorithm is distinct
from Murdoch and Green’s method discussed above (Murdoch and Green (1998)), inthat it does not use Doeblin’s condition, and does not require explicit knowledge of thetransition density.2 As long as one can simulate the overall Markov process, one cansample exactly from the stationary distribution using the algorithm.3

2

Preliminaries

In 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.4
2 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.
3 As usual, exactness is modulo the errors associated with floating point arithmetic, which cannot be
avoided.
4 F n (· | φ ) is the conditional distribution for productivity after n periods, given current productivity
φ.

4

R
We let P denote the corresponding stochastic kernel. That is, P(φ, A) := A F (dφ0 | φ)
for φ ∈ S and A ∈ B, where B represents the Borel sets on [0, 1]. Incumbent firms
exit the industry whenever their current productivity falls below a reservation value
xt . Letting Mt be the mass of entrants at time t and ν be the Borel probability measurefrom which the productivity of entrants is drawn, the sequence of firm distributions
R{µt } on S satisfies µt+1 ( A) = 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 probability5 measure µ satisfyingµ( A) =Z

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

(2)

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

Q(φ, A)µ(dφ)

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

(3)

(4)

Equation (3) 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.
It is straightforward to produce an ergodic Markov process suitable for simulation suchthat its stationary distribution coincides with the distribution µ in (3). In essence, weneed a method for sampling from the stochastic kernel Q. The first step is to simulate
from the conditional distribution P(φ, ·) = F (· | φ). In particular, we seek a randomvariable U and a function g such that 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 | φ).6 Now consider the process {Φt } defined byΦt+1 = g(Φt , Ut+1 )1{Φt ≥ x } + Zt+1 1{Φt < x }(5)where {(Ut , Zt )} is IID with D( Zt ) = ν and D(Ut ) = Uniform[0, 1]. Comparing (4)and (5), it can be seen that {Φt } is a Markov process with stochastic kernel Q.
5 We focus only on normalized measures, since other cases are just scalar multiples.
6 Here F −1 (· | φ ) is the generalized inverse of F (·| φ ). That is, F −1 ( u | φ ) : = inf{ z : F ( z | φ ) ≥ u }.

5

3

Exact SamplingLet {(Ut , Zt )}t∈Z be an infinite sequence of IID shocks indexed on Z and with each 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 (5) becomes Φt+1 = ht+1 Φt := ht+1 (Φt ).
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 }(6)

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Λ T := {h0 · · · h−T +k+1 Z−T +k : k = 1, . . . , σT + 1}.

6

(7)

1

x

0 10

9

8

7

5

6

4

3

2

1

0

Figure 1: Tracking processes with T = 10 and σT = 8
The proof of lemma 3.1 is given in section 4. The intuition for the result can be obtainedby considering figure 1. 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
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 2 and 3 show simulations with successful and unsuccessful coalescence re7

1

0 60

50

40

30

20

10

0

50

40

30

20

10

0

1

0 60

Figure 2: Successful coalescence from T = 50

1

0 60

50

40

30

20

10

0

50

40

30

20

10

0

1

0 60

Figure 3: Failure of coalescence8

spectively. In each figure, the top panel shows only the tracking processes. (As with
figure 1, 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 2, there is only one final, coalesced value, and Φ0 must take this value. Hence Φ0 is observed. On the other hand,in figure 3, Φ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
the seed, which determines the values of the shock processes, and then taking T larger
and larger until coalescence occurs. Algorithm 1 gives details. The algorithm terminates with an exact draw from µ. Replication with independent shocks will generate
independent draws.
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 Φ0 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 Φ0 ;At this stage we do not know that the algorithm will terminate with probability one.
This issue is central to the correctness of the algorithm because, as discussed above,
the way we avoid conditioning is by revealing the value of Φ0 for every random seed.
We now show that probability one termination in finite time holds under the following
9

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 4: Stationary density
condition, which is satisfied by many standard distributions.
Assumption 3.1. If G ⊂ S is a nonempty open set, then ν( G ) > 0.
Proposition 3.1. Let T ∗ := min{ T ∈ N : Σ T is nonempty and Λ T is a singleton }. If assumption 3.1 holds, then there exists a γ ∈ (0, 1) such that P{ T ∗ > t} ≤ γt . In particular,
P{ T ∗ < ∞} = 1.
Note that proposition 3.1 not only gives probability one occurrence, but also provides
the geometric rate P{ T ∗ > t} = O(γt ). The proof of proposition 3.1 is given in section 4.
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.7
Figure 4 shows the density computed from 36,000 observations combined with a standard nonparametric kernel density estimator (using a Gaussian kernel). Figure 5 shows
a 95% confidence set for the cumulative distribution function corresponding to µ, basedon the same observations and calculated using the Kolmogorov distribution of the sup
7 Our workstation has a 2.67GHz Intel CPU and 4 gigabytes of RAM.

10

1.0
0.8
0.6
0.4
0.2
0.00.0

0.2

0.4

0.6

0.8

1.0

Figure 5: 95% confidence set for the stationary distribution
norm 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.

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 ,11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 .
Since j ∈ {0, . . . , σ}, the right-hand side is an element of Λ T . This completes the proof.

Proof of proposition 3.1. Let n be an integer such that F n ( x | 1) > 0, existence of which isdue 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}.

(8)

To prove that (8) 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
12

and so on. Together, these equalities give (8). As a consequence, we conclude that Σn+1
is nonempty and Λn+1 is a singleton whenever E1 occurs, and, more generally, Σnj+1
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 result stated in the proposition.

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.
Cotroneo, Domenico, Roberto Natella, Roberto Pietrantuono, and Stefano Russo (2014)“A survey of software aging and rejuvenation studies,” ACM Journal on Emerging
Technologies in Computing Systems (JETC), Vol. 10, p. 8.
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.

13

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.
Schumpeter, Joseph (1942) “Creative destruction,” Capitalism, socialism and democracy.

14