Geometrically stopped Markovian random growth
processes and Pareto tails
Brendan K. Beare and Alexis Akira Toda

arXiv:1712.01431v1 [econ.EM] 5 Dec 2017

Department of Economics, University of California, San Diego
December 6, 2017

Abstract
Many empirical studies document power law behavior in size distributions of economic interest such as cities, firms, income, and wealth. One
mechanism for generating such behavior combines independent and identically distributed Gaussian additive shocks to log-size with a geometric
age distribution. We generalize this mechanism by allowing the shocks
to be non-Gaussian (but light-tailed) and dependent upon a Markov state
variable. Our main results provide sharp bounds on tail probabilities and
simple formulas for Pareto exponents. We present two applications: (i) we
show that the tails of the wealth distribution in a heterogeneous-agent dynamic general equilibrium model with idiosyncratic endowment risk decay
exponentially, unlike models with investment risk where the tails may be
Paretian, and (ii) we show that a random growth model for the population
dynamics of Japanese prefectures is consistent with the observed Pareto
exponent but only after allowing for Markovian dynamics.
Keywords: exponential tails, Gibrat’s law, Pareto tails, power law,
random growth, Tauberian theorem.
JEL codes: C46, C65, D30, D52, D58, R12.

1

Introduction

In this paper we are concerned with the tail behavior of geometric sums of the
form
T
X
WT =
Xt ,
(1.1)
t=1

We thank Alberto Bisin, Xavier Gabaix, Émilien Gouin-Bonenfant, Makoto Nirei, Werner
Ploberger, Jonathan Weinstein, and seminar participants at Queen’s University, UCSB, UCSD,
the University of Toronto, Washington University in St. Louis, and the 2017 Summer Workshop
in Economic Theory at the Cowles Foundation for helpful comments and suggestions.

1

where {Xt }∞
t=1 is a sequence of random innovations and T is a geometrically distributed random variable. The innovations {Xt }∞
t=1 may depend contemporaneously on a time-homogeneous Markov state variable (i.e., Xt is a “hidden Markov
process”; see Definition 3.1 below). We may view WT in (1.1) as a geometrically
stopped random walk with possibly serially dependent innovations. The main result of our paper shows that under quite general conditions the distribution of WT
has exponential tails. Further, we provide a simple formula characterizing the tail
decay rates. For instance, when {Xt }∞
t=1 is independent and identically distributed
(iid) and T has mean 1/p, the upper tail exponent α of WT solves the equation
E[eαX ] =

1
.
1−p

(1.2)

In more general settings, the left-hand side of (1.2) is replaced with the spectral
radius of a matrix determined by the law of the innovations conditional on the
state variable and its transition probabilities; see equation (3.2) below.
Our problem is motivated by the observed characteristics of size distributions
in economics and other fields, where variables of interest are often known to exhibit
power law behavior in the upper tail.1 A variety of explanations for the emergence
of these power laws have been proposed in the economics literature. One particular
mechanism for generating power laws proposed by Reed (2001), and to some extent
anticipated by Wold and Whittle (1957) and Simon and Bonini (1958), combines
two main ingredients: Gibrat’s law of proportional growth (Gibrat, 1931) and an
exponential age distribution.2 Suppose that the size St of an individual unit at
age t ≥ 0 follows a geometric Brownian motion initialized at some fixed S0 , so
that
dSt = µSt dt + σSt dBt ,

(1.3)

where Bt is a standard Brownian motion and µ and σ are the drift and volatility
parameters. If the distribution of the age T of a unit randomly selected at a
given point in time is exponential with parameter η > 0, so that it has probability
density function (pdf) fT (t) = ηe−ηt for t ≥ 0, then the size of a randomly selected
unit is given by ST , our geometric Brownian motion evaluated at an exponentially
1

Such variables include city sizes (Gabaix, 1999; Reed, 2002; Giesen et al., 2010), firm
sizes (Axtell, 2001; Luttmer, 2007), and household income (Pareto, 1896; Reed, 2003;
Reed and Jorgensen, 2004; Toda, 2011, 2012), consumption (Toda and Walsh, 2015; Toda,
2017b), and wealth (Klass et al., 2006; Vermeulen, 2017).
2
There is a large empirical literature documenting that Gibrat’s law of proportional growth
is a good first approximation; see Sutton (1997) and the references therein. The exponential
age distribution has drawn much less attention but there is some evidence both for firms (Coad,
2010) and cities (Giesen and Suedekum, 2014).

2

distributed time. Reed (2001) showed that the pdf of this quantity is given by

fST (s) =





αβ
S −β sβ−1
α+β 0

for 0 ≤ s < S0

αβ
S α s−α−1
α+β 0

for S0 ≤ s < ∞,

where ζ = α, −β are the positive and negative roots of the quadratic equation


σ2 2
σ2
ζ − η = 0.3
ζ + µ−
2
2

(1.4)

He dubbed this distribution the double Pareto distribution. Reed’s mechanism is
notable in that it generates power law behavior in not only the upper tail of a
distribution, but also the lower tail: we have
lim Asα P(ST > s) = 1 and

s↑∞

lim Bs−β P(ST < s) = 1
s↓0

(1.5)

for some positive constants A and B. This is the case even though the distribution
of a geometric Brownian motion evaluated at a fixed point in time is lognormal and
therefore has tails decaying more rapidly than a power law. Reed’s mechanism
has recently been applied in economics to characterize the tail behavior of size
distributions in heterogeneous-agent models in continuous-time with Brownian
shocks.4
Does Reed’s mechanism also generate power law tails when applied to more
general stochastic processes? Given that power law distributions are empirically
so common, and that realistic alternatives to Brownian motion may involve nonGaussian, non-independent increments, it is natural to conjecture that more general random growth processes with geometric age distributions give rise to power
law tails, or exponential tails after taking the logarithm. Some affirmative results in this direction are available in the physics literature for the geometric sum
(1.1) in the special case where {Xt }∞
t=1 is iid. Manrubia and Zanette (1999) give
a heuristic derivation of (1.2), which corresponds to equation (10) in their paper,
and provide supporting evidence from numerical simulations. Reed and Hughes
(2002) observe that the tails of their “killed discrete multiplicative process” are
characterized by (1.2). In the economics literature, (1.2) appears in Proposition 5
Note that assuming µ > 0 and letting σ → 0, we obtain α = µη , which is exactly the result
of Wold and Whittle (1957). (Their paper contains a typographical error in the formula for α
on the first page, in which the numerator and the denominator are flipped.)
4
See, for example, Toda (2014), Toda and Walsh (2015, 2017), Arkolakis (2016),
Benhabib et al. (2016), and Gabaix et al. (2016).
3

3

of Nirei and Aoki (2016), who appeal to Manrubia and Zanette (1999) to characterize the tail exponent of the wealth distribution in a heterogeneous-agent model
with iid shocks.
The main results of our paper, Theorems 3.1–3.4 below, provide a formal justification for the formula (1.2) in the iid case and extend it to a wider class of
processes. The key to proving them is a Tauberian theorem due to Nakagawa
(2007),5 which we discuss in Section 2 and Appendix A. Nakagawa’s theorem provides sharp bounds on the tail probabilities of a random variable whose Laplace
transform has a pole at its abscissa of convergence. In the case of a simple pole,
which is the relevant one for our purposes, those bounds depend on the residue of
the pole and the distance to other singularities along the axis of convergence. To
obtain the residue in the non-iid case, we rely on the Perron-Frobenius theory of
nonnegative irreducible matrices (Appendix B) and some results on simple poles
of matrix pencil inverses (Appendix C). When the distribution of Xt conditional
on the current Markov state is not concentrated on an evenly spaced grid, S = eWT
has Pareto tails as in (1.5). Without this condition, the limits in (1.5) may not
exist, but we have sharp bounds for the corresponding limits superior and inferior.
Those bounds imply (Remark 2.3) that the tail probabilities satisfy
log P(S > s)
= −α and
s↑∞
log s
lim

lim
s↓0

log P(S < s)
= β,
log s

which is a weaker property than (1.5). In Theorem 3.5, we derive comparative
statics for the Pareto exponent with respect to perturbations in lifespan, growth,
volatility, and persistence.
As one application of our results, in Section 4 we characterize the tail behavior of the wealth distributions in heterogeneous-agent macroeconomic models. We
consider two cases, economies with idiosyncratic endowment and investment risks,
both with Markovian shocks. In the former case, we prove that the wealth distribution has exponential tails, which are thin. This theoretical result explains the numerical findings in the literature that the so-called Bewley (1983)-Huggett (1993)Aiyagari (1994) models with labor income risk alone have difficulty in matching
the (empirically fat) upper tail of the wealth distribution (see Benhabib et al.,
5

Reed and Hughes (2003, p. 588) suggested that Tauberian theorems may be useful to characterize tail probabilities in the related class of Galton-Watson branching processes. The Online
Appendix of Gabaix et al. (2016) appeals to a Tauberian theorem of Mimica (2016, Corollary
1.4), slightly more general than a result of Nakagawa (2007, Theorem 3), to characterize the
tail decay rate of the wealth distribution. We appeal to a sharper result (Theorem 2.1 below)
essentially obtained by Nakagawa (2007, Theorem 5*).

4

2017, for an excellent discussion of this issue). The latter case with idiosyncratic
investment risk has already be known to generate Pareto tails under specific conditions (typically iid Gaussian shocks). We extend those results by permitting
non-Gaussian shocks dependent on a Markov state variable and thus considerably
expand the toolkit for applied economists.
As a second application of our results, in Section 5 we ask whether the historical population dynamics of Japanese prefectures are consistent with a Pareto
cross-sectional size distribution. Using the population of 46 Japanese prefectures
(excluding Okinawa) since 1873, we document that the power law emerged only
around 1920 and the Pareto exponent has decreased over the century from about
3 in 1920 to about 1.3 in 2015 (except for a brief break during World War II).
We fit a random growth model to the panel of prefecture populations. Using an
iid random growth model, the implied Pareto exponent is 56.7, which is far larger
than the value directly estimated from the cross-section (1.3). However, we show
that by estimating a Markov switching model with three states, the implied Pareto
exponent becomes 1.6, which is close to what we see in the data. This shows that
incorporating Markovian dynamics into random growth models can be crucial for
explaining tail exponents observed in data. Our theorems provide a technical tool
to study such models.

1.1

Related literature

Important early contributions on generative mechanisms for power laws were made
by Champernowne (1953), who proposed a model of income dynamics in which
a power law emerges as a steady state equilibrium. Wold and Whittle (1957)
found that a constant growth rate of wealth together with an exponential age
distribution generates a Pareto upper tail, which anticipates the mechanism of
Reed (2001). Simon and Bonini (1958) observed that Gibrat’s law in combination
with exponential growth in the number of firms could lead to a power law in
the upper tail of the firm size distribution. Kesten (1973) studied the random
difference equation Xt = At Xt−1 + Bt and showed that the stationary distribution
exhibits a power law tail.6
6

da Saporta (2005), Roitershtein (2007), and others have extended the Kesten (1973) result
to Markovian environments using renewal theory. Our results are different primarily because
(i) growth is proportional (Bt = 0), consistent with Gibrat’s law, and (ii) we study the tails
of a geometrically stopped nonstationary process, not the tails of the invariant distribution
of a stationary (Kesten) process. Furthermore, our results are more generally applicable in
economics because proportional growth arises in a wide variety of dynamic models, whereas the
Kesten process is more rare.

5

More recently, a paper by Gabaix (1999) providing an explanation for Zipf’s law
in city sizes sparked renewed interest in power laws amongst economists. Gabaix
showed that augmenting the geometric Brownian motion law for city size evolution
with a lower reflecting boundary generates a power law exponent slightly above one
in the upper tail of the city size distribution, consistent with Zipf’s law. Note that
this generative mechanism is distinct from that of Reed (2001). Survey papers by
Gabaix (2009, 2016) discuss much of the subsequent economics literature. Another
survey paper by Mitzenmacher (2004) discusses generative mechanisms for power
laws proposed across a wider range of disciplines, including biology, computer
science, networks, and physics. A more recent survey by Benhabib and Bisin
(2017) focuses on the wealth distribution.
Given the long tradition of random growth models, a little historical digression
may be justified. By the early 20th century, it was recognized as a puzzle amongst
statisticians that variables of empirical interest that were considered to be accumulating numerous small independent shocks frequently exhibited skewness or excess
kurtosis, despite the central limit theorem. A simple explanation was suggested
by Kapteyn (1903): if small shocks accumulate to a function F (x) of a quantity
of interest x, rather than x itself, then one can generate skewed distributions by
applying the change of variable formula to the normal distribution. Gibrat (1931,
ch. 5) applied this argument to F (x) = log x and arrived at his celebrated “la loi
de l’effect proportionnel” (the law of proportional effect). In Chapter 6, Gibrat
elaborates further in the context of the firm size distribution and mentions that
necessary and sufficient conditions for obtaining “formula A” (which is essentially
the lognormal distribution) are:
(i) “Les causes de fluctuation du personnel sont nombreuses” (there is a large
number of shocks to the fluctuation of employees),
(ii) “Leur effet relatif sur le nombre d’ouvriers (ou leur effet absolut sur le logarithme), ne dépend pas de ce nombre d’ouvriers.” (the relative effect on the
number of workers (or the absolute effect on the logarithm) does not depend
on the number of workers),
(iii) “L’effet de chaque cause de fluctuation est petit vis-à-vis de l’effet de toutes”
(the effect of each shock is small relative to the total effect).
Essentially, he is applying the central limit theorem to the shocks to the logarithm
of the firm size measured by the number of employees (and cites Lindeberg, 1922,

6

for a justification).7
Gibrat’s argument has one pitfall, as pointed out by Kalecki (1945): with a random growth model with infinitely lived agents, the cross-sectional distribution is
approximately lognormal but the log mean and variance increase linearly overtime,
and hence a stationary distribution does not exist. One of the very first solutions
to this problem was to introduce birth and death. In particular, Rutherford (1955)
assumes an exponential age distribution. Since the income shocks are Gaussian
in his model, it is a discrete-time version of the model of Reed (2001). However,
Rutherford (1955) did not characterize the tail behavior; for this we had to wait
until Wold and Whittle (1957) and Reed (2001).

2

Exponential tails via the Laplace transform

Theorems that allow us to deduce limit properties of a probability distribution or
other function from limit properties of its Laplace transform are called Tauberian
theorems. In this section we state a version of a Tauberian theorem of Nakagawa
(2007), which will be used to prove our main results on the tail probabilities of
geometric sums of hidden Markov processes. First we briefly review the Laplace
transform. For more details see Widder (1941) and Lukacs (1970, ch. 7).
Given a cumulative distribution function (cdf) F , let
M(s) =

Z

∞

esx dF (x) ∈ R ∪ {∞}
−∞

be its moment generating function (mgf). Since esx is convex in s, the domain
I = {s ∈ R : M(s) < ∞} is convex, and hence an interval. Clearly 0 ∈ I, so there
exist boundary points −β ≤ 0 ≤ α of I (which may be 0 or ±∞). The numbers
α, −β are called the right and left abscissae of convergence.
We obtain the Laplace transform of F by extending the domain of M into the
complex plane. Suppose that −β < α and let z = s + it ∈ C. By the definition of
the Lebesgue-Stieltjes integral,
M(z) =

Z

∞

ezx dF (x)

(2.1)

−∞

7

Note that nowhere in Gibrat’s original argument is the assumption of Gaussian shocks. Using
the geometric Brownian motion to represent a random growth process is merely a mathematical
convenience, which may not be supported empirically. For example, Arata (2015, ch. 2) provides
evidence that the changes in log firm size are not Gaussian and argues that the more general
class of Lévy processes may be empirically preferable. However, such processes were introduced
by Lévy (1937) subsequent to the publication of Gibrat (1931).

7

exists if and only if
Z

∞
zx

|e | dF (x) =

−∞

Z

∞
(s+it)x

e

dF (x) =

−∞

Z

∞

esx dF (x) = M(s) < ∞.

−∞

Therefore by the definition of α, −β, the value M(z) is well-defined for z ∈ C with
Re z ∈ I. Let S = {z ∈ C : −β < Re z < α}. Using the dominated convergence
theorem, it is easy to see that M(z) is holomorphic on S, which is called the strip
of holomorphicity (Figure 2.1). The lines Re z = α, −β comprising the boundary
of S are called the right and left axes of convergence. In this paper we refer to
(2.1) as the (two-sided) Laplace transform of F , or of any real random variable X

Im z

S = Strip of
holomorphicity

0

Right axis of convergence

Left axis of convergence

with cdf F . We also use “Laplace transform” and “moment generating function”
interchangeably.

Re z

Left abscissa of
convergence

Right abscissa of
convergence

Figure 2.1: Region of convergence of the Laplace transform M (z) =

R∞

zx
−∞ e dF (x).

There is a close relationship between the tail probabilities of a cdf and the
abscissae of convergence of its Laplace transform. In general, we have
lim sup
x→∞

1
log P(X > x) = −α,
x

lim sup
x→∞

1
log P(X < −x) = −β
x

(2.2)

whenever the relevant abscissa is nonzero (Widder, 1941, pp. 42–43, 241; Lukacs,
1970, p. 194). Moreover, each abscissa is a singularity of the Laplace transform
(Widder, 1941, p. 58). A Tauberian theorem of Nakagawa (2007) establishes a
tighter relationship between the abscissae of convergence and the tail probabilities
that depends on the nature of the singularities at the abscissae: if either singularity
is a pole then the relevant limit superior in (2.2) may be replaced with an ordinary
limit. In addition, bounds for the limits superior and inferior of eαx P(X > x) and
8

eβx P(X < −x) become available which depend upon the order of the poles at the
abscissae, the leading Laurent coefficients, and the location of other singularities
along the axes of convergence. For our purposes, it will be enough to consider the
case where the singularities at the abscissae are simple poles. The following result
is similar to Theorem 5* of Nakagawa (2007), specialized to the case of a simple
pole. We provide a proof in Appendix A (see Remark 2.5 below).
Theorem 2.1. Let X be a real random variable and M(z) = E[ezX ] its Laplace
transform, with right abscissa of convergence α satisfying 0 < α < ∞, and strip
of holomorphicity S. Fix A > 0, and let B be the supremum of all b > 0 such that
the function M(z) + A(z − α)−1 can be continuously extended from S to the set
Sb+ = S ∪ {z ∈ C : z = α + it, |t| < b} .

(2.3)

Suppose that B > 0. Then we have
2πA/B
2πA/B
≤ lim inf eαx P(X > x) ≤ lim sup eαx P(X > x) ≤
, (2.4)
x→∞
−1
1 − e−2πα/B
x→∞

e2πα/B

where the bounds should be read as A/α if B = ∞. These bounds are sharp.
Remark 2.1. Applying Theorem 2.1 to the random variable −X can yield bounds
analogous to (2.4) for the lower tail probabilities. Specifically, given A > 0, if the
left abscissa −β satisfies −∞ < −β < 0 and B is the supremum of all b > 0 such
that M(z) − A(z + β)−1 can be continuously extended from S to the set
Sb− = S ∪ {z ∈ C : z = −β + it, |t| < b} ,
then, if B > 0, we have the bounds
2πA/B
2πA/B
≤ lim inf eβx P(X < −x) ≤ lim sup eβx P(X < −x) ≤
.
x→∞
−1
1 − e−2πβ/B
x→∞

e2πβ/B

Remark 2.2. To see why the bounds in (2.4) are sharp, let X be a geometrically
distributed random variable with mean 1/p ∈ (0, ∞). The Laplace transform of
X is given by
zX

M(z) = E[e

]=

∞
X

p(1 − p)n−1 ezn =

n=1

pez
1 − (1 − p)ez

on its strip of holomorphicity S. The zeros of the denominator 1−(1−p)ez are zk =
α + 2πik, k ∈ Z, where α = − log(1 − p) > 0. The right abscissa of convergence
9

is therefore α and the strip of holomorphicity is S = {z ∈ C : Re z < α}. The
zeros zk are singularities of M(z) lying along its axis of convergence Re z = α. By
l’Hôpital’s rule we have
pez + (z − α)pez
p
(z − α)pez
=
−
lim
=−
,
z
z
z→α
z→α 1 − (1 − p)e
(1 − p)e
1−p
lim

which shows that the singularity at α is a simple pole with residue −p/(1 − p).
The singularities along the axis of convergence are separated by gaps of 2πi, so
the assumptions of Theorem 2.1 are satisfied with A = p/(1 − p) and B = 2π, and
thus (2.4) holds with the lower and upper bounds
2πA/B
= 1,
−1

e2πα/B

2πA/B
1
=
.
−2πα/B
1−e
1−p

On the other hand, if n ≤ x < n + 1 for some n ∈ N, then
P(X > x) =

∞
X

p(1 − p)k−1 = (1 − p)n = (1 − p)⌊x⌋

k=n+1

by direct computation. Therefore
eαx P(X > x) = (1 − p)−x P(X > x) = (1 − p)⌊x⌋−x ,
which oscillates between 1 and 1/(1 − p) as x → ∞. The bounds in (2.4) are
therefore sharp.
Remark 2.3. When the bounds in (2.4) are satisfied, for any ǫ > 0 there exists
x0 < ∞ such that
2πA/B
2πA/B
− ǫ ≤ eαx P(X > x) ≤
+ǫ
−1
1 − e−2πα/B

e2πα/B

for all x ≥ x0 . Taking logarithms, dividing by x, and letting x → ∞, we obtain
1
log P(X > x) = −α,
x→∞ x
lim

(2.5)

which improves upon the characterization of right tail probabilities in (2.2).
Remark 2.4. In the case B = ∞, which corresponds to the simple pole at α
being the unique singularity on the axis of convergence, we may rewrite (2.4) as
lim sα P(S > s) =

s→∞

10

A
,
α

(2.6)

where S = eX . In this sense, the right tail probabilities of S are Paretian. Property
(2.6) implies regular variation of the right tail probabilities of S with index −α,
under which sα P(S > s) is required to be slowly varying at infinity. Property
(2.6) also implies (2.5), which will be satisfied for any B > 0. On the other hand,
(2.5) does not imply (2.6). For instance, if X is a geometric random variable as
in Remark 2.2, then (2.5) holds but (2.6) fails to hold because
1 = lim inf sα P(S > s) < lim sup sα P(S > s) =
s→∞

s→∞

1
.
1−p

Remark 2.5. Aside from the fact that it deals only with nonnegative random
variables (which is unimportant), Theorem 5* of Nakagawa (2007) is more general
than Theorem 2.1 above, as it provides bounds on tail probabilities that apply
when the pole at the right abscissa of convergence is of arbitrary order. This
additional generality makes the proof quite complicated—and in fact it is given
in detail only for the case of a second order pole. For the reader’s convenience, in
Appendix A we provide a simpler proof of Theorem 2.1 that applies in our more
restrictive setting. Following Nakagawa, we draw on a technique introduced by
Graham and Vaaler (1981) to prove a refinement of the Wiener-Ikehara Tauberian
theorem, a helpful account of which may be found in the monograph of Korevaar
(2004, ch. 5). Theorem 2.1 is in fact very similar to the Tauberian theorem proved
by Graham and Vaaler (1981), with the only real difference being that it concerns
exponentially decaying tail probabilities rather than exponential growth.

3
3.1

Main results
Tail behavior of geometric sums

In this section we characterize the tail behavior of the geometric sum (1.1) when
{Xt }∞
t=1 is a hidden Markov process.
Definition 3.1 (Hidden Markov process). A sequence of real-valued random
variables {Xt }∞
t=1 is called a hidden Markov process if there exists a natural
number N, a time-homogeneous Markov process {Jt }∞
t=0 taking values in N =
{1, . . . , N }, and an iid sequence of RN -valued random variables {Yt }∞
t=1 (where
⊤
Yt = (Y1t , . . . , YN t ) ) such that Xt = YJt ,t for all t ∈ N.
Example 3.1. If {Xt }∞
t=1 is iid then it is a hidden Markov process with N = 1.

11

Example 3.2. If {Xt }∞
t=1 is a finite-state time-homogeneous Markov chain taking the values x1 , . . . , xN then it is a hidden Markov process: just set Yt ≡
(x1 , . . . , xN )⊤ (a constant vector).
Throughout this section, let {Xt }∞
t=1 be a hidden Markov process with the
∞
underlying states denoted by {Jt }t=0 . Let Π = (πnn′ ) be the transition probability
matrix, where πnn′ = P(J1 = n′ | J0 = n). We assume that Π is irreducible.
Define the set



I := s ∈ R : (∀n ∈ N ) E esX1 J1 = n < ∞ ,

(3.1)

and for s ∈ I, let D(s) be the N × N diagonal matrix with n-th diagonal entry


equal to E esX1 J1 = n . Let Wt be the cumulative sum of Xt , so W0 = 0 and

Wt = Wt−1 + Xt for t ∈ N. Finally, let T be a geometrically distributed random
variable with mean 1/p and independent of the X, J processes.
The following theorem is our main result. It shows that under empirically
plausible conditions the geometric sum WT has exponential tails, and characterizes
their decay rates in terms of Π, D and p. That characterization is illustrated in
Figure 3.1. For a square matrix A, we let ρ(A) denote the spectral radius of A.
Theorem 3.1. The set I is convex and contains zero. There can be at most one
positive value of s ∈ I that solves the equation
ρ(ΠD(s)) =

1
.
1−p

(3.2)

Suppose that a positive solution α to (3.2) exists in the interior of I. Then
lim

w→∞

1
log P(WT > w) = −α.
w

(3.3)

Similarly, there can be at most one negative value of s ∈ I that solves equation
(3.2). Suppose that a negative solution −β exists in the interior of I. Then
1
log P(WT < −w) = −β.
w→∞ w
lim

(3.4)

Proof. We divide our proof into five steps, stated in italics as claims to be proved.
Step 1. The spectral radius ρ(ΠD(s)) is convex in s ∈ I. The number α in the
statement of Theorem 3.1, if it exists, is unique.

12

ρ(ΠD(s))

1
1−p

1

0

−β
Left endpoint of I

α

s

Right endpoint of I

Figure 3.1: Determination of α and β.

For each n ∈ N define In =
Hölder’s inequality, we obtain




s ∈ R : E esX1 J1 = n < ∞ .

Applying



E e((1−θ)s1 +θs2 )X1 J1 = n


1−θ

θ
≤ E es1 X1 J1 = n
E es2 X1 J1 = n
< ∞ (3.5)

for all θ ∈ (0, 1), all s1 , s2 ∈ In , and all n ∈ N . This shows that each In is convex
T
and contains zero, implying that I = n∈N In is convex and contains zero. It

also shows that each diagonal entry of the matrix D(s) is a log-convex function of
s ∈ I. Letting ⊙ denote the Hadamard (entry-wise) product, collecting (3.5) into
a matrix we obtain
ΠD((1 − θ)s1 + θs2 ) ≤ ΠD(s1 )(1−θ) D(s2 )(θ) = (ΠD(s1 ))(1−θ) ⊙ (ΠD(s2 ))(θ)
entry-wise for all θ ∈ (0, 1) and all s1 , s2 ∈ I, where for a nonnegative matrix
A = (amn ) we define A(θ) = (aθmn ). By Propositions B.2(ii) and B.1, we obtain
ρ(ΠD((1 − θ)s1 + θs2 )) ≤ ρ (ΠD(s1 ))(1−θ) ⊙ (ΠD(s2 ))(θ)
≤ ρ(ΠD(s1 ))1−θ ρ(ΠD(s2 ))θ ,



which shows that ρ(ΠD(s)) is log-convex (and hence convex).
To show the uniqueness of α, suppose on the contrary that there are two
numbers 0 < α1 < α2 such that λ(s) := ρ(ΠD(s)) equals 1/(1 − p) for s = α1 , α2 .
Let θ = α1 /α2 ∈ (0, 1). Since λ(0) = ρ(ΠD(0)) = ρ(Π) = 1 (a nonnegative matrix
whose rows sum to one has spectral radius one—see e.g. Horn and Johnson, 1985,

13

p. 547) and 1 < 1/(1 − p), it follows from the convexity of λ that
1
= λ(α1 ) = λ(θα2 + (1 − θ)0)
1−p
1
θ
+1−θ <
,
≤ θλ(α2 ) + (1 − θ)λ(0) =
1−p
1−p
which is a contradiction.
Define the strip S = {z ∈ C : Re z ∈ I} (Figure 3.2). For each z ∈ S and




n ∈ N we have E ezX1 J1 = n = E e(Re z)X1 J1 = n < ∞, and may therefore

extend the domain of definition of D(s) from I to S by letting D(z) be the N × N


diagonal matrix with n-th diagonal entry equal to E ezX1 J1 = n , a holomorphic
function of z on the interior of S.
ρ(ΠD(z))

Im z

1
1−p

1

−β

S = Strip of
holomorphicity
for D(z)

α

Re z

Figure 3.2: Definition of S.

Step 2. Let e be the N × 1 vector of ones. For t ∈ {0} ∪ N and z ∈ S, the
conditional mgf of Wt given J0 is given by


E ezWt J0 = 1


..
 = (ΠD(z))t e.

.
 

zWt
E e
J0 = N


(3.6)

Consequently, the mgf of Wt is given by E[ezWt ] = ω0⊤ (ΠD(z))t e for z ∈ S, where
ω0 is the N × 1 vector of probabilities P(J0 = n), n = 1, . . . , N.
Equality (3.6) is trivially true for t = 0 since W0 = 0 by definition. To deal
with the case t ≥ 1 we use the law of iterated expectations and the conditional
independence of Wt+1 and J0 given J1 to write
N
 zWt+1
 X


E e
J0 = n =
P(J1 = n′ | J0 = n) E ezWt+1 J1 = n′ .
n′ =1

14

Since J0 , J1 , . . . is a time-homogeneous Markov chain, the law of Wt+1 − X1 given
J1 = n′ is the same as the law of Wt given J0 = n′ . Thus, using the conditional
independence of X1 and Wt+1 − X1 given J1 we have



 

E ezWt+1 J1 = n′ = E ezX1 J1 = n′ E ez(Wt+1 −X1 ) J1 = n′

 

= E ezX1 J1 = n′ E ezWt J0 = n′ .

It follows that


zWt+1

E e

N
 X

 

J0 = n =
P(J1 = n′ | J0 = n) E ezX1 J1 = n′ E ezWt J0 = n′ .
n′ =1



If (3.6) is true for some t ≥ 0, this last equation tells us that E ezWt+1 J0 = n

is equal to the n-th row of ΠD(z)(ΠD(z))t e. Thus (3.6) is true for all t ≥ 0 by
induction.
Define the set


1
Ip = s ∈ I : ρ(ΠD(s)) <
.
1−p

The set Ip is convex and contains zero, due to the previously established convexity
of ρ(ΠD(s)) and the fact that ρ(ΠD(0)) = 1 < 1/(1 − p). Define the strip Sp =
{z ∈ C : Re z ∈ Ip } (Figure 3.3). Let I be the N ×N identity matrix, and for z ∈ S
define the matrix-valued complex function (matrix pencil) A(z) = I−(1−p)ΠD(z).
ρ(ΠD(z))

1
1−p

Im z

1

−β

Sp = Strip of
holomorphicity
for WT

α

Re z

Figure 3.3: Definition of Sp .

Step 3. For z ∈ Sp , A(z) is invertible and the mgf of WT is given by
E[ezWT ] = pω0⊤ A(z)−1 ΠD(z)e.

(3.7)

For z ∈ Sp , we have ρ(ΠD(z)) ≤ ρ(ΠD(Re z)) < 1/(1 − p), where the first
15

inequality follows from Proposition B.2(ii) and the second from the definition
of Sp . Hence by Proposition B.2(i), (1 − p)k (ΠD(z))k decays to zero at an
exponential rate as k → ∞ for any matrix norm k·k and any z ∈ Sp . Hence the
P
k
k
geometric series ∞
k=0 (1 − p) (ΠD(z)) converges on Sp . However,
(I − (1 − p)ΠD(z))

K
X

(1 − p)k (ΠD(z))k = I − (1 − p)K+1 (ΠD(z))K+1 → I

k=0

as K → ∞, so A(z) = I − (1 − p)ΠD(z) is invertible on Sp and
A(z)−1 =

∞
X

(1 − p)k (ΠD(z))k .

k=0

Therefore, using the fact that E[ezWt ] = ω0⊤ (ΠD(z))t e for z ∈ S and hence for
z ∈ Sp , we obtain
zWT

E[e

]=

∞
X

p(1 − p)k−1ω0⊤ (ΠD(z))k e = pω0⊤ A(z)−1 ΠD(z)e

k=1

for z ∈ Sp , which is (3.7).
Step 4. Suppose there exists α > 0 in the interior of I such that ρ(ΠD(α)) =
1/(1 − p). Then A(z)−1 is well-defined and holomorphic on a punctured neighborhood of the singularity at z = α, which is a simple pole.
Suppose there exists α > 0 in the interior of I such that ρ(ΠD(α)) = 1/(1−p).
By assumption Π is nonnegative and irreducible, and D(α) is a diagonal matrix
with positive diagonal elements, so ΠD(α) is also nonnegative and irreducible.
The Perron-Frobenius theorem (Appendix B) therefore implies that the spectral
radius 1/(1 − p) is an eigenvalue of ΠD(α), called the Perron root. Let x, y be
the right and left Perron vectors corresponding to the Perron root. Since x, y are
positive, we can normalize them so that the entries of x, y sum to 1.
Since 1/(1 − p) is an eigenvalue of ΠD(α), we have det A(α) = 0. On the other
hand, since A(z) is invertible for z ∈ Sp ⊂ S, we have det A(z) 6= 0 for some z ∈ S.
Moreover, A(z) and det A(z) inherit from D(z) the property of holomorphicity on
the interior of S. Since det A(z) is nonconstant and holomorphic on the interior
of S with a zero at z = α, that zero must be isolated (Ahlfors, 1979, p. 127),
and therefore A(z) is holomorphic on a neighborhood of z = α with an isolated
point of noninvertibility at z = α. Consequently, Proposition C.1 tells us that
the inverse A(z)−1 is holomorphic on a punctured neighborhood of the singularity
16

at z = α, and that this singularity is a simple pole if and only if the algebraic
and geometric multiplicities of the eigenvalue 1/(1 − p) of ΠD(α) are equal. We
know from the Perron-Frobenius theorem that both multiplicities are one; A(z)−1
therefore has a simple pole at z = α.
Step 5. The right-hand side of (3.7) is well-defined and holomorphic on a punctured neighborhood of the singularity at z = α, which is a simple pole.
Proposition C.1 implies that the residue R of the simple pole α of A(z)−1 is
given by the N × N matrix
R = x(y ⊤ A′ (α)x)−1 y ⊤ =: cxy ⊤ ,

(3.8)

where A′ (z) is the matrix of derivatives of A(z), x, y are the Perron vectors introduced above, and c = (y ⊤ A′ (α)x)−1 is a nonzero scalar. Since x, y are right and
left eigenvectors of ΠD(α) with corresponding eigenvalue 1/(1 − p), it follows that
A(α)R = RA(α) = 0.
By assumption α belongs to the interior of S, and we observed earlier that
D(z) is holomorphic on the interior of S, so we may use the right-hand side of
(3.7) to holomorphically extend E[ezWT ] to a punctured neighborhood of α, which
is a singularity. Let us show that this singularity is a simple pole by showing that
lim (z − α)pω0⊤ A(z)−1 ΠD(z)e

z→α

(3.9)

exists and is nonzero. Since A(z)−1 has a simple pole at z = α with residue R we
know that
lim (z − α)A(z)−1 = R,

z→α

and since D(z) is continuous at z = α we know that limz→α D(z) = D(α). Therefore the limit in (3.9) exists and is equal to pω0⊤ RΠD(α)e. Since RA(α) = 0, we
have
RΠD(α)e = −

1
1
1
R(I − (1 − p)ΠD(α))e +
Re =
Re.
1−p
1−p
1−p

Using (3.8) again, it follows that
(3.9) = pω0⊤

1
p
p
Re = c
(ω0⊤ x)(y ⊤ e) = c
ω0⊤ x 6= 0
1−p
1−p
1−p

(3.10)

since the entries of y sum to 1, ω0 > 0, and x ≫ 0.8 Therefore our holomorphic
8

Given real vectors a = (an ) and b = (bn ), we say that a ≥ b if an ≥ bn for all n, that a > b

17

extension of E[ezWT ] has a simple pole at z = α.
The limit (3.3) now follows from Theorem 2.1 and Remark 2.3. A symmetric
argument establishes the analogous result for the lower tail.
Remark 3.1. Equation (3.2) is the most important equation in our paper. It
provides an implicit characterization of the tail exponents α and β in terms of the
parameters Π, D, and p. In Section 3.3 we use the implicit function theorem to
investigate the response of the tail exponents to perturbations in these parameters.
sX1
], and the
Remark 3.2. If N = 1, so that {Xt }∞
t=1 is iid, then ρ(ΠD(s)) = E[e

tail exponents are determined by the equations
E[eαX1 ] =

1
1−p

and

E[e−βX1 ] =

1
,
1−p

whenever such α and/or −β exist in the interior of I. This explains (1.2).
Remark 3.3. One of the assumptions of Theorem 3.1 is that there exists α > 0
in the interior of I such that ρ(ΠD(α)) = 1/(1 − p). This assumption may not
always be satisfied. As a counterexample, fix α > 0 and suppose that {Xt }∞
t=1 is
iid (N = 1) with an mgf M(s) such that M(α) < ∞ and M(s) = ∞ for s > α.9
Then the mgf of WT is infinite for s > α and satisfies
E[esWT ] =

∞
X

p(1 − p)k−1 M(s)k =

k=1

pM(s)
<∞
1 − (1 − p)M(s)

for s ∈ [0, α] if p ∈ (0, 1) is sufficiently close to 1 because M(α) < ∞. Thus the
right abscissa of convergence of WT is α and recalling (2.2) we have
lim sup
w→∞

1
log P(WT > w) = −α,
w

but there is no s > 0 that satisfies ρ(ΠD(s)) = M(s) = 1/(1 − p), so we may not
appeal to Theorem 3.1 to strengthen the above limit superior to an ordinary limit.
Remark 3.4. A sufficient condition for there to exist α > 0 in the interior of I


such that ρ(ΠD(α)) = 1/(1 − p) is that (i) E esX1 J1 = n < ∞ for all s > 0

and all n ∈ N and (ii) πnn > 0 and P(X1 > 0 | J1 = n) > 0 for some n ∈ N .
To see this, note that under condition (i) the matrix D(s) is well-defined for all

if a ≥ b and a 6= b, and that a ≫ b if an > bn for all n.
9
One such example is providedRby the density f (x) = 1(x ≥ 1)Cx−κ−1 e−αx , where κ, α > 0
and C > 0 is a number such that f (x) dx = 1. The corresponding mgf is finite for s ≤ α and
infinite for s > α.

18



s. Let M(s) be a matrix whose n-th diagonal entry is πnn E esX1 J1 = n for one
n satisfying condition (ii) and 0 for all other entries. Clearly ΠD(s) ≥ M(s). By
Proposition B.2(ii), we have


∞ > ρ(ΠD(s)) ≥ ρ(M(s)) = πnn E esX1 J1 = n → ∞

as s → ∞ by condition (ii), so ρ(ΠD(s)) crosses 1/(1 − p).10 A similar sufficient
condition applies to the lower tail exponent.
Remark 3.5. Theorem 3.1 is somewhat related to the (asymmetric) Laplace distribution, which is the logarithm of the double Pareto distribution. When {Xt }∞
t=1
is iid Laplace and T is a geometric random variable with mean 1/p, then the gePT
ometric sum WT =
t=1 Xt is also Laplace (Kotz et al., 2001, p. 151). When

{Xt }∞
t=1 is iid with a general distribution with finite variance, then WT (properly
scaled) weakly converges to the Laplace distribution as p → 0 (Kotz et al., 2001,
pp. 152–155). Toda (2014, Theorem 15) proves the same for the non-iid case
provided that the central limit theorem holds for {Xt }∞
t=1 .

Example 3.3 (Gaussian distribution). Consider the geometric Brownian motion
(1.3). Using Itô’s lemma, we obtain
d log St =



σ2
µ−
2



dt + σ dBt .

The discrete-time analog of this process is
Wt = Wt−1 + Xt ,

Xt ∼ iidN((µ − σ 2 /2)∆, σ 2 ∆),

where the choice of ∆ > 0 depends on the unit of time. Suppose that in the
continuous-time model, the age distribution is exponential with parameter η. Then
in the discrete-time analog, the birth/death probability is p = 1−e−η∆ . By Remark
3.2, WT has exponential tails with exponents ζ = α, −β solving the equation
E[eζX1 ] = 1/(1 − p). But with X1 ∼ N((µ − σ 2 /2)∆, σ 2 ∆) this equation can be
rewritten as
e(µ−σ

2 /2)∆ζ+σ 2 ∆ζ 2 /2

= eη∆ ,

which is equivalent to (1.4). Therefore with iid Gaussian shocks, the tail exponents
10

There may not be α such that ρ(ΠD(α)) = 1/(1 − p) if condition (ii) fails. A counterexample
is the two-state Markov chain that takes the values ±1 with transition probability matrix Π =
[(0, 1)⊤ , (1, 0)⊤ ]. Then ΠD(s) = [(0, e−s )⊤ , (es , 0)⊤ ], whose spectral radius is identically equal
to 1. In this case WT ∈ {−1, 0, 1} (bounded), so it does not have exponential tails.

19

are identical in the discrete-time and continuous-time models.
Example 3.4 (Two-state Markov chain). Suppose that {Xt }∞
t=1 is a two-state
time-homogeneous Markov chain, as in Example 3.2 with N = 2. This is the
simplest example of a non-iid process satisfying the assumptions of Theorem 3.1.
The cumulated process Wt is the “Markov trend in levels” studied by Hamilton
(1989). By the Perron-Frobenius theorem, ρ(ΠD(s)) is the maximum eigenvalue
of ΠD(s), and so we compute

ρ(ΠD(s)) = ρ

#!
"
π11 esx1 π12 esx2

π21 esx1 π22 esx2


q
1
2
sx1
sx2
sx
sx
s(x
+x
)
=
π11 e + π22 e + (π11 e 1 − π22 e 2 ) + 4π12 π21 e 1 2 . (3.11)
2

Setting this quantity equal to 1/(1 − p) and solving for s gives a unique positive
solution α provided that πnn > 0 and xn > 0 for some n ∈ {1, 2}, and a unique
negative solution −β provided that πnn > 0 and xn < 0 for some n (Remark 3.4).

3.2

Refinements

In the proof of Theorem 3.1 we applied Theorem 2.1 in conjunction with Remark
2.3 to establish the limits (3.3) and (3.4). By exploiting the sharp bounds in (2.4),
we can improve Theorem 3.1 as follows. (A similar statement holds for the lower
tail, which we omit.)
Theorem 3.2. Let everything be as in Theorem 3.1, and let x, y be the right and
left Perron vectors of ΠD(α) whose entries sum to 1. Let z = s + it and suppose
that for some b > 0 the matrix pencil A(z) = I − (1 − p)ΠD(z) is invertible on
the axis of convergence Re z = α for t ∈ (−b, b) except at t = 0. Then
2πC/b
2πC/b
≤ lim inf eαw P(WT > w) ≤ lim sup eαw P(WT > w) ≤
,
w→∞
−1
1 − e−2πα/b
w→∞
(3.12)
where
pω0⊤ x
.
(3.13)
C=
(1 − p)2 y ⊤ ΠD ′ (α)x
e2πα/b

In particular, if z = α is the unique point of noninvertibility of A(z) on the axis
of convergence Re z = α, then
lim eαw P(WT > w) =

w→∞

20

C
.
α

(3.14)

Proof. Let Iα,b = {z = α + it : |t| < b}. We showed in the proof of Theorem 3.1
that A(z) is holomorphic on the interior of S with an isolated point of noninvertibility at z = α. If also A(z) is invertible on Iα,b \ {α} then we may deduce that
A(z) is holomorphic on an open set containing Iα,b with a unique point of noninvertibility at z = α. Therefore, Proposition C.1 implies that A(z)−1 and thus
pω0⊤ A(z)−1 ΠD(z)e are holomorphic on an open set containing Iα,b \ {α}. In view
of (3.7) the latter function constitutes a continuous extension of the mgf of WT
to the union of its strip of holomorphicity and Iα,b \ {α}; moreover, it was shown
in the proof of Theorem 3.1 that the singularity at α is a simple pole with residue
−C. By applying Theorem 2.1 with A = C (note that the lower and upper bounds
in (2.4) are increasing and decreasing respectively in B, so that if they are valid
for B > b then they are valid for b) we obtain (3.12).
If z = α is the unique point of noninvertibility of A(z) on the axis Re z = α,
then we can take b arbitrarily large. Letting b → ∞, both sides of (3.12) converge
to C/α and we obtain (3.14).
The following theorem characterizes the upper tail behavior of a geometrically
stopped random growth process. It is more generally applicable than (a discretetime reformulation of) the main result of Reed (2001) because the growth rate
process is permitted to be non-Gaussian and dependent on a Markov state variable;
on the other hand, we only characterize the upper tail of the stopped process, not
its entire distribution. A similar statement holds for the lower tail, which we omit.
Theorem 3.3. Let everything be as in Theorem 3.1. Let S0 > 0 be a random
variable independent of WT satifying E[S0α+ǫ ] < ∞ for some ǫ > 0, and define the
random variable S = S0 eWT . Then
log P(S > s)
= −α.
s→∞
log s
lim

(3.15)

If in addition z = α is the unique point of noninvertibility of A(z) on the axis of
convergence Re z = α, then
lim sα P(S > s) =

s→∞

C
E[S0α ],
α

(3.16)

where C is defined as in Theorem 3.2.
Proof. Since S0 is independent of WT , the mgf of log S is the product of the mgfs
of log S0 and WT . The moment condition E[S0α+ǫ ] < ∞ ensures that the mgf
of log S0 is holomorphic on the strip 0 < Re z < α + ǫ, and clearly it cannot
21

have a zero at α. The proof of Theorem 3.1 establishes that the mgf of WT
is holomorphic on the strip 0 < Re z < α with a simple pole at α. Therefore
the mgf of log S is also holomorphic on the strip 0 < Re z < α with a simple
pole at α, and applying Theorem 2.1 in conjunction with Remark 2.3 we find
that limx→∞ x−1 log P(log S > x) = −α. The limit (3.15) follows by substituting
x = log s.
If z = α is the unique point of noninvertibility of A(z) on Re z = α then from
Theorem 3.2 we have lims→∞ sα P(eWT > s) = C/α. The upper tail of the random
variable eWT is therefore regularly varying with index −α, and so by applying the
well-known lemma of Breiman (1965) we obtain
α

lim s P(S > s) =

s→∞



α

WT

lim s P(e

s→∞

> s)



P(S0 eWT > s)
lim
s→∞ P(eWT > s)



=

C
E[S0α ],
α

as claimed.
Theorem 3.3 provides a stronger characterization of the tail behavior of S when
z = α is the unique point of noninvertibility of A(z) on the axis of convergence
Re z = α. The following theorem shows that this is the case except when the
support of X1 is a subset of an evenly-spaced grid. For a scalar c, let cZ =
{cm : m ∈ Z}, where Z is the set of integers.
Theorem 3.4. Let everything be as in Theorem 3.3. If A(α + iτ ) is noninvertible
for some τ 6= 0, then there exist a1 , . . . , aN ∈ R such that
supp(X1 | J1 = n) ⊂ an +

2π
Z
τ

for all n ∈ N , with an = 0 if πnn > 0. Conversely, if
supp(X1 | J1 = n) ⊂

2π
Z
τ

for some τ 6= 0 and all n ∈ N , then A(α + iτ ) is noninvertible.
Proof. We divide our proof into three steps.
Step 1. Let X be a random variable and suppose that E[e(α+iτ )X ] = E[eαX ] < ∞
for some τ > 0. Then supp X ⊂ a + (2π/τ )Z for some a ∈ R. If E[e(α+iτ )X ] =
E[eαX ], then we can take a = 0.
Using the triangle inequality, we obtain


E[eαX ] = E[e(α+iτ )X ] ≤ E e(α+iτ )X = E[eαX ].
22

Fix any support point a ∈ supp X. Since the triangle inequality holds with equality, it must be the case that 0 ≤ e(α+iτ )x /e(α+iτ )a = e(α+iτ )(x−a) for all x ∈ supp X.
Therefore, for each x ∈ supp X, there exists m ∈ Z such that τ (x − a) = 2πm, so
supp X ⊂ a + (2π/τ )Z.
If E[e(α+iτ )X ] = E[eαX ], letting pm = P(X = a + 2πm/τ ) we obtain
X

pm eα(a+2πm/τ ) = E[eαX ] = E[e(α+iτ )X ]

m∈Z

=

X

pm e(α+iτ )(a+2πm/τ ) = eiτ a

m∈Z

X

pm eα(a+2πm/τ ) .

m∈Z

Therefore eiτ a = 1, so there exists m0 ∈ Z such that a = 2πm0 /τ . We thus obtain
supp X ⊂ a +

2π
2π
2π
Z=
(m0 + Z) =
Z.
τ
τ
τ

Step 2. If A(α + iτ ) is noninvertible, then there exist a1 , . . . , aN ∈ R such that
supp(X1 | J1 = n) ⊂ an + (2π/τ )Z for all n ∈ N , with an = 0 if πnn > 0.
Without loss of generality we may assume τ > 0. If A(α + iτ ) is noninvertible then 1/(1 − p) is an eigenvalue of ΠD(α + iτ ). Therefore 1/(1 − p) ≤
ρ(D(α + iτ )). Since D(α + iτ ) is a diagonal matrix whose n-th diagonal entry


is E e(α+iτ )X1 J1 = n , by the triangle inequality we obtain |D(α + iτ )| ≤ D(α)
entry-wise. Since Π is a nonnegative matrix, by Proposition B.2(ii) we obtain
1
1
≤ ρ(ΠD(α + iτ )) ≤ ρ(Π |D(α + iτ )|) ≤ ρ(ΠD(α)) =
.
1−p
1−p
Since all inequalities hold with equality and ΠD(α) is irreducible, and noting
that 1/(1 − p) is an eigenvalue of ΠD(α + iτ ), by Proposition B.2(iii) we have
ΠD(α + iτ ) = ΘΠD(α)Θ−1 , where Θ = diag(eiθ1 , . . . , eiθN ) for some θ1 , . . . , θN ∈
R. Comparing the (m, n)-th entry, we have




πmn E e(α+iτ )X1 J1 = n = πmn ei(θm −θn ) E eαX1 J1 = n

(3.17)

for all m, n ∈ N .

Since Π is irreducible, for each n ∈ N there exists m ∈ N such that πmn >
0. Taking the absolute value of (3.17) and dividing by πmn > 0, we obtain




E e(α+iτ )X1 J1 = n = E eαX1 J1 = n for all n ∈ N . It now follows from
the previous step that there exist a1 , . . . , aN ∈ R such that supp(X1 | J1 = n) ⊂
an + (2π/τ )Z for all n ∈ N .
23

If πnn > 0, setting m = n in (3.17) and dividing by πnn > 0, we obtain
 (α+iτ )X



1
E e
J1 = n = E eαX1 J1 = n . Hence by the previous step we have
supp(X1 | J1 = n) ⊂ (2π/τ )Z, so we can take an = 0.
Step 3. If supp(X1 | J1 = n) ⊂ (2π/τ )Z for all n ∈ N , then A(α + iτ ) is
noninvertible.
Fix any n ∈ N and let pm = P(X1 = 2πm/τ | J1 = n) for m ∈ Z. Then
X

 X


E e(α+iτ )X1 J1 = n =
pm e(α+iτ )2πm/τ =
pm e2παm/τ = E eαX1 J1 = n .
m∈Z

m∈Z

Since n is arbitrary, we have D(α + iτ ) = D(α), and consequently A(α + iτ ) =
A(α). We know that A(z) is noninvertible at α, so it must also be noninvertible
at α + iτ .
Remark 3.6. The condition supp(X1 | J1 = n) ⊂ (2π/τ )Z for all n ∈ N is
sufficient but not necessary for A(α + iτ ) to be noninvertible. To see this, note
that if supp(X1 | J1 = n) ⊂ an + (2π/τ )Z for all n, then (3.17) (which is sufficient
for noninvertibility) is equivalent to πmn eiτ an = πmn ei(θm −θn ) . This equation holds,
for example, if N = 2, π11 = π22 = 0, π12 = π21 = 1, τ = 2π, θ1 = 1, θ2 = −1,
a1 = −1/π, and a2 = 1/π. Then supp(X1 | J1 = n) ⊂ ±1/π + Z.
The following example shows that if A(z) is noninvertible at multiple points on
the axis of convergence Re z = α, the upper tail of S is not necessarily Paretian.
Example 3.5 (Deterministic growth). Let Xt = µ > 0, a constant. Wold and Whittle
(1957) used a continuous-time model with deterministic growth and mortality
rates to investigate the tail behavior of wealth distributions. In this case we have
A(z) = 1 − (1 − p)eµz , and the mgf of WT = µT is given by
zWT

E[e

]=

∞
X

p(1 − p)

t−1 µzt

t=1

e

peµz
=
.
1 − (1 − p)eµz

(3.18)

The right abscissa of convergence is α = − log(1 − p)/µ > 0. Setting the denominator in (3.18) equal to zero, we obtain poles at z = α + 2πik/µ for k ∈ Z.
These poles are the points of noninvertibility of A(z) on the axis of convergence
Re z = α. Let S = eWT . If t ≤ w/µ < t + 1, then
P(WT > w) = P(T > w/µ) = (1 − p)t = (1 − p)⌊w/µ⌋ .

24

Therefore
P(S > s) = P(log S > log s) = P(WT > log s) = (1 − p)⌊(log s)/µ⌋ .
Clearly log P(S > s)/ log s → log(1 − p)/µ = −α as s → ∞, consistent with
Theorem 3.3. However,
sα P(S > s) = s− log(1−p)/µ (1 − p)⌊(log s)/µ⌋ = (1 − p)⌊(log s)/µ⌋−(log s)/µ
oscillates between 1 and 1/(1 − p) as s → ∞, so the tail is not Paretian. These
limits are precisely the bounds we obtain in (3.12) by setting C = p/(µ(1 − p))
and b = 2π/µ.

3.3

Comparative statics

Theorems 3.1–3.3 show that geometrically stopped random growth processes have
Pareto tails under empirically plausible conditions. Our characterization of the
Pareto exponents is implicit, in the sense that they are given by the positive and
negative solutions to (3.2). It may be desirable to provide an explicit characterization of how the exponents vary when we perturb the exogenous parameters Π,
D, and p. In this section, we derive such comparative statics.
We consider perturbations that involve the stopping probability p, locationscale transformations of the increments Xt , and the overall persistence of the
hidden Markov state Jt . We interpret increases (decreases) in p as decreases
(increases) in lifespan. Regarding the location-scale transformations of Xt , let
Y1 , . . . , YN be zero mean random variables and suppose that the distribution of Xt
conditional on Jt = n is parametrized as µn + σn Yn , where µn ∈ R and σn > 0. Let
µ = (µ1 , . . . , µN )⊤ and σ = (σ1 , . . . , σN )⊤ , and let D(s; µ, σ) be the N ×N diagonal
matrix with n-th diagonal entry equal to E[es(µn +σn Yn ) ]. We interpret increases
(decreases) in µn and σn as increases (decreases) in growth and volatility in state
n. Regarding persistence, given our fixed irreducible transition probability matrix
Π we define Π(τ ) = τ I + (1 − τ )Π for τ ∈ [0, 1]. We interpret increases (decreases)
in τ as increases (decreases) in persistence. Note that since Π is irreducible, so is
Π(τ ) for τ ∈ (0, 1).
The following theorem shows that increasing the lifespan, growth, volatility,
or persistence makes the upper Pareto exponent smaller (the upper tail heavier),
which is intuitive. Its proof is an application of the implicit function theorem.
Analogous results hold for the lower tail.
25

Theorem 3.5. Let θ = (p, µ, σ, τ ) ∈ (0, 1) × RN × RN
++ × (0, 1) be the parameters
we perturb. Suppose that the conditions of Theorem 3.1 are satisfied for some
parameter value θ0 . Then there exists an open neighborhood Θ of θ0 and a unique
continuously differentiable function α : Θ → (0, ∞) such that
ρ(Π(τ )D(α(θ); µ, σ)) =

1
.
1−p

(3.19)

The partial derivatives of α satisfy the following inequalities at θ0 .
(i) ∂α/∂p > 0: Longer lifespan implies a smaller Pareto exponent.
(ii) ∂α/∂µn ≤ 0: Higher growth implies a smaller Pareto exponent.
(iii) ∂α/∂σn ≤ 0: Higher volatility implies a smaller Pareto exponent.
(iv) ∂α/∂τ ≤ 0: Higher persistence implies a smaller Pareto exponent.
Proof. Define F : (0, 1) × RN × RN
++ × (0, 1) × (0, ∞) → R ∪ {∞} by
F (θ, s) = ρ(Π(τ )D(s; µ, σ)) −

1
,
1−p

where the spectral radius is understood to be ∞ whenever one or more elements of
D(s; µ, σ) are infinite. Since the conditions of Theorem 3.1 are satisfied at the parameter value θ0 = (p0 , µ0, σ 0 , τ 0 ), there exists α0 > 0 in the interior of I such that
F (θ0 , α0 ) = 0. Let us show that there exists an open neighborhood U of (θ0 , α0 )
such that F is finite on U. Clearly F is finite if and only if E[es(µn +σn Yn ) ] < ∞ for
all n. Since E[es(µn +σn Yn ) ] = esµn E[esσn Yn ], it suffices to show that E[esσn Yn ] < ∞
if (s, σn ) is sufficiently close to (α0 , σn0 ). However, this follows trivially from the
definition of I in (3.1) and the fact that α0 lies in the interior of I.
To prove our claim we apply the implicit function theorem. Proposition B.3
tells us that the spectral radius of a nonnegative irreducible square matrix varies
holomorphically with local perturbations to the elements of that matrix. Therefore, noting that Π(τ 0 )D(α0; µ0 , σ 0 ) is nonnegative and irreducible and that the
elements of D(s; µ, σ) = diag((E[es(µn +σn Yn ) ])N
n=1 ) and Π(τ ) are continuously differentiable with respect to µ, σ, τ, s, we find that F is continuously differentiable
on U. In Step 1 of the proof of Theorem 3.1 we established that ρ(ΠD(s)) is
convex in s. It follows that F (θ0 , s) is convex in s, and so
0>−

1
∂F 0 0
p
=1−
= F (θ0 , 0) − F (θ0 , α0 ) ≥
(θ , α )(0 − α0 ),
1−p
1−p
∂s
26

implying that ∂F/∂s > 0 on a neighborhood of (θ0 , α0 ). We may therefore apply
the implicit function theorem, which guarantees the existence of an open neighborhood Θ of θ0 and a unique continuously differentiable function α : Θ → (0, ∞)
such that (3.19) holds. The partial derivatives of α on Θ are then given by
∇α(θ) = −

1
∇θ F.
∂F/∂s

It remains only to show that the partial derivatives of F with respect to
p, µn , σn , τ have the appropriate sign at (θ0 , α0). We consider them in turn.
(i) Noting that ∂F/∂p = −1/(1 − p)2 < 0, we deduce that ∂α/∂p > 0 at θ0 .
(ii) F (θ, s) depends on µn only through the n-th diagonal element of D(s; µ, σ),
which we denote by Mn (s; µn , σn ) = E[es(µn +σn Yn ) ]. Since for s > 0 we have
∂Mn
= s E[es(µn +σn Yn ) ] > 0,
∂µn
Mn is increasing in µn . Since by Proposition B.2(ii) the spectral radius of a
nonnegative matrix is increasing in its elements, it follows that ∂F/∂µn ≥ 0
at (θ0 , α0). Therefore ∂α/∂µn ≤ 0 at θ0 .
(iii) Since by assumption Yn has zero mean,
∂Mn
∂σn

= s E[es(µn +σn Yn ) Yn ]
σn =0

σn =0

= 0.

Since the exponential function is convex, Mn is convex in σn . Therefore
∂Mn /∂σn ≥ 0 for s, σn > 0. By the same argument as the previous case, we
have ∂F/∂σn ≥ 0 at (θ0 , α0 ) and hence ∂α/∂σn ≤ 0 at θ0 .
(iv) Proposition B.4 guarantees that
∂
∂F 0 0
(θ , α ) =
ρ((τ I + (1 − τ )Π)D(α0 ))
∂τ
∂τ

≥ 0,
τ =τ 0

implying that ∂α/∂τ ≤ 0 at θ0 .
Example 3.6 (Two-state Markov chain). Consider the random growth model
where the growth rate process {Xt }∞
t=1 is a two-state Markov chain as in Example
3.4. Let Xt = µn > 0 if Jt = n, where n = 1, 2. One can interpret this model
as the deterministic growth model of Wold and Whittle (1957) except that the
growth rate has two regimes.
27

Let η = − log(1 − p) be the death rate. Suppose that the transition probability
matrix is
"
# "
#
0 1
τ
1−τ
Π(τ ) = τ I + (1 − τ )
=
1 0
1−τ
τ
and let α(τ ) be the Pareto exponent corresponding to τ . By Theorem 3.5, α(τ )
is decreasing in τ . Using (3.11), we obtain
α(0) =

η
µ̄

and

lim α(τ ) = min
n

τ ↑1

η
,
µn

where µ̄ = (µ1 + µ2 )/2.
As a numerical example, suppose that (µ1 , µ2 ) = (0.03, 0.01) and η = 0.04,
so units grow by 3% and 1% in each state and the average life span is about 25
periods. The blue curve in Figure 3.4 shows the Pareto exponent. In this case
µ̄ = 0.02, so α(0) = η/µ̄ = 2 and α(1) = η/µ1 = 4/3. The red line corresponds to
the case where {Xt }∞
t=1 is iid with Xt = µn with probability 1/2. We can see that
the two Pareto exponents are similar when the persistence τ of the Markov chain
is not so high, say τ ∈ [0, 0.8]. However, increasing τ beyond 0.8 quickly makes
the Pareto exponent small. In the limit as τ → 1, the exponent converges to the
case of highest growth maxn µn .
2.2
Markov
iid

Pareto exponent

2

1.8

1.6

1.4

1.2
0

0.2

0.4

0.6

0.8

1

persistence

Figure 3.4: Pareto exponent of the two-state random growth model.

28

4

Tails of wealth in heterogeneous-agent models

As an application of Theorems 3.1–3.3, in this section we characterize the tails of
the wealth distribution in heterogeneous-agent models with non-Gaussian, Markovian shocks. We consider two cases, models with idiosyncratic endowment and
investment risks. All omitted proofs in this section may be found in Appendix D.

4.1

Wealth distribution in CARA Huggett economies

The first application is to a standard Huggett (1993) economy, where agents are
subject to uninsurable endowment risk and trade a risk-free asset. We present
an analytical solution to this model with Markovian shocks that exploits constant
absolute risk aversion (CARA) preferences. We first consider the single agent
problem (partial equilibrium), which we subsequently embed into a general equilibrium model.
4.1.1

Single agent problem

Consider an agent with additive CARA utility
E0

∞
X

β t u(ct ),

(4.1)

t=0

where ct ∈ R is consumption at period t, the discount factor is β ∈ (0, 1), and
the period utility function u(c) = −e−γc /γ exhibits constant absolute risk aversion
γ > 0.11 The agent can borrow or save at an exogenous gross risk-free rate R > 1
and is subject to uninsurable idiosyncratic income risk.12 Letting yt be the income
at t and wt the financial wealth at the beginning of period t excluding the current
income yt , the budget constraint of the agent is
wt+1 = R(wt − ct + yt ).
11

(4.2)

We focus on CARA preferences because they are tractable with additive shocks. For example, the income process is iid Gaussian in Calvet (2001), non-Gaussian AR(1) in Wang (2003),
and non-Gaussian VAR(1) in Toda (2017a). Also, since the CARA utility is defined on the entire
real line, as is common with these models we assume that consumption can be negative.
12
The assumption R > 1 is without loss of generality: since the expected present value of
future income is infinite when R ≤ 1, without borrowing constraints (as in the current setup)
the income fluctuation problem does not have a solution. In the subsequent general equilibrium
analysis, we prove that a stationary equilibrium with R > 1 always exists.

29

Given the initial wealth w0 and the income process {yt }∞
t=0 , the agent’s objective
is to maximize the utility (4.1) subject to the budget constraint (4.2) for all t ≥ 0.
Suppose that the income process {yt }∞
t=0 is a time-homogeneous finite-state
Markov chain. Let s ∈ S = {1, . . . , S} denote the states, P = (pss′ ) be the
(irreducible) transition probability matrix, and y = (ȳs )Ss=1 be the vector of incomes. By merging redundant states, without loss of generality we may assume
ȳ1 < · · · < ȳS . Due to time-homogeneity, the state variables of the utility maximization problem are the current wealth w = wt and the exogenous state s = st ,
so let Vs (w) be the value function in state s. The Bellman equation is
Vs (w) = max {u(c) + β E [Vs′ (R(w − c + ȳs )) | s]} ,
c

(4.3)

where s′ is next period’s state.
The following proposition characterizes the optimal consumption rule in closedform up to a system of S nonlinear equations in S unknowns.
Proposition 4.1. Suppose that R > 1. Then the utility maximization problem
has a unique solution. The value function and optimal consumption rule are
Vs (w) = −

1 −γ(aw+bs )
e
,
γa

cs (w) = aw + bs ,

(4.4a)
(4.4b)

where a = 1 − 1/R and b = (bs )Ss=1 uniquely solves
b=
4.1.2



1
1−
R



y−


1
log βRP e−γb .13
γR

(4.5)

General equilibrium and wealth distribution

Next we embed the above single agent problem into the general equilibrium. Consider a Huggett (1993) economy, where there are infinitely many independent
agents receiving income and trading a risk-free asset in zero net supply. We restrict the analysis to a stationary equilibrium, where the risk-free rate is constant
over time.
Before formally stating the model, we note that we cannot obtain stationarity
just by replicating agents, and therefore we need to introduce additional assumptions. To see this, use the budget constraint (4.2), the consumption rule (4.4b),
For notational simplicity, we apply functions entry-wise. For example, for b = (bs )Ss=1 , we
have e−γb = (e−γbs )Ss=1 .
13

30

and the equation a = 1 − 1/R to write
wt+1 = R(wt − (awt + bst ) + ȳst ) = wt + R(ȳst − bst ).

(4.6)

Since wealth is a random walk in levels, with infinitely lived agents there is no
stationary wealth distribution. In order to obtain a stationary distribution, we
assume that agents enter/exit the economy at constant probability p as in Yaari
(1965) and Blanchard (1985). Because agents survive each period with probability
1 − p, the effective discount factor is β̃ = β(1 − p).
The formal model works as follows. There are countably infinite agents indexed
by i ∈ N. Agent i’s consumption, income, and wealth at time t are denoted by
cit , yit , and wit . All agents have the same utility function (4.1). When an agent
is born, their initial wealth is 0 and their income is ȳs , where the initial state s
is drawn from the stationary distribution π = (πs )Ss=1 of the Markov chain with
transition probability matrix P = (pss′ ) independently across agents. By the
strong law of large numbers, the limits of average quantities such as
I

1X
cit ,
I i=1

I

1X
yit ,
I i=1

I

1X
wit
I i=1

exist almost surely as I → ∞. We refer to such quantities as (per capita) aggregate
consumption, income, and wealth and denote by capital letters Ct , Yt , Wt , etc.
Because agents die with probability p each period, we need to specify what
happens to the risk-free asset position of deceased agents. Suppose that there are
perfectly competitive insurance companies that offer annuities and life insurances.
Let R be the gross risk-free rate and R̃ be the “effective” risk-free rate that agents
face. If an agent saves or borrows 1, the position grows to R̃ next period if the
agent survives, and 0 if he dies (an agent who dies with positive assets surrender
to the insurance company; the debt of an agent who dies with negative assets
is covered by life insurance). Letting At be aggregate savings, by accounting we
obtain RAt = (1 − p)R̃At + p0, and thus deduce that R̃ = R/(1 − p).
Now we can formally define the stationary equilibrium.
Definition 4.1 (Stationary equilibrium of CARA Huggett economy). A stationary equilibrium of a CARA Huggett economy consists of a gross risk-free rate
R > 0 and a sequence of consumption and wealth {(cit , wit )i∈N }∞
t=0 such that
(i) (Optimization) for all i ∈ N, {(cit , wit )}∞
t=0 maximizes utility (4.1) subject
to the budget constraint (4.2), where β, R are replaced with β̃ = β(1 − p)
31

and R̃ = R/(1 − p),
(ii) (Market clearing) the risk-free asset market clears, i.e., the aggregate savings
is 0, and
(iii) (Stationarity) all aggregate quantities are constant over time almost surely.
To derive the market clearing condition, let C, W be aggregate consumption
and wealth in the stationary equilibrium. Since there is no aggregate savings in
a Huggett economy, we have W = 0. Aggregating the optimal consumption rule
(4.4b) across agents, the aggregate consumption is
C = aW + π ⊤ b = π ⊤ b,
where π = (πs )Ss=1 is the stationary distribution of the Markov chain. Furthermore, aggregate consumption must equal aggregate income, so C = Y = π ⊤ y.
Combining these two equations, the equilibrium condition is
π ⊤ b = π ⊤ y.

(4.7)

The following theorem shows the existence of a stationary equilibrium.
Theorem 4.2. There exists a stationary equilibrium in the CARA Huggett economy. The gross risk-free rate satisfies 1 − p < R ≤ 1/β.
The following theorem shows that the stationary wealth distribution in a
CARA Huggett economy has exponential tails.
Theorem 4.3. Let P = (pss′ ) be the transition probability matrix with positive
diagonal entries, y = (ȳ1 , . . . , ȳS )⊤ be the vector of incomes, b = (b1 , . . . , bS )⊤ be
R
be the effective risk-free rate in equilibrium.
as in Proposition 4.1, and R̃ = 1−p
For any z ∈ R, let λ(z) > 0 be the Perron root of the matrix P (diag ez R̃(y−b) ). Then
1
. The stationary
there exist unique α1 , α2 > 0 such that λ(α1 ) = λ(−α2 ) = 1−p
wealth distribution has exponential tails with upper tail exponent α1 and lower tail
exponent α2 .

Proof. By (4.6) with R replaced by R̃, individual wealth increases by R̃(ȳs − bs ) in
state s. Since P has positive diagonal entries and the vector of moment generating
functions is M(z) = ez R̃(y−b) , assuming y − b has a positive element, by Theorem
3.1 (see Remark 3.4) the wealth distribution has an exponential upper tail with
1
, where λ(z) is the Perron root of
exponent α1 > 0 that satisfies λ(α1 ) = 1−p
P (diag M(z)). The argument for the lower tail is similar.
32

Therefore to complete the proof it remains to show that y − b has both positive
and negative elements. Since by the equilibrium condition we have π ⊤ (y − b) = 0,
it suffices to show that y 6= b. Suppose on the contrary that y = b. Then by (4.5)
we obtain γy = − log(βRP e−γy ), or equivalently
Pv =

1
v
βR

with v = e−γy ≫ 0. Therefore v is an eigenvector of P with eigenvalue 1/(βR).
Since P is an irreducible nonnegative matrix with spectral radius 1, and βR ≤ 1
by Theorem 4.2, it must be the case that βR = 1 and P v = v. By the PerronFrobenius theorem, v must be a scalar multiple of the Perron vector, which is 1.
Therefore ȳ1 = · · · = ȳS , which contradicts the assumption ȳ1 < · · · < ȳS .
In the quantitative macro literature, it is well known that heterogeneous-agent
models with idiosyncratic income risk alone have difficulty in matching the wealth
distribution (Huggett, 1996; Castañeda et al., 2003). Theorem 4.3 provides a theoretical explanation of these numerical findings: idiosyncratic income risk (additive
shocks) alone can only generate exponential tails—not Pareto tails—when the income process is light-tailed. Benhabib et al. (2017) show that when the income
process is fat-tailed, then the tail exponent of the wealth distribution is identical
to that of the income process, which is counterfactual. These theoretical results
suggest that researchers need to go beyond models with idiosyncratic income risk
in order to understand the upper tail of the wealth distribution.
4.1.3

Numerical example

As a numerical example, let the discount factor be β = 0.96, absolute risk aversion
γ = 3, death probability p = 0.025, and suppose that log income xt = log yt is a
centered stationary Gaussian AR(1) process,
xt = ρxt−1 + εt ,

εt ∼ N(0, σ 2 ),

where ρ = 0.9 and σ = 0.1. Since our theory concerns a finite-state Markov chain,
instead of the original AR(1) process we take a discrete-state analog as the true
process. For this purpose we apply the method in Farmer and Toda (2017) to
discretize the AR(1) process {xt } using an even-spaced grid with S = 9 points.
Letting x̄s be the grid point in state s, the income in state s is defined by ȳs = ex̄s .
We then solve for the equilibrium risk-free rate R = 1.0350 using (4.5) and (4.7).

33

Finally, we compute the tail exponents of the wealth distribution using Theorem
4.3. The theoretical values are α1 = 0.2857 and α2 = 0.3669.
Figure 4.1a shows the histogram of the wealth from a simulation with 100,000
agents. The wealth distribution seems to have exponential tails. Letting F (x) be
the empirical cumulative distribution function, Figure 4.1b shows the tail prob(x)
for x < 0, and 1−F
for x > 0) in a semi log plot. Since Figure
ability ( FF (x)
(0)
1−F (0)
4.1b shows straight line patterns, the tails are exponential. Estimating the tail
exponents by maximum likelihood (using the largest 5% of observations in each

tail), we obtain α
b1 = 0.2805 and α
b2 = 0.3657, close to the theoretical values.
10000

0

log tail probability

8000

6000

4000

2000

0
-20

-10

0

10

-2

-4

-6

-8
-20

20

wealth

-10

0

10

20

wealth

(a) Histogram.

(b) Log tail probability.

Figure 4.1: Wealth distribution of the CARA Huggett economy.

4.2

Wealth distribution with idiosyncratic investment risk

Heterogeneous-agent models with idiosyncratic investment risk have recently become quite popular for studying size distributions, but most papers assume iid
shocks.14 Benhabib et al. (2011) consider an overlapping generations model where
the returns on wealth and income remain constant over the life cycle of the agent
but change when the generation turns over. This way they are able to solve for the
optimal consumption rule during the lifetime but allow for Markovian dynamics
across generations using the results in da Saporta (2005). Toda (2014) considers
a heterogeneous-agent model with Markovian shocks but the Pareto tail results
hold only in the continuous-time limit, where shocks are Gaussian. Cao and Luo
(2017) consider a continuous-time Markov switching model with two states. Using our main results (Theorems 3.1–3.4), applied researchers may characterize the
14

See, for example, Luttmer (2007), Nirei and Souma (2007), Benhabib et al. (2015, 2016),
Acemoglu and Cao (2015), Toda and Walsh (2015), Gabaix et al. (2016), Nirei and Aoki (2016),
and Aoki and Nirei (2017).

34

tail behavior of the size distribution provided that the agents solve a homogeneous problem (i.e., maximizing a homothetic function subject to proportional
constraints such as an optimal consumption-portfolio problem) in a Markovian
setting. Due to space considerations we do not discuss a concrete example, but
interested readers may refer to the above cited papers.

5

Power law in Japanese prefectures

As an empirical illustration, we apply the random growth model to a panel of
Japanese prefecture populations. The main question is whether the time series
properties of population dynamics estimated from the panel are consistent with a
Pareto index estimated from the cross-section.15

5.1

Cross-sectional estimation

We use annual population data for 46 Japanese prefectures (excluding Okinawa)
since 1920. Only for the cross-sectional estimation, we also use the census data in
1873, 1884, 1893, 1908, and 1913. See Appendix E.1 for a more detailed description
of the data.
Figure 5.1 displays estimates of the upper Pareto exponent of the size distribution of prefectures in different years. The blue curve is the maximum likelihood estimate of the Pareto-lognormal (PlN) distribution, which is the product
of independent lognormal and Pareto random variables (see Appendix E.2 for details). The red curve is obtained by running the log(Rank − 1/2) regression as in
Gabaix and Ibragimov (2011) using the largest 15 prefectures in each year.
We can make a few observations. First, the two curves in Figure 5.1 are similar
after 1920 except during WWII, which suggests the robustness of the estimation.
Since the maximum likelihood estimation of PlN using the whole sample is more
efficient than the log(Rank − 1/2) regression using the tail observations and does
not require a cutoff selection, below we focus on the former, although the latter
may be more robust to misspecification.
Second, the Pareto exponent estimated by maximum likelihood is very large
before 1920 and during WWII. Since PlN nests the lognormal distribution as a
limiting case (α → ∞), this suggests that the prefecture distribution was closer
to lognormal before 1920 and during WWII. Figure 5.2a shows the p-value of the
15

Davis and Weinstein (2002) examine historical population data from Stone Age Japanese
regions and reject the random growth hypothesis, but this is because they narrowly interpret
“random growth” as “iid growth”.

35

Pareto exponent

5

4

3

2
Pareto-lognormal
log(Rank-1/2) regression

1
1880 1900 1920 1940 1960 1980 2000

Year

Figure 5.1: Pareto exponent of Japanese prefectures.

likelihood ratio test of the lognormal distribution (α = ∞) against the Paretolognormal (α < ∞). The lognormal distribution is rejected at significance level
0.05 for the years 1933–1943 and all years since 1960. It is not surprising that
the power law breaks down during WWII since much of the urban population
migrated to rural areas in order to avoid the air raid by the United States.
1

0.3
0.25

0.8

p-value

p-value

0.2
0.6

0.4

0.15
0.1

0.2

0.05

0

0
1880 1900 1920 1940 1960 1980 2000

1880 1900 1920 1940 1960 1980 2000

Year

Year

(a) Test of lognormality (α = ∞).

(b) Test of Zipf’s law (α = 1).

Figure 5.2: p-values for likelihood ratio tests.

Third, except for the break during WWII, the Pareto exponent is roughly
monotonically decreasing and approaches one, which is consistent with Zipf’s law.
Figure 5.2b shows the p-value of the likelihood ratio test of Zipf’s law (α = 1)
against the unrestricted Pareto-lognormal. Zipf’s law fails to be rejected at significance level 0.05 for all years since 1968.
36

Figure 5.3 displays log-log plots of the rank size distribution of prefectures for
several years, together with the values predicted from the Pareto-lognormal and
lognormal distributions using the maximum likelihood estimates. In 1873, the
size distribution was approximately lognormal. The power law pattern starts to
emerge in 1920, and the Pareto exponent drifts downward toward one until 1944.
Dislocation and destruction caused by WWII halved the population of Tokyo in
1945, which made the size distribution close to lognormal again. The power law
pattern reemerges around 1950 and stabilizes by 1970.
10 2

1873

1920

10 2

Data
PlN
lognormal

Rank

Rank

Data
PlN
lognormal

10 1

10 0
10 5

10 6

10 1

10 0
10 5

10 7

Prefecture population
10 2

10 6

10 7

Prefecture population

1944

1945

10 2

Data
PlN
lognormal

Rank

Rank

Data
PlN
lognormal

10 1

10 0
10 5

10 6

10 1

10 0
10 5

10 7

Prefecture population
10 2

10 6

10 7

Prefecture population

1950

1970

10 2

Data
PlN
lognormal

Rank

Rank

Data
PlN
lognormal

10 1

10 0
10 5

10 6

10 1

10 0
10 5

10 7

Prefecture population

10 6

10 7

Prefecture population

Figure 5.3: Log-log plot of rank size distribution of prefectures.

37

10 8

5.2

Panel estimation

Are the time series properties of the population dynamics consistent with a power
law in the cross-section? To answer this question, we estimate a random growth
model for the population dynamics and compute the theoretical Pareto exponent.
Due to data quality concerns and the structural break during WWII, we only use
post war data (1946–2015). We assume that the relative size (population of a
prefecture divided by the total population) Sit of prefecture i in year t follows the
random growth process
Si,t+1 = Gi,t+1 Sit ,
where Gi,t+1 is the gross growth rate of the relative size of prefecture i from year t
to t+1. Following the setting in Section 3, suppose that there are N states indexed
by n = 1, . . . , N and in state n the growth rate is lognormally distributed, so
log Gi,t+1 | nit = n ∼ N(µn , σn2 ),
where nit is the state of prefecture i in year t. For the number of states, we consider
N = 1, 2, 3 and estimate the model using the Hamilton (1989) filter (see Appendix
E.3 for details). After estimating the transition probability matrix Π = (πnn′ ),
conditional mean µ = (µ1 , . . . , µN )⊤ , and standard deviation σ = (σ1 , . . . , σN )⊤ ,
we compute the implied Pareto exponent α > 0 by solving the equation
2 2 /2

ρ(Π diag(eµ1 s+σ1 s

2

, . . . , eµN s+σN s

2 /2

)) =

1
,
1−p

(5.1)

which corresponds to (3.2) in Theorem 3.1. Since some regions in Japan were
developed much earlier than others (for example, the region around Kyoto and
Osaka were well-developed by the 6th century, while Hokkaido was developed only
in the late 19th century), the tail in the age distribution is long. We therefore set
p = 0 when computing the upper tail exponent.16 Table 5.1 shows the results.
When the model is iid (N = 1), the implied Pareto exponent α = 56.7 is
wildly inconsistent with the value estimated from the cross-section (about 1.3–
1.5 in 2015 according to Figure 5.1). When we expand the model to two states
(N = 2), the log-likelihood increases dramatically, so the model with a single
state is clearly misspecified. The expected growth rate does not change much
(only from −0.0035 to −0.0030) but there are now two volatility regimes with
identical growth rate. The implied Pareto exponent α = 26.8 decreases by half
16

Empirically, we find that the spectral radius in (5.1) is equal to one at s = 0 and at a positive
value of s; we take the latter value to be the implied upper Pareto exponent.

38

Table 5.1: Estimation of random growth model for Japanese prefectures.

Parameter

Number of states (N)
2



0.9439
0.9754 0.0246

0.0145
1
0.0283 0.9717
0.0210

 
−0.0030
−0.0122
−0.0035 −0.0030



0.0029 0.0169
0.0053
0.0111
1

Π
µ⊤
σ⊤
log L
α

9,925
56.7

11,638
26.8

3


0.0561 0.0000
0.9671 0.0184
0.0141 0.9649 
−0.0022 0.0084

0.0026 0.0199
12,388
1.61

Note: The structural parameters are Π: transition probability matrix, µ: conditional mean,
σ: conditional standard deviation. log L is the log-likelihood. α is the upper Pareto exponent
obtained by solving (5.1) for s.

compared to the case N = 1 but is still very large. With three states (N = 3),
the log-likelihood increases dramatically again and the Pareto exponent α = 1.61
becomes similar to the cross-sectional estimate. Therefore the emergence of a
power law in the Japanese prefecture sizes seems to be consistent with a random
growth model, but only by allowing for Markovian dynamics. Note that in each
case the estimated transition probability matrix has large diagonal elements, so
the process is very persistent. This is another reason why our non-iid results in
Section 3 are useful.

6

Conclusion

It has long been conjectured that random growth models robustly generate power
law tails. For example, in his well-known paper on Zipf’s law, Gabaix (1999,
footnote 13) writes:
Let each city grow at an arbitrary mean rate. [. . . ] [I]t does not matter
if this mean rate is time varying[, which] is a conjecture that we firmly
believe to be true. [. . . ] However, we could not find any argument in
the mathematical literature—here we deal with Markov chains with
time-varying transition matrices—to help us establish this rigorously.
In this paper we have presented results that affirmatively resolve the robustness
conjecture of power law tails. It is not essential that the growth rate process is
Gaussian or iid. While the details of the underlying process matter quantitatively
39

in the sense that they affect the magnitude of the Pareto exponent, they do not
matter qualitatively: the tail behavior is Pareto under mild conditions.
Applying a random growth model to historical Japanese prefecture population
data, we have shown that the time series properties of the growth process and
the tail behavior of the cross-sectional distribution of levels can only be reconciled
by introducing a hidden Markov state variable with several regimes. Because
our main theorems provide an analytical tool to study such non-iid processes, we
believe they have potentially a wide range of applications.

References
Acemoglu, D. and D. Cao (2015): “Innovation by Entrants and Incumbents,”
Journal of Economic Theory, 157, 255–294.
Ahlfors, L. V. (1979): Complex Analysis, New York: McGraw-Hill, third ed.
Aiyagari, S. R. (1994): “Uninsured Idiosyncratic Risk and Aggregate Saving,”
Quarterly Journal of Economics, 109, 659–684.
Aoki, S. and M. Nirei (2017): “Zipf’s Law, Pareto’s Law, and the Evolution
of Top Incomes in the U.S.” American Economic Journal: Macroeconomics, 9,
36–71.
Arata, Y. (2015): “Essays on Stochastic Macroeconomics,” Ph.D. thesis, University of Tokyo.
Arkolakis, C. (2016): “A Unified Theory of Firm Selection and Growth,” Quarterly Journal of Economics, 131, 89–155.
Axtell, R. L. (2001): “Zipf Distribution of U.S. Firm Sizes,” Science, 293,
1818–1820.
Beare, B. K., J. Seo, and W.-K. Seo (2017): “Cointegrated Linear Processes
in Hilbert Space,” Journal of Time Series Analysis, 38, 1010–1027.
Benhabib, J. and A. Bisin (2017): “Skewed Wealth Distributions: Theory and
Empirics,” Journal of Economic Literature, forthcoming.
Benhabib, J., A. Bisin, and M. Luo (2017): “Earnings Inequality and Other
Determinants of Wealth Inequality,” American Economic Review Papers and
Proceedings, 107, 593–597.
Benhabib, J., A. Bisin, and S. Zhu (2011): “The Distribution of Wealth
and Fiscal Policy in Economies with Finitely Lived Agents,” Econometrica, 79,
123–157.

40

——— (2015): “The Wealth Distribution in Bewley Economies with Capital Income Risk,” Journal of Economic Theory, 159, 489–515.
——— (2016): “The Distribution of Wealth in the Blanchard-Yaari Model,”
Macroeconomic Dynamics, 20, 466–481.
Bewley, T. F. (1983): “A Difficulty with the Optimum Quantity of Money,”
Econometrica, 51, 1485–1504.
Blackwell, D. (1965): “Discounted Dynamic Programming,” Annals of Mathematical Statistics, 36, 226–235.
Blanchard, O. J. (1985): “Debt, Deficits, and Finite Horizons,” Journal of
Political Economy, 93, 223–247.
Breiman, L. (1965): “On Some Limit Theorems Similar to the Arc-Sin Law,”
Theory of Probability & Its Applications, 10, 323–331.
Calvet, L.-E. (2001): “Incomplete Markets and Volatility,” Journal of Economic
Theory, 98, 295–338.
Cao, D. and W. Luo (2017): “Persistent Heterogeneous Returns and Top End
Wealth Inequality,” Review of Economic Dynamics, 26, 301–326.
Castañeda, A., J. Dı́az-Giménez, and J.-V. Rı́os-Rull (2003): “Accounting for the U.S. Earnings and Wealth Inequality,” Journal of Political Economy,
111, 818–857.
Champernowne, D. G. (1953): “A Model of Income Distribution,” Economic
Journal, 63, 318–351.
Coad, A. (2010): “The Exponential Age Distribution and the Pareto Firm Size
Distribution,” Journal of Industry, Competition and Trade, 10, 389–395.
da Saporta, B. (2005): “Tail of the Stationary Solution of the Stochastic Equation Yn+1 = an Yn + bn with Markovian Coefficients,” Stochastic Processes and
their Applications, 115, 1954–1978.
Davis, D. R. and D. E. Weinstein (2002): “Bones, Bombs, and Break Points:
The Geography of Economic Activity,” American Economic Review, 92, 1269–
1289.
Elsner, L., C. R. Johnson, and J. A. Dias da Silva (1988): “The Perron
Root of a Weighted Geometric Mean of Nonneagative Matrices,” Linear and
Multilinear Algebra, 24, 1–13.
Engle, R. F. and C. W. J. Granger (1987): “Co-integration and Error
Correction: Representation, Estimation and Testing,” Econometrica, 55, 251–
276.

41

Farmer, L. E. and A. A. Toda (2017): “Discretizing Nonlinear, Non-Gaussian
Markov Processes with Exact Conditional Moments,” Quantitative Economics,
8, 651–683.
Gabaix, X. (1999): “Zipf’s Law for Cities: An Explanation,” Quarterly Journal
of Economics, 114, 739–767.
——— (2009): “Power Laws in Economics and Finance,” Annual Review of Economics, 1, 255–293.
——— (2016): “Power Laws in Economics: An Introduction,” Journal of Economic Perspectives, 30, 185–206.
Gabaix, X. and R. Ibragimov (2011): “Rank−1/2: A Simple Way to Improve
the OLS Estimation of Tail Exponents,” Journal of Business and Economic
Statistics, 29, 24–39.
Gabaix, X., J.-M. Lasry, P.-L. Lions, and B. Moll (2016): “The Dynamics
of Inequality,” Econometrica, 84, 2071–2111.
Gibrat, R. (1931): Les Inégalités Économiques, Paris: Librairie du Recueil Sirey.
Giesen, K. and J. Suedekum (2014): “City Age and City Size,” European
Economic Review, 71, 193–208.
Giesen, K., A. Zimmermann, and J. Suedekum (2010): “The Size Distribution Across All Cities—Double Pareto Lognormal Strikes,” Journal of Urban
Economics, 68, 129–137.
Graham, S. W. and J. D. Vaaler (1981): “A Class of Extremal Functions for
the Fourier Transform,” Transactions of the American Mathematical Society,
265, 283–302.
Hamilton, J. D. (1989): “A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle,” Econometrica, 57, 357–384.
Horn, R. A. and C. R. Johnson (1985): Matrix Analysis, New York: Cambridge University Press.
Howland, J. S. (1971): “Simple Poles of Operator-Valued Functions,” Journal
of Mathematical Analysis and Applications, 36, 12–21.
Huggett, M. (1993): “The Risk-Free Rate in Heterogeneous-Agent IncompleteInsurance Economies,” Journal of Economic Dynamics and Control, 17, 953–
969.
——— (1996): “Wealth Distribution in Life-Cycle Economies,” Journal of Monetary Economics, 38, 469–494.
Johansen, S. (1991): “Estimation and Hypothesis Testing of Cointegration Vectors in Gaussian Vector Autoregressive Models,” Econometrica, 59, 1551–1580.
42

——— (1995): Likelihood-Based Inference in Cointegrated Vector Autoregressive
Models, Oxford, UK: Oxford University Press.
Kalecki, M. (1945): “On the Gibrat Distribution,” Econometrica, 13, 161–170.
Kamihigashi, T. (2014): “Elementary Results on Solutions to the Bellman Equation of Dynamic Programming: Existence, Uniqueness, and Convergence,” Economic Theory, 56, 251–273.
Kapteyn, J. C. (1903): Skew Frequency Curves in Biology and Statistics,
Groningen: P. Noordhoff.
Karlin, S. (1982): “Classifications of Selection-Migration Structures and Conditions for a Protected Polymorphism,” in Evolutionary Biology, ed. by M. K.
Hecht, B. Wallace, and G. T. Prance, New York: Plenum, vol. 14, chap. 2,
61–204.
Kesten, H. (1973): “Random Difference Equations and Renewal Theory for
Products of Random Matrices,” Acta Mathematica, 131, 207–248.
Klass, O. S., O. Biham, M. Levy, O. Malcai, and S. Solomon (2006):
“The Forbes 400 and the Pareto Wealth Distribution,” Economics Letters, 90,
290–295.
Korevaar, J. (2004): Tauberian Theory: A Century of Developments, no. 329
in Gruncllehren der mathematischen Wissenschaften, Berlin: Springer.
Kotz, S., T. J. Kozubowski, and K. Podgórski (2001): The Laplace Distribution and Generalizations, Boston: Birkhäuser.
Lévy, P. (1937): Théorie de l’Addition des Variables Aléatoires, Paris: GauthierVillars.
Lindeberg, J. W. (1922): “Eine Neue Herleitung des Exponentialgesetzes in der
Wahrscheinlichkeitsrechnung,” Mathematische Zeitschrift, 15, 211–225.
Lukacs, E. (1970): Characteristic Functions, London: Griffin, second ed.
Luttmer, E. G. J. (2007): “Selection, Growth, and the Size Distribution of
Firms,” Quarterly Journal of Economics, 122, 1103–1144.
Manrubia, S. C. and D. H. Zanette (1999): “Stochastic Multiplicative Processes with Reset Events,” Physical Review E, 59, 4945–4948.
Mimica, A. (2016): “Exponential Decay of Measures and Tauberian Theorems,”
Journal of Mathematical Analysis and Applications, 40, 266–285.
Mitzenmacher, M. (2004): “A Brief History of Generative Models for Power
Law and Lognormal Distributions,” Internet Mathematics, 1, 226–252.

43

Nakagawa, K. (2007): “Application of Tauberian Theorem to the Exponential
Decay of the Tail Probability of a Random Variable,” IEEE Transactions on
Information Theory, 53, 3239–3249.
Nirei, M. and S. Aoki (2016): “Pareto Distribution of Income in Neoclassical
Growth Models,” Review of Economic Dynamics, 20, 25–42.
Nirei, M. and W. Souma (2007): “A Two Factor Model of Income Distribution
Dynamics,” Review of Income and Wealth, 53, 440–459.
Pareto, V. (1896): Cours d’Économie Politique, Lausanne: F. Rouge.
Reed, W. J. (2001): “The Pareto, Zipf and Other Power Laws,” Economics
Letters, 74, 15–19.
——— (2002): “On the Rank-Size Distribution for Human Settlements,” Journal
of Regional Science, 42, 1–17.
——— (2003): “The Pareto Law of Incomes—An Explanation and an Extension,”
Physica A, 319, 469–486.
Reed, W. J. and B. D. Hughes (2002): “From Gene Families and Genera to
Incomes and Internet File Sizes: Why Power Laws Are So Common in Nature,”
Physical Review E, 66, 067103.
——— (2003): “On the Distribution of Family Names,” Physica A, 319, 579–590.
Reed, W. J. and M. Jorgensen (2004): “The Double Pareto-Lognormal
Distribution—A New Parametric Model for Size Distribution,” Communications in Statistics—Theory and Methods, 33, 1733–1753.
Roitershtein, A. (2007): “One-Dimensional Linear Recursions with MarkovDependent Coefficients,” Annals of Applied Probability, 17, 572–608.
Rutherford, R. S. G. (1955): “Income Distributions: A New Model,” Econometrica, 23, 277–294.
Schumacher, J. M. (1986): “Residue Formulas for Meromorphic Matrices,”
in Computational and Combinatorial Methods in Systems Theory, ed. by C. I.
Byrnes and A. Lindquist, North-Holland, 97–111.
——— (1991): “System-Theoretic Trends in Econometrics,” in Mathematical System Theory: The Influence of R. E. Kalman, ed. by A. C. Antoulas, Berlin:
Springer, 559–577.
Shakai Kougaku Kenkyujo (1974): Nihon Rettou Ni Okeru Jinkou Bunpu
No Chouki Jikeiretsu Bunseki (Long Time Series Analysis of Population Distribution in Japan), Shakai Kougaku Kenkyujo (Social Engineering Research
Institute).

44

Simon, H. A. and C. P. Bonini (1958): “The Size Distribution of Business
Firms,” American Economic Review, 48, 607–617.
Steinberg, S. (1968): “Meromorphic Families of Compact Operators,” Archive
for Rational Mechanics and Analysis, 31, 372–379.
Sutton, J. (1997): “Gibrat’s Legacy,” Journal of Economic Literature, 35, 40–
59.
Toda, A. A. (2011): “Income Dynamics with a Stationary Double Pareto Distribution,” Physical Review E, 83, 046122.
——— (2012): “The Double Power Law in Income Distribution: Explanations
and Evidence,” Journal of Economic Behavior and Organization, 84, 364–381.
——— (2014): “Incomplete Market Dynamics and Cross-Sectional Distributions,”
Journal of Economic Theory, 154, 310–348.
——— (2017a): “Huggett Economies with Multiple Stationary Equilibria,” Journal of Economic Dynamics and Control, 84, 77–90.
——— (2017b): “A Note on the Size Distribution of Consumption: More Double
Pareto than Lognormal,” Macroeconomic Dynamics, 21, 1508–1518.
Toda, A. A. and K. Walsh (2015): “The Double Power Law in Consumption
and Implications for Testing Euler Equations,” Journal of Political Economy,
123, 1177–1200.
Toda, A. A. and K. J. Walsh (2017): “Fat Tails and Spurious Estimation of
Consumption-Based Asset Pricing Models,” Journal of Applied Econometrics,
32, 1156–1177.
Vahrenkamp, R. (1976): “Derivatives of the Dominant Root,” Applied Mathematics and Computation, 2, 29–39.
Vermeulen, P. (2017): “How Fat is the Top Tail of the Wealth Distribution?”
Review of Income and Wealth.
Wang, N. (2003): “Caballero Meets Bewley: The Permanent-Income Hypothesis
in General Equilibrium,” American Economic Review, 93, 927–936.
Widder, D. V. (1941): The Laplace Transform, Princeton, NJ: Princeton University Press.
Wold, H. O. A. and P. Whittle (1957): “A Model Explaining the Pareto
Distribution of Wealth,” Econometrica, 25, 591–595.
Yaari, M. E. (1965): “Uncertain Lifetime, Life Insurance, and the Theory of the
Consumer,” Review of Economic Studies, 32, 137–150.

45

Appendices
A

Proof of Tauberian theorem

The proof of Theorem 2.1 uses a technique developed in Graham and Vaaler
(1981). Define

e−x for x ≥ 0,
E(x) =
0
for x < 0.

Fix λ > 0, let ω = 2π/λ, and define the continuous functions Kλ (x), kλ (x) by
2 "X

#
∞
∞
X
sin λx/2
1
1
e−nω
Kλ (x) =
,
−
e−nω
−
2
λ/2
(x
−
nω)
x
−
nω
x
n=1
n=0
2

sin λx/2
,
kλ (x) = Kλ (x) −
λx/2


for all real x that are not integer multiples of ω. Note that this validly defines
Kλ (x) and kλ (x) for all real x since sin λx/2 ∼ (−1)n λ(x − nω)/2 around x = nω
for each integer n. We have the following lemma.
Lemma A.1. Kλ and kλ satisfy the following properties.
(i) kλ (x) ≤ E(x) ≤ Kλ (x) for all real x.
b λ and b
(ii) Kλ and kλ are absolutely integrable, with Fourier transforms K
kλ
supported on [−λ, λ] and satisfying
b λ (0) =
K

ω
,
1 − e−ω

Proof. See Korevaar (2004, Proposition 5.2).

b
kλ (0) =

eω

ω
.
−1

The following proof of Theorem 2.1 is an adaptation of arguments appearing
in Graham and Vaaler (1981), Korevaar (2004, ch. 5), and Nakagawa (2007).
Proof of Theorem 2.1. Let F be the cdf of X, and fix σ ∈ (0, α). Since E(σ(y−
x))eσy = eσx 1(y ≥ x), we have
σx

e P(X ≥ x) =

Z

∞

E(σ(y − x))eσy dF (y).

−∞

46

Noting that |E(σ(y − x))eσy | ≤ 1 ∨ eαx and |E(α(y − x))eσy | ≤ eαx , two applications of the dominated convergence theorem with dominating function 1 ∨ eαx
(constant as a function of y) reveal that
lim
σ↑α

Z

Z

∞
σy

∞

E(α(y − x))eαy dF (y)
−∞
Z ∞
E(α(y − x))eσy dF (y).
= lim

E(σ(y − x))e

dF (y) =

−∞

σ↑α

Therefore,
αx

e P(X ≥ x) = lim
σ↑α

Z

−∞

∞

E(α(y − x))eσy dF (y).

(A.1)

−∞

For any λ > 0, Lemma A.1(i) implies that
Z

∞
σy

E(α(y − x))e

dF (y) ≤

−∞

Z

∞

Kλ (α(y − x))eσy dF (y).

(A.2)

−∞

Lemma A.1(ii) implies that Kλ (α(y−x))eσy , viewed as a function of x, is absolutely
integrable with Fourier transform
Z

∞

σy −itx

Kλ (α(y − x))e e

−∞

Z
1 ∞
dx =
Kλ (w)eσy e−it(y−w/α) dw
α −∞
1b
(σ−it)y
.
= K
λ (t/α)e
α

Therefore, applying Fubini’s theorem, we find that the Fourier transform of the
upper bound in (A.2) is given by
Z

∞
−∞

1b
1b
Kλ (t/α)e(σ−it)y dF (y) = K
λ (t/α)M(σ − it).
α
α

This Fourier transform has support [−αλ, αλ] by Lemma A.1(ii) so, applying the
Fourier inversion theorem, we find that the upper bound in (A.2) satisfies
Z

∞
σy

Kλ (α(y − x))e
−∞

1
dF (y) =
2π

Z

αλ

−αλ

1b
Kλ (t/α)M(σ − it)eitx dt.
α

(A.3)

In deriving (A.3) we have only assumed that the Laplace transform of F has right
abscissa of convergence α ∈ (0, ∞). It must therefore be valid for the cdf

G(x) =


1 − e−αx
0

47

for x ≥ 0
for x < 0,

R∞

whose Laplace transform is
Z

∞
σy

Kλ (α(y − x))e
−∞

ezx dG(x) =

−∞

1
dG(y) =
2π

Z

αλ

−αλ

α
.
α−z

In this case (A.3) specializes to

α
1 b
Kλ (t/α)
eitx dt.
α
α − σ + it

(A.4)

Combining (A.2), (A.3) and (A.4), we obtain
Z

∞
σy

E(α(y − x))e

−∞



Z αλ
1
A
1 b
dF (y) ≤
Kλ (t/α) M(σ − it) +
eitx dt
2π −αλ α
σ − it − α
Z
A ∞
Kλ (α(y − x))eσy dG(y)
+
α −∞
=: I1 (σ, x) + I2 (σ, x).

Let H(z) denote the continuous extension of M(z) + A(z − α)−1 to SB+ , which
exists due to the definition of B. Fix b ∈ (0, B). Since H(z) is uniformly continuous on the compact set
{z ∈ C : 0 ≤ Re z ≤ α, −b ≤ Im z ≤ b} ,
it must be the case that
lim sup |H(σ − it) − H(α − it)| = 0.
σ↑α t∈[−b,b]

Consequently, setting λ = b/α we obtain
1
lim I1 (σ, x) =
σ↑α
2π

Z

b

−b

1b
Kb/α (t/α)H(α − it)eitx dt,
α

and the Riemann-Lebesgue lemma then yields limx→∞ limσ↑α I1 (σ, x) = 0. Next,
applying the dominated convergence theorem we obtain
A
lim I2 (σ, x) =
σ↑α
α

Z

∞

Kλ (α(y − x))eαy dG(y)
Z −∞
Z
∞
A ∞
=A
Kλ (α(y − x)) dy =
Kλ (w) dw.
α −αx
0

Letting x → ∞ and noting that
A.1(ii), we obtain

R∞

−∞

b λ (0) = ω/(1 − e−ω ) by Lemma
Kλ (w) dw = K

lim lim I2 (σ, x) =

x→∞ σ↑α

A ω
2πA/b
=
.
−ω
α1−e
1 − e−2πα/b
48

Recalling (A.1), this establishes that
lim sup eαx P(X > x) ≤
x→∞

2πA/b
.
1 − e−2πα/b

Since b can be chosen arbitrarily close to B, the claimed upper bound in (2.4)
follows. The proof of the lower bound is similar. Sharpness of the bounds was
demonstrated in Remark 2.2.

B

Properties of nonnegative matrices

In this appendix we collect properties of nonnegative matrices that are used
throughout the paper. For a square (complex) matrix A, let ρ(A) denote its spectral radius, i.e., the largest modulus of all eigenvalues, |A| be the matrix obtained
by taking the modulus of each element of A, and kAk be any matrix norm.
For matrices A, B of the same size, let A ⊙ B denote the Hadamard (entrywise) product, so for A = (amn ) and B = (bmn ), we have A ⊙ B = (amn bmn ). For
a nonnegative matrix A = (amn ), let A(α) = (aαmn ) denote the matrix of entry-wise
power. The following proposition shows that the spectral radius has a convexity
property with respect to the Hadamard product.
Proposition B.1 (Theorem 1, Elsner et al., 1988). Let A, B be nonnegative square
matrices of the same size and θ ∈ (0, 1). Then
ρ(A(1−θ) ⊙ B (θ) ) ≤ ρ(A)(1−θ) ρ(B)θ .
We call an N × N square matrix A = (amn ) irreducible if for any m, n,
there exist numbers m = k1 , k2 , . . . , kp = n (1 ≤ p ≤ N) such that the entries ak1 k2 , ak2 k3 , . . . , akp−1 kp are all nonzero. There are many different ways to
characterize irreducibility: see Theorem 6.2.24 of Horn and Johnson (1985).
Proposition B.2. For N × N complex matrices A, B, the followings are true.
(i) ρ(A) = limk→∞ Ak

1
k

.

(ii) If |B| ≤ A, then ρ(B) ≤ ρ(|B|) ≤ ρ(A).
(iii) If A is nonnegative and irreducible, |B| ≤ A, ρ(A) = ρ(B), and λ =
eiφ ρ(B) is an eigenvalue of B, then there exist θ1 , . . . , θN ∈ R such that
B = eiφ DAD −1 , where D = diag(eiθ1 , . . . , eiθN ).
49

Proof. Property (i) (which is known as Gelfand’s spectral radius formula) is Corollary 5.6.14 of Horn and Johnson (1985). Property (ii) is Theorem 8.1.18 of Horn and Johnson
(1985). Property (iii) is Theorem 8.4.5 of Horn and Johnson (1985).
Perron-Frobenius theorem. Let A be a square, nonnegative, and irreducible
matrix. Then
(i) ρ(A) > 0,
(ii) ρ(A) is an eigenvalue of A (which is called the Perron root),
(iii) there exist positive vectors x, y (called the right and left Perron vectors) such
that Ax = ρ(A)x and y ⊤ A = ρ(A)y ⊤ , and
(iv) ρ(A) is an algebraically (and hence geometrically) simple eigenvalue of A.
Proof. See Theorem 8.4.4 of Horn and Johnson (1985).
Proposition B.3. Let A be a square, nonnegative, and irreducible matrix, and B
be a complex matrix of the same dimension. Then ρ(A + zB) depends holomorphically on z in a neighborhood of z = 0.
Proof. See Vahrenkamp (1976).
Proposition B.4 (Theorem 5.2, Karlin, 1982). Let Π be a nonnegative, irreducible
matrix with rows summing to one, D be a positive diagonal matrix, and f (t) =
ρ((tI + (1 − t)Π)D) for t ∈ [0, 1]. Then f ′ (t) ≥ 0.

C

Simple poles of matrix pencil inverses

In this appendix we collect some properties of matrix pencil inverses that are used
to prove the main results of our paper. A matrix pencil is a square (complex)
matrix-valued function of a complex variable.
Proposition C.1. Let A(z) be an n × n matrix pencil depending holomorphically
on z ∈ Ω, where Ω is some open and connected subset of the complex plane.
Suppose that A(z) is invertible for some z ∈ Ω. Then A(z) has a meromorphic
inverse A(z)−1 on Ω, with poles at the points of noninvertibility of A(z). If A(z0 )
has rank r < n for some z0 ∈ Ω, so that A(z)−1 has a pole at z = z0 , then the
following four conditions are equivalent.
(i) The pole at z0 is simple.
50

(ii) The geometric and algebraic multiplicities of the unit eigenvalue of I − A(z0 )
are equal.
(iii) The r × r matrix y ⊤ A′ (z0 )x is invertible, where x and y can be any n × r
matrices of full column rank such that A(z0 )x = 0 and y ⊤ A(z0 ) = 0.
(iv) The complex vector space Cn is the direct sum of the column space of A(z0 )
and the image of the null space of A(z0 ) under A′ (z0 ).
Under any of these equivalent conditions, the residue of A(z)−1 at the simple pole
at z = z0 is equal to x(y ⊤ A′ (z0 )x)−1 y ⊤ .
Proof. Meromorphicity of A(z)−1 when A(z) is somewhere invertible was proved
by Steinberg (1968). The equivalence of (i), (ii), and (iv) was proved by Howland
(1971). Note that both of the authors just cited worked in a more general Banach
space setting. The equivalence of (i) and (iii) and the residue formula were proved
by Schumacher (1986); see also Schumacher (1991, pp. 562–563).
Remark C.1. Proposition C.1 is closely related to the Granger-Johansen representation theorem (Engle and Granger, 1987; Johansen, 1991). The connection
was first commented upon by Schumacher (1991), who observed that the condition
used by Johansen (1991) to guarantee that the solution to a vector autoregressive equation is integrated of order one, or I(1), corresponded to condition (iii).
Johansen (1995, Corollary 4.3) gave a reformulation of the I(1) condition corresponding to condition (ii), while Beare et al. (2017) gave a reformulation of the
I(1) condition corresponding to condition (iv), and used it to develop an extension of the Granger-Johansen representation theorem to a general Hilbert space
setting. Condition (iv) is not used to prove any of the results in this paper but
we have included it for the sake of completeness.

D

Proofs of results in Section 4

Proof of Proposition 4.1. We divide the proof into three steps.
Step 1. If (4.5) has a solution, then the policy (4.4) satisfies the Bellman equation.
We prove by guess-and-verify. Substituting (4.4a) into the Bellman equation
(4.3), we obtain


1 −γ(aw+bs )
1 −γc
β  −γ(aR(w−c+ȳs )+bs′ ) 
− e
= max − e
−
E e
s .
c
γa
γ
γa
51

(D.1)

The first-order condition with respect to c is


e−γc − βR E e−γ(aR(w−c+ȳs )+bs′ ) s = 0.

(D.2)

Substituting (D.2) into (D.1), we obtain

1
1
− e−γ(aw+bs ) = −
γa
γa



1
a+
e−γc .
R

(D.3)

Comparing the coefficients, (D.3) trivially holds if a = 1 − 1/R and c = aw + bs .
In this case (D.1) holds, and so does the Bellman equation (4.3). To determine
b = (bs )Ss=1 , by the first-order condition we have


(D.2) ⇐⇒ e−γ(aw+bs ) = βR E e−γ(aw+(R−1)(−bs +ȳs )+bs′ ) s


⇐⇒ 1 = βR E e−γ(−Rbs +(R−1)ȳs +bs′ ) s




1
1
ȳs −
log E βRe−γbs′ s
⇐⇒ bs = 1 −
R
γR
!


S
X
1
1
⇐⇒ bs = 1 −
ȳs −
log βR
pss′ e−γbs′ .
R
γR
′
s =1

Expressing in matrix form, we obtain (4.5).
Step 2. The system of equations (4.5) (for s = 1, . . . , S) admits a unique solution
b = (bs )Ss=1 .
For b ∈ RS , define T : RS → RS by



1
1
Tb = 1 −
y−
log βRP e−γb .
R
γR
Then (4.5) is equivalent to T b = b, i.e., b is a fixed point of T . To show that T has
a unique fixed point, we show that T is a contraction using the sufficient condition
in Blackwell (1965, Theorem 5). Clearly b1 ≤ b2 =⇒ T b1 ≤ T b2 , so monotonicity
holds. Furthermore, for any t ≥ 0, we have


1
T (b + t1) = 1 −
y−
R


1
y−
= 1−
R
t
= T b + 1.
R


1
log βRP e−γ(b+t1)
γR
 t
1
log βRP e−γb + 1
γR
R

52

Since R > 1, T satisfies discounting, and hence by Blackwell’s condition T is a
contraction. By the contraction mapping theorem, T has a unique fixed point.
Step 3. The policy (4.4) characterizes the solution to the utility maximization
problem.
Since the value function (4.4a) satisfies the Bellman equation, in order to
demonstrate the optimality of the consumption rule (4.4b), it remains to verify
the transversality condition
lim β t E0 [Vst (wt )] = 0,

t→∞

(D.4)

where wt is the wealth at time t determined by the budget constraint (4.2) and
the consumption rule (4.4b).17 Using the formula for the value function (4.4a) and
the budget constraint (4.2), the continuation value in the next period is
Vs′ (w ′ ) = −

1
1 −γ(aw′ +bs′ )
e
= − e−γ(aR(w−c+ȳs )+bs′ ) .
γa
γa

Taking the expectation conditional on s, by (D.2) we obtain
E [Vs′ (w ′) | s] = −

1
1 −γ(aw+bs )
1
e−γc = −
e
=
Vs (w).
γaβR
γaβR
βR

Iterating this equation from t = 0 to t, we obtain
β t E0 [Vst (wt )] =

1
Vs (w0 ).
Rt 0

Since R > 1, letting t → ∞ we obtain the transversality condition (D.4).
Proof of Theorem 4.2. For notational simplicity, assume p = 0. The case p > 0
R
instead of β, R. We
is completely analogous by using β̃ = β(1 − p) and R̃ = 1−p
divide the proof into three steps.
Step 1. If an equilibrium exists, then βR ≤ 1.
By (4.5), we have
!


S
X
1
1
pss′ e−γbs′ .
bs = 1 −
ȳs −
log βR
R
γR
s′ =1
17

See Kamihigashi (2014) and references therein for the sufficiency of the transversality condition.

53

Since (pss′ )Ss′ =1 are conditional probabilities, they are nonnegative and sum to 1.
Since − log(·) is strictly convex, it follows that


S
1 X
1
ȳs −
pss′ log(βRe−γbs′ )
bs ≤ 1 −
R
γR ′
s =1
!


S
X
1
1
= 1−
ȳs −
log(βR) +
pss′ (−γbs′ )
R
γR
s′ =1


1
1
⇐⇒ b ≤ 1 −
y−
(log(βR)1 − γP b) ,
R
γR
where the last inequality is element-by-element. Left-multiplying the stationary
distribution π ≫ 0 as an inner product, it follows that


1
1
π b≤ 1−
π⊤y −
(log(βR) − γπ ⊤ P b)
R
γR


1
1
π⊤y −
(log(βR) − γπ ⊤ b)
= 1−
R
γR
1
⇐⇒ π ⊤ b ≤ π ⊤ y −
log(βR).
γ(R − 1)
⊤

(∵ π ⊤ = π ⊤ P )
(D.5)

By the equilibrium condition (4.7), we have π ⊤ b = π ⊤ y, so
0≤−

1
log(βR) ⇐⇒ βR ≤ 1
γ(R − 1)

since R > 1 is necessary for equilibrium existence.
Step 2. Let b(R) be the value of b = (bs )Ss=1 implied by (4.5). Then mins bs (R) →
∞ as R ↓ 1.
Let b = mins bs (R) and y = mins ȳs . Using R > 1 and the monotonicity of T ,
¯
¯
it follows from (4.5) that


S
1
1 X
b≥ 1−
y−
pss′ log(βRe−γb¯)
¯
R ¯ γR ′
s =1


1
1
= 1−
y−
(log(βR) − γb)
¯
R ¯ γR
1
log(βR) → ∞
⇐⇒ b ≥ y −
¯ ¯ γ(R − 1)
as R ↓ 1 because β < 1.
Step 3. An equilibrium exists. The gross risk-free rate satisfies 1 − p < R ≤ 1/β.
54

Let b(R) be the value of b = (bs )Ss=1 implied by (4.5). By (4.7), the equilibrium
condition is g(R) := π ⊤ b(R) − π ⊤ y. By Proposition 4.1, b(R) is well-defined for
R > 1, and it is smooth by the implicit function theorem. Therefore g(R) is welldefined and continuous for R > 1. By the previous step, we have limR↓1 g(R) = ∞.
Letting R = 1/β in (D.5), we obtain
g(1/β) ≤ −

1
log(1) = 0.
γ(1/β − 1)

By the intermediate value theorem, there exists R ∈ (1, 1/β] such that g(R) = 0.
If p > 0, by the same argument as above we obtain 1 < R̃ ≤ 1/β̃ ⇐⇒ 1 − p <
R ≤ 1/β.

E

Japanese population data and estimation

E.1

Historical background and data

During the Edo era (1603–1868), Japan was divided into provinces called han,
which were controlled by feudal lords called daimyō. The movement of people
across regions was severely restricted during this period. In 1868, the government power was transferred from the Tokugawa shogunate to the Emperor (Meiji
Restoration) and Japan transitioned to a modern nation with a market-based
economy. In 1871 the Meiji government abolished the han system (haihan-chiken).
Hans were reorganized into prefectures. Initially (August 1871) there were 305
prefectures (302 kens plus Tokyo, Osaka, Kyoto), but they were consolidated into
75 by November 1871. The prefectures were further consolidated into 72 in 1872,
63 in 1873, 62 in 1875, and 38 in 1876. In 1889, some of them were divided and
the total number became 46 plus Hokkaido, which is the current number.
There are 47 modern prefectures in Japan. Since 1920, the census has been
conducted every five years. The early census was irregular, conducted in 1873,
1884, 1893, 1908, and 1913. Special census was conducted in 1944, 1945, 1946,
and 1947 due to WWII and its aftermath. Statistics Japan18 computes prefecture
population data for all years (as of October 1) since 1920 by aggregating the
residence registration data in each city. Since completing residence registration is
a prerequisite for receiving government service, the quality of this data is high.
The population data come in spreadsheets for 1920-2000, 2000-2010, and each year
thereafter. Combining these spreadsheets, we can construct a balanced panel of
18

http://www.stat.go.jp/data/jinsui/2.htm

55

annual prefecture population data from 1920 to 2015 except for the prefecture of
Okinawa, which was occupied by U.S. from 1945 to 1972.
Before 1920, the boundaries of prefectures changed several times. Using the
census data in 1873, 1884, 1893, 1908, and 1913, the Japanese urban planning
book Shakai Kougaku Kenkyujo (1974) has imputed the prefecture population
data with the 1920 boundaries. Due to data quality concerns, we use this data
only for the cross-sectional estimation and not the panel one (which is our focus).
Figure E.1 plots the population of selected prefectures over time.

Population

10 8

10 7

Hokkaido
Niigata
Tottori
Tokyo

10 6

10 5

1880 1900 1920 1940 1960 1980 2000

Year

Figure E.1: Population of selected prefectures.

The prefectures are Hokkaido (smallest in 1873), Niigata (largest in 1873),
Tottori (smallest in 2015), and Tokyo (largest in 2015). Hokkaido (the northern
island) experienced high growth before WWII, which is largely due to migration.
(Before the Meiji Restoration in 1868, Hokkaido was an undeveloped territory
inhabited by the Ainu people, who were hunter-gatherers. After 1868 the Japanese
government encouraged migration and development.) Tokyo’s population has two
dips, one in 1923 (due to the earthquake on September 1 that killed more than
100,000) and another in 1945 (due to refugees and death during WWII). Niigata
is a north-central prefecture whose main industry is agriculture (rice). It was the
largest prefecture in 1873 but has experienced little growth (except in 1945 when
it received war refugees from other urban areas). The graph for Tottori (a rural
western prefecture) is almost parallel to Niigata.

56

E.2

Cross-sectional estimation

A Pareto-lognormal (PlN) random variable is the product of a lognormal variable
and an independent Pareto variable with minimum size 1 (normalization) and
Pareto exponent α > 0 (Reed and Jorgensen, 2004). The distribution of the
logarithm of a PlN variable is the convolution of the normal and exponential
distributions, which is known as the normal-exponential distribution. It has three
parameters, the mean µ and standard deviation σ > 0 of the normal component
and the exponent α > 0 of the exponential component. The density of the normalexponential distribution is given by
−α(x−µ)+α2 σ2 /2

f (x; µ, σ, α) = αe




x−µ
Φ
− ασ ,
σ

where Φ(·) is the cdf of N(0, 1). Maximum likelihood estimation is straightforward
using this density.

E.3

Panel estimation

We estimate the hidden Markov model using the Hamilton (1989) filter. The
following algorithm computes the log-likelihood. (The algorithm considers the
case of one time series, but the panel case is similar.)
(i) Given the transition probability matrix Π, compute the stationary distribution π as the right Perron vector of Π that sums to 1. Set ξ0|0 = π.
(ii) For t = 1, . . . , T , do the following:
(a) Let ξt|t−1 = Π⊤ ξt−1|t−1 be the one period ahead prior forecast.
(b) Define the vector of conditional likelihoods ηt = (ηnt )N
n=1 by
2

ηnt

(log Gt −µn )
1
−
2
2σn
e
,
=p
2
2πσn

where log Gt is the realized log growth rate.

(c) Update the one period ahead forecast by ξt|t = (ηt ⊙ ξt|t−1 )/(ηt⊤ ξt|t−1 ),
where ⊙ denotes the Hadamard (element-wise) product.
(iii) Compute the log likelihood as log L =

PT

t=1

log(ηt⊤ ξt|t−1 ).

We then estimate the parameters Π, µ = (µ1 , . . . , µN )⊤ , and σ = (σ1 , . . . , σN )⊤
by maximizing the log-likelihood.
57

