$\displaystyle=\psi(y)\quad\text{for all }y\in\mathsf{Y}$ (3.21)

The constraints in (3.20)–(3.21) require that

  1. (i)

for any $x\in\mathsf{X}$, the total amount of probability mass flowing out of
$x$ is $\varphi(x)$ and

  2. (ii)

for any $y\in\mathsf{Y}$, the total amount of probability mass flowing into
$y$ is $\psi(y)$.

In the present setting, a coupling is also called a transport plan. The
constraints (3.20)–(3.21), which define a coupling $\pi$ of $\varphi$ and
$\psi$, generalize the Monge constraint in (3.18). To see this, let $T$ be a
map satisfying (3.18) and set

$\pi(x,y)=\varphi(x)\mathbbm{1}\\{T(x)=y\\}$

on $\mathsf{X}\times\mathsf{Y}$, so that $\pi$ is the joint distribution that
puts all mass on the image of $T$. Prove that (3.20)–(3.21) both hold. Let
$\Pi(\varphi,\psi)$ be the set of all couplings of $\psi$ and $\varphi$.
Taking $\varphi,\psi$ and the cost function $c$ as given, the general
Monge–Kantorovich problem, also called the optimal transport problem, is to
solve

$P:=\min_{\pi}\left\langle{c,\pi}\right\rangle\quad\text{subject
to}\quad\pi\in\Pi(\varphi,\psi).$ (3.22)

where

$\left\langle{c,\pi}\right\rangle:=\sum_{x}\sum_{y}c(x,y)\pi(x,y)$

is the Frobenius inner product of $c$ and $\pi$, treated as
$|\mathsf{X}|\times|\mathsf{Y}|$ matrices. The sum in
$\left\langle{c,\pi}\right\rangle$ measures the total cost of transporting
$\varphi$ into $\psi$ under the plan $\pi$. There is linearity embedded in
this cost formulation, since doubling the amount sent from $x$ to $y$ scales
the associated cost at rate $c(x,y)$. We call any $\pi$ solving (3.22) an
optimal plan. Since we are maximizing over a finite set, at least one such
plan exists.

######  Remark 3.3.1.

The problem (3.22) is sometimes expressed in terms of random variables, as
follows. In this setting, a coupling $\pi$ in $\Pi(\varphi,\psi)$ is
identified with a pair of random elements $(X,Y)$ such that
$X\stackrel{{\scriptstyle d}}{{=}}\varphi$ and $Y\stackrel{{\scriptstyle
d}}{{=}}\psi$. We can then write

$P=\min_{(X,Y)}\mathbbm{E}\,c(X,Y)\quad\text{subject
to}\quad(X,Y)\in\mathscr{D}(\mathsf{X}\times\mathsf{Y})\text{ with
}X\stackrel{{\scriptstyle d}}{{=}}\varphi\text{ and }Y\stackrel{{\scriptstyle
d}}{{=}}\psi.$

One of the most important features of the Kantorovich relaxation is that, for
given $\varphi$ and $\psi$, the constraint set is convex. To verify this, we
let $n=|\mathsf{X}|$ and $m=|\mathsf{Y}|$, associate each $x\in\mathsf{X}$
with some $i\in[n]$, associate each $y\in\mathsf{Y}$ with some $j\in[m]$, and
treat $c$ and $\pi$ as $n\times m$ matrices, with typical elements $c_{ij}$
and $\pi_{ij}$. The constraints are

$\pi\mathbbm{1}_{m}=\varphi\quad\text{and}\quad\pi^{\top}\mathbbm{1}_{n}=\psi,$
(3.23)

where $\mathbbm{1}_{k}$ is an $k\times 1$ vector of ones. With this notation,
prove that the set $\Pi(\varphi,\psi)$ of $\pi\in\mathbbm{M}^{n\times m}$
satisfying the constraints is convex, in the sense that

$\pi,\hat{\pi}\in\Pi(\varphi,\psi)\text{ and
}\alpha\in[0,1]\;\implies\;\alpha\pi+(1-\alpha)\hat{\pi}\in\Pi(\varphi,\psi).$

######  3.3.1.4 Optimal Transport as a Linear Program

With some relatively simple manipulations, the general optimal transport
problem can be mapped into a standard equality form linear program. This
provides two significant benefits. First, we can apply duality theory, which
yields important insights. Second, on the computational side, we can use
linear program solvers to calculate optimal plans. To map the optimal
transport problem into a linear program, we need to convert matrices into
vectors. Will use the $\operatorname{vec}$ operator, which takes an arbitrary
$A\in\mathbbm{M}^{n\times m}$ and maps it to a vector in $\mathbbm{R}^{nm}$ by
stacking its columns vertically. For example,

$\operatorname{vec}\begin{pmatrix}a_{11}&a_{12}\\\
a_{21}&a_{22}\end{pmatrix}=\begin{pmatrix}a_{11}\\\ a_{21}\\\ a_{12}\\\
a_{22}\end{pmatrix}.$

In this section we adopt the notational conventions in Exercise 3.3.1.3. The
objective function $\left\langle{c,\pi}\right\rangle$ for the optimal
transport problem can now be expressed
$\operatorname{vec}(c)^{\top}\operatorname{vec}(\pi)$. To rewrite the
constraints in terms of $\operatorname{vec}(\pi)$, we use the Kronecker
product, which is denoted by $\otimes$ and defined as follows. Suppose $A$ is
an $m\times s$ matrix with entries $(a_{ij})$ and that $B$ is an $n\times t$
matrix. The Kronecker product $A\otimes B$ of $A$ and $B$ is the $mn\times st$
matrix defined, in block matrix form, by

$A\otimes B=\begin{pmatrix}a_{11}B&a_{12}B&\dots&a_{1s}B\\\
a_{21}B&a_{22}B&\dots&a_{2s}B\\\ &&\vdots&\\\ a_{m1}B&a_{m2}B&\dots&a_{ms}B\\\
\end{pmatrix}.$

It can be shown that Kronecker products and the $\operatorname{vec}$ operator
are connected by the following relationship: for conformable matrices $A$, $B$
and $M$, we have

$\operatorname{vec}(AMB)=(B^{\top}\otimes A)\operatorname{vec}(M).$ (3.24)

Using (3.24) and the symbol $I_{k}$ for the $k\times k$ identity matrix, we
can rewrite the first constraint in (3.23) as

$\varphi=I_{n}\pi\mathbbm{1}_{m}=\operatorname{vec}(I_{n}\pi\mathbbm{1}_{m})=(\mathbbm{1}_{m}^{\top}\otimes
I_{n})\operatorname{vec}(\pi).$ (3.25)

Show that the second constraint in (3.23) can be expressed as

$\psi=(I_{m}\otimes\mathbbm{1}_{n}^{\top})\operatorname{vec}(\pi).$ (3.26)

Applying (3.24) on page 3.24, we have

$\psi=\operatorname{vec}(\psi^{\top})=\operatorname{vec}(\mathbbm{1}_{n}^{\top}\pi
I_{m})=(I_{m}\otimes\mathbbm{1}_{n}^{\top})\operatorname{vec}(\pi).$

Now, using block matrix notation and setting

$A:=\begin{pmatrix}\mathbbm{1}_{m}^{\top}\otimes I_{n}\\\
I_{m}\otimes\mathbbm{1}_{n}^{\top}\end{pmatrix}\quad\text{and}\quad
b:=\begin{pmatrix}\varphi\\\ \psi\end{pmatrix},$

the optimal transport problem can be expressed as the standard equality form
linear program

$\min_{x}\operatorname{vec}(c)^{\top}x\quad\text{ over
}x\in\mathbbm{R}^{nm}_{+}\text{ such that }Ax=b.$ (3.27)

Finally, for a given solution $x$, the transport plan is recovered by
inverting $x=vec(\pi)$.

######  3.3.1.5 Implementation

Listing 3.13 is a function that implements the above steps, given flat (one-
dimensiona) arrays phi and psi representing the distributions over the source
and target locations, plus a two-dimensional array c representing transport
costs. (The method argument highs-ipm tells linprog to use a particular
interior point method, details of which can be found in the linprog
documentation. Simplex and other methods give similar results.)

⬇

import numpy as np

from scipy.optimize import linprog

\pardef ot_solver(phi, psi, c, method=’highs-ipm’):

"""

Solve the OT problem associated with distributions phi, psi

and cost matrix c.

\parParameters

———-

phi : 1-D array

Distribution over the source locations.

psi : 1-D array

Distribution over the target locations.

c : 2-D array

Cost matrix.

"""

n, m = len(phi), len(psi)

\par# vectorize c

c_vec = c.reshape((m * n, 1), order=’F’)

\par# Construct A and b

A1 = np.kron(np.ones((1, m)), np.identity(n))

A2 = np.kron(np.identity(m), np.ones((1, n)))

A = np.vstack((A1, A2))

b = np.hstack((phi, psi))

\par# Call solver

res = linprog(c_vec, A_eq=A, b_eq=b, method=method)

\par# Invert the vec operation to get the solution as a matrix

pi = res.x.reshape((n, m), order=’F’)

return pi

Figure 3.13:  Function to solve a transport problem via linear programming

Notice that in Listing 3.13, the reshape order is specifed to F. This tells
NumPy to reshape with Fortran column-major order, which coincides with the
definition of the vec operator described in §3.3.1.4. (The Python vectorize
operation defaults to row-major order, which concatenates rows rather than
stacking columns. In contrast, Julia uses column-major by default.) Let’s call
this function for the very simple problem

$\varphi=\begin{pmatrix}0.5\\\ 0.5\end{pmatrix},\quad\psi=\begin{pmatrix}1\\\
0\end{pmatrix}\quad\text{and}\quad c=\begin{pmatrix}1&1\\\ 1&1\end{pmatrix}$

With these primitives, all mass from $\varphi_{1}$ and $\varphi_{2}$ should be
sent to $\psi_{1}$. To implement this problem we set

⬇

phi = np.array((0.5, 0.5))

psi = np.array((1, 0))

c = np.ones((2, 2))

and then call ot_solver via

⬇

ot_solver(phi, psi, c)

The output is as expected:

⬇

array([[0.5, 0. ],

[0.5, 0. ]])

######  3.3.1.6 Python Optimal Transport

In the case of Python, the steps above have been automated by the Python
Optimal Transport package, due to Flamary et al., (2021). For the simple
problem from §3.3.1.5 we run

⬇

import ot

ot.emd(phi, psi, c) # Use simplex method via the emd solver

The output is again equal to

⬇

array([[0.5, 0. ],

[0.5, 0. ]])

Figure 3.14 shows an example of an optimal transport problem solved using the
Python Optimal Transport package. The interpretation is similar to Figure
3.12, although the number of vertices is larger. In addition, the edges show
the optimal transport configuration, in the sense that $\pi^{*}$, the optimal
transport plan, is treated as the adjacency matrix for the graph. (The figure
shows the unweighted graph, with an arrow drawn from $\varphi_{i}$ to
$\psi_{j}$ whenever $\pi^{*}_{ij}>0$.) The optimal transport plan is obtained
by converting the transport problem into a linear program, as just described,
and applying the simplex method. Although there are 32 nodes of each type, the
problem is solved by the simplex routine in less than one millisecond.

Figure 3.14:  An optimal transport problem solved by linear programming

######  3.3.1.7 Kantorovich Relaxation and Linear Assignment

We showed in §3.3.1.2 that the linear assignment problem studied in §3.2.1 is
a special case of optimal transport, where $\varphi$ and $\psi$ become
discrete uniform distributions on $[n]$. Moreover, as we have just seen,
Kantorovich’s relaxation method allows us to apply linear programming. This
leads to fast solutions. In the discussion in §3.2.1 we used $n=40$. Let’s
start here with $n=4$, to illustrate the method, and then try with $n=40$. The
matrix $c(i,j)$ of costs, which is the only other primitive, will be generated
randomly as an independent array of uniform random variables: Then we apply
the Python Optimal Transport (POT) library, as in §3.3.1.6. Here is our set
up:

⬇

import numpy as np

import ot

\parn = 4

phi = np.ones(n)

psi = np.ones(n)

We have broken the rule that $\varphi$ and $\psi$ should sum to one. This
could be fixed easily by using np.ones(n)/n instead of np.ones(n), but the POT
library does not care (as long as np.sum(phi) equals np.sum(psi)) and,
moreover, the idea of putting unit mass everywhere is natural, since each
element of $\varphi$ represents one worker, and each element of $\psi$
represents one job. Now we build the cost matrix:

⬇

c = np.random.uniform(size=(n, n))

The output is

⬇

array([[0.03841, 0.32896, 0.55989, 0.41682],

[0.91527, 0.24566, 0.26022, 0.64510],

[0.96275, 0.44089, 0.79274, 0.93065],

[0.40454, 0.87307, 0.43555, 0.54903]])

(For example, the cost of retraining worker 1 for job 2 is 0.32896.) Finally,
we call the solver:

⬇

ot.emd(phi, psi, c)

The output is

⬇

array([[1., 0., 0., 0.],

[0., 0., 1., 0.],

[0., 1., 0., 0.],

[0., 0., 0., 1.]])

This is a permutation matrix, which provides another way to express a
permutation of $[n]$. The first row tells us that worker 1 is assigned to job
1, the second tells us that worker 2 is assigned to job 3, and so on. If we
now set $n=40$ and rerun the code, the line ot.emd(phi, psi, c), which calls
the simplex-based solver, runs in less than 1 millisecond on a mid-range
laptop. This is a remarkable improvement on the $2.5\times 10^{30}$ year
estimate for the brute force solver we obtained in §3.2.1.

######  3.3.1.8 Tight Relaxation

Notice that the solution we obtained for the linear assignment problem using
the simplex method does not split mass, as permitted by the Kantorovich
relaxation. For example, we do not send half of a worker to one job and the
other half to another. This is convenient but why does it hold? While we omit
the details, the basic idea is that a general Kantorovich transport plan is a
bistochastic matrix, and all such matrices can be formed as convex
combinations of permutation matrices. (This is called Birkhoff’s Theorem.) In
other words, the permutation matrices are extreme points of the set of
bistochastic matrices. Moreover, Theorem 3.2.3 tells us that any optimizer of
a linear program will be an extreme point—in this case, a permutation matrix.

#####  3.3.2 Kantorovich Duality

One of the greatest achievements of Kantorovich was to show that the optimal
transport problem can be connected to a dual problem, and how that dual
problem can be used to characterize solutions. This work anticipated much of
the later development of duality theory for arbitrary linear programs.
Throughout this section, in stating the main results, we use the notation

$\left\langle{f,\varphi}\right\rangle=\sum_{x}f(x)\varphi(x)\quad\text{for
}f\in\mathbbm{R}^{\mathsf{X}}\text{ and }\varphi\in\mathscr{D}(\mathsf{X}).$

This is just the usual inner product, when we think of $f$ and $\varphi$ as
vectors in $\mathbbm{R}^{|\mathsf{X}|}$. Also, given a cost function $c$ on
$\mathsf{X}\times\mathsf{Y}$, let $\mathscr{F}_{c}$ be all pairs $(w,p)$ in
$\mathbbm{R}^{\mathsf{X}}\times\mathbbm{R}^{\mathsf{Y}}$ such that

$p(y)\leqslant c(x,y)+w(x)\text{ on }\mathsf{X}\times\mathsf{Y}.$ (3.28)

One part of Kantorovich’s duality results runs as follows.

######  Theorem 3.3.1.

For all $\varphi\in\mathscr{D}(\mathsf{X})$ and
$\psi\in\mathscr{D}(\mathsf{Y})$, we have $P=D$, where

$D:=\max_{(w,p)}\left\\{\left\langle{p,\psi}\right\rangle-\left\langle{w,\varphi}\right\rangle\right\\}\quad\text{subject
to}\quad(w,p)\in\mathscr{F}_{c}.$ (3.29)

Theorem 3.3.1 can now be understood as a special case of the more general
result that strong duality holds for linear programs, which we stated in
Theorem 3.2.1. In that spirit, let us verify Theorem 3.3.1 using Theorem
3.2.1, by working with the linear programming formulation of the optimal
transport problem provided in §3.3.1.4. To do this, we take that formulation,
which is stated in (3.27), and apply the dual formula in (3.13), which yields

$D=\max_{\theta\in\mathbbm{R}^{n+m}}\begin{pmatrix}\varphi\\\
\psi\end{pmatrix}^{\top}\theta\;\text{ subject to
}\;\begin{pmatrix}\mathbbm{1}_{m}^{\top}\otimes I_{n}\\\
I_{m}\otimes\mathbbm{1}_{n}^{\top}\end{pmatrix}^{\top}\theta\leqslant\operatorname{vec}(c).$

If we write the argument $\theta\in\mathbbm{R}^{n+m}$ as $(-w,p)$, so that we
now maximize over the two components $-w\in\mathbbm{R}^{n}$ and
$p\in\mathbbm{R}^{m}$, as well as transposing the constraint, we get

$\max_{w,\,p}\left\\{p^{\top}\psi-w^{\top}\varphi\right\\}\;\text{ subject to
}\;p^{\top}(I_{m}\otimes\mathbbm{1}_{n}^{\top})-w^{\top}(\mathbbm{1}_{m}^{\top}\otimes
I_{n})\leqslant\operatorname{vec}(c)^{\top}$

where $w\in\mathbbm{R}^{n}$ and $p\in\mathbbm{R}^{m}$. By using the definition
of the Kronecker product and carefully writing out the individual terms, it
can be shown that the constraint in this expression is equivalent to requiring
that $p_{j}-w_{i}\leqslant c_{ij}$ for all $(i,j)\in[n]\times[m]$. Recalling
that $\mathsf{X}$ has been mapped to $[n]$ and $\mathsf{Y}$ has been mapped to
$[m]$, this is exactly the same restriction as (3.28). At this point it is
clear that (3.29) is nothing but the dual of the linear program formed from
the optimal transport problem. The claims in Theorem 3.3.1 now follow directly
from the strong duality of linear programs (Theorem 3.2.1). Show that

$\left\langle{c,\pi}\right\rangle\geqslant\left\langle{p,\psi}\right\rangle-\left\langle{w,\varphi}\right\rangle\text{
whenever }\pi\in\Pi(\varphi,\psi)\text{ and }(w,p)\in\mathscr{F}_{c}.$ (3.30)

Use this fact to provide a direct proof that weak duality holds for the
optimal transport problem, in the sense that $D\leqslant P$. (Here $P$ is
defined in the primal problem (3.22) and $D$ is defined in the dual problem
(3.29).) Let $\pi$ be feasible for the primal problem and let $(w,p)$ be
feasible for the dual. By dual feasibility, we have $c(x,y)\geqslant
p(y)-w(x)$ for all $x,y$, so

$\left\langle{\pi,c}\right\rangle\geqslant\sum_{x}\sum_{y}\pi(x,y)[p(y)-w(x)]=\sum_{x}\sum_{y}\pi(x,y)p(y)-\sum_{x}\sum_{y}\pi(x,y)w(x)$

Rearranging and using primal feasibility now gives
$\left\langle{\pi,c}\right\rangle\geqslant\left\langle{p,\psi}\right\rangle-\left\langle{w,\varphi}\right\rangle$.
This proves the first claim. To see that $D\leqslant P$ follows from the last
inequality, just fix $\pi\in\Pi(\varphi,\psi)$ and maximize over all feasible
dual pairs to obtain $\left\langle{\pi,c}\right\rangle\geqslant D$. Now
minimize over $\pi\in\Pi(\varphi,\psi)$.

#####  3.3.3 Optimal Transport and Competitive Equilibria

The other major achievement of Kantorovich in the context of duality theory
for optimal transport was to connect optimality of transport plans with the
existence of functions $w,p$ from the dual problem such that a version of the
complementary slackness conditions hold. Here we present this result, not in
the original direct formulation, but rather through the lens of a competitive
equilibrium problem. In doing so, we illustrate some of the deep connections
between prices, decentralized equilibria and efficient allocations.

######  3.3.3.1 The Advisor’s Problem

We imagine the following scenario. Iron is mined at a finite collection of
sites, which we denote by $\mathsf{X}$. We identify an element
$x\in\mathsf{X}$ with a point $(a,b)\in\mathbbm{R}^{2}$, which can be
understood as the location of the mine in question on a map. At the wish of
the queen, who seeks to defend the empire from greedy rivals, this iron is
converted to swords by blacksmiths. There are a number of talented blacksmiths
in this country, located at sites given by $\mathsf{Y}$. As for $\mathsf{X}$,
each $y\in\mathsf{Y}$ indicates a point in $\mathbbm{R}^{2}$. Henceforth, we
refer to “mine $x$” rather than “the mine at $x$” and so on. Each month, mine
$x$ produces $\varphi(x)$ ounces of iron ore, while blacksmith $y$ consumes
$\psi(y)$ ounces. We take these quantities as fixed. We assume that total
supply equals total demand, so that $\sum_{x}\varphi(x)=\sum_{y}\psi(y)$. For
convenience, we normalize this sum to unity. As a result, $\varphi$ and $\psi$
are elements of $\mathscr{D}(\mathsf{X})$ and $\mathscr{D}(\mathsf{Y})$
respectively. The cost of transporting from $x$ to $y$ is known and given by
$c(x,y)$ per ounce. The king’s chief advisor is tasked with allocating and
transporting iron from the mines to the blacksmiths, such that each blacksmith
$y$ receives their desired quantity $\psi(y)$, at minimum cost. A small amount
of thought will convince you that the advisor’s problem is a version of the
optimal transport problem (3.22). We call this the primal problem in what
follows. Operating in the days before Kantorovich, Dantzig and the electronic
computer, the advisor employs a large team of bean counters, instructing them
to find the allocation with least cost by trying different combinations.
However, after a few days, she realizes the futility of the task. (With
infinite divisibility, which corresponds to our mathematical model, the number
of allocations is infinite. If we replace infinite divisibility with a finite
approximation, the scale can easily be as large as that of the matching
problem discussed in §3.2.1, with only a moderate number of mines and
blacksmiths.)

######  3.3.3.2 The Guild’s Problem

At this point she has another idea. There is a guild of traveling salesmen,
who buy goods in one town and sell them in another. She seeks out the guild
master and asks him to bid for the project along the following lines. The
guild will pay the queen’s treasury $w(x)$ per ounce for iron ore at mine $x$.
It will then sell the iron at price $p(y)$ per ounce to the queen’s
representative at blacksmith $y$. The difference can be pocketed by the guild,
as long as all blacksmiths are provided with their desired quantities. The
guild master is asked to propose price functions $w$ and $p$. The guild master
sees at once that $p$ and $w$ must satisfy $p(y)-w(x)\leqslant c(x,y)$ at each
$x,y$, for otherwise the advisor, who is no ones fool, will see immediately
that money could be saved by organizing the transportation herself. Given this
constraint, the guild master seeks to maximize aggregate profits, which is
$\sum_{y}p(y)\psi(y)-\sum_{x}w(x)\varphi(x)$. At this point it will be clear
to you that the problem of the guild master is exactly that of Kantorovich’s
dual problem, as given in Theorem 3.3.1. Since the advisor has given up on her
team of bean counters, the guild master employs them, and asks them to produce
the optimal pair of prices. The bean counters set to work, trying different
combinations of prices that satisfy the constraints. However, without a
systematic methodology to follow or fast computers to turn to, their progress
is slow. The advisor begins to fear that the coming war will be over before
the guild master replies.

######  3.3.3.3 Decentralized Equilibrium

At this point, it occurs to the advisor that yet another approach exists:
privatize the mines, abolish the guild, and let the traveling salesmen, mine
owners and blacksmiths make individual choices in order to maximize their
profits. Purchase and sales prices, as well as quantities transported from
each mine to each blacksmith, will be determined by the free market. Although
the advisor predates Kantorovich, she reasons that competition will prevent
each salesman from profiteering, while the desire for profits will encourage
high levels of transportation and minimal waste. It turns out that this idea
works amazingly well, in the sense that we now describe. For the record, we
define a competitive equilibrium for this market as pair of price vectors
$(w,p)$ in $\mathbbm{R}^{\mathsf{X}}\times\mathbbm{R}^{\mathsf{Y}}$ and a set
of quantities $\pi\colon\mathsf{X}\times\mathsf{Y}\to\mathbbm{R}_{+}$ such
that the following three conditions hold: For all $(x,y)$ in
$\mathsf{X}\times\mathsf{Y}$,

$\displaystyle\sum_{v\in\mathsf{Y}}\pi(x,v)=\varphi(x)\text{ and
}\sum_{u\in\mathsf{X}}\pi(u,y)=\psi(y)$ (RE) $\displaystyle p(y)\leqslant
c(x,y)+w(x)$ (NA) $\displaystyle p(y)=c(x,y)+w(x)\text{ whenever }\pi(x,y)>0.$
(IC)

Condition (RE) is a resource constraint that builds in the assumption no ore
is wasted or disposed. Condition (NA) imposes no arbitrage. If it is violated
along route $(x,y)$, then another salesman, of which we assume their are many,
will be able to gain business without suffering losses by offering a slightly
higher purchase prices at $x$ or a slightly lower sales prices at $y$.
Finally, condition (IC) is an incentive constraint, which says that, whenever
a route is active (in the sense that a nonzero quantity is transported),
prices are such that the salesmen do not lose money. We do not claim that a
competitive equilibrium will hold immediately and at every instant in time.
However, we reason, as the advisor does, that competitive equilibrium has
natural stability properties, as described in the previous paragraph. As such,
we predict it as a likely outcome of decentralized trade, provided that
private property rights are enforced (e.g., bandits are eliminated from the
routes) and noncompetitive behaviors are prevented (e.g., collusion by mine
owners are met by suitably painful punishments). Taking $c$, $\varphi$ and
$\psi$ as given, we can state the following key theorem, which states that any
competitive equilibrium simultaneously solves both the advisor’s quantity
problem _and_ the guild master’s price problem.

######  Theorem 3.3.2.

If prices $(w,p)$ and $\pi\in\mathbbm{M}^{n\times m}$ form a competitive
equilibrium, then

  1. (i)

$\pi$ is an optimal transport plan, solving the primal problem (3.22), and

  2. (ii)

$(w,p)$ solves the Kantorovich dual problem (3.29).

To prove the theorem we will use the results from the next exercise. Let $A$
and $B$ be nonempty sets. Let $f$ and $g$ be real-valued on $A$ and $B$ such
that $f(a)\geqslant g(b)$ for all $(a,b)\in A\times B$ and, in addition,
$\min_{a\in A}f(a)=\max_{b\in B}g(b)$. Prove the following statement: If there
exist $\bar{a}\in A$ and $\bar{b}\in B$ such that $f(\bar{a})=g(\bar{b})$,
then $\bar{a}$ is a minimizer of $f$ on $A$ and $\bar{b}$ is a maximizer of
$g$ on $B$. From $\min_{a\in A}f(a)=\max_{b\in B}g(b)$ we have $f(a)\geqslant
g(b)$ for all $(a,b)\in A\times B$. Taking $(\bar{a},\bar{b})$ with
$f(\bar{a})=g(\bar{b})$, we have $f(\bar{a})=g(\bar{b})\leqslant f(a)$ for any
given $a\in A$. In particular, $\bar{a}$ minimizes $f$ on $A$. The argument
for $\bar{b}$ is similar. We will use Exercise 3.3.3.3 in the following way.
Let $A=\Pi(\varphi,\psi)$ and $B=\mathscr{F}_{c}$. Let $f$ be the value of the
primal and $g$ be the value of the dual. By (3.30), the ordering
$f(\pi)\geqslant g(w,p)$ holds over all feasible $\pi\in A$ and $(w,p)\in B$
pairs. By strong duality, we also have $\min_{\pi\in A}f(\pi)=\max_{(w,p)\in
B}g(w,p)$. Hence we need only show that, when $(w,p)$ and $\pi$ form a
competitive equilibrium, we have $f(\pi)=g(w,p)$.

###### Proof of Theorem 3.3.2.

Suppose that $(w,p)$ and $\pi$ form a competitive equilibrium. From (RE) we
know that $\pi$ is feasible for the primal problem. From (NA) we know that
$(w,p)$ is feasible for the dual. Since the equality in the (IC) condition
holds when $\pi(x,y)>0$, we can multiply both sides of this equality by
$\pi(x,y)$ and sum over all $x,y$ to obtain

$\sum_{x,y}c(x,y)\pi(x,y)=\sum_{y}p(y)\psi(y)-\sum_{x}w(x)\psi(x).$ (3.31)

The result of Exercise 3.3.3.3 now applies, so $\pi$ attains the minimum in
the primal problem and $(w,p)$ attains the maximum in the dual. ∎

We also have the following converse:

######  Theorem 3.3.3.

If $\pi$ is an optimal transport plan, then there exists a pair
$(w,p)\in\mathbbm{R}^{\mathsf{X}}\times\mathbbm{R}^{\mathsf{Y}}$ such that the
quantities determined by $\pi$ and the prices in $(w,p)$ form a competitive
equilibrium.

###### Proof.

Let $\pi$ be an optimal plan. To be optimal, $\pi$ must be feasible, so
$\pi\in\Pi(\varphi,\psi)$, which implies that (RE) holds. By Kantorovich’s
duality theorem (Theorem 3.3.1, we can obtain $(w,p)\in\mathscr{F}_{c}$ such
that (3.31) holds. Since $(w,p)\in\mathscr{F}_{c}$ (NA) holds. From (NA) we
have $c(x,y)+w(x)-p(y)\geqslant 0$ for all $x,y$. From this and (3.31) we see
that (IC) must be valid. We conclude that $\pi$ and $(w,p)$ form a competitive
equilibrium. ∎

#####  3.3.4 The General Flow Problem

We now describe a general network flow problem that can be used to analyze a
large range of applications, from international trade to communication and
assignment. This general problem includes optimal transport as a special case.
Once we have introduced the problem, we show two results. First, the problem
can easily be formulated as a linear program and solved using standard linear
programming methods. Second, even though optimal transport is a strict subset
of the general flow problem, every general flow problem can be solved using a
combination of optimal transport and shortest path methods.

######  3.3.4.1 Problem Statement

We are interested in flow of a good or service across a network with $n$
vertices. This network can be understood as a weighted directed graph
$(V,E,c)$. To simplify notation, we label the nodes from $1$ to $n$ and let
$V=[n]$. Existence of an edge $e=(i,j)\in E$ with weight $c(i,j)$ indicates
that the good can be shipped from $i$ to $j$ at cost $c(i,j)$. We recall from
§1.4.1 that $\mathscr{I}(i)$ is the set of direct predecessors of vertex $i$
(all $u\in V$ such that $(u,i)\in E$) and $\mathscr{O}(i)$ is the set of
direct successors (all $j\in V$ such that $(i,j)\in E$). A classic example is
the famous Silk Road of antiquity, part of which is illustrated in Figure
3.15. Silk was produced in eastern cities such as Loyang and Changan, and then
transported westward to satisfy final demand in Rome, Constantinople and
Alexandria. Towns such as Yarkand acted as trade hubs. Rather than covering
the whole route, traders typically traveled backward and forward between one
pair of hubs, where they knew the language and customs.555Our use of the Silk
Road as an example of a network flow problem is not original. Galichon, (2018)
provides a highly readable treatment in the context of optimal transport.

Figure 3.15:  The Silk Road

Returning to the model, we allow for both initial supply of and final demand
for the good at every node (although one or both could be zero). Let $s(i)$
and $d(i)$ be supply and demand at node $i$ respectively. Aggregate supply and
demand over the network is assumed to be equal, so that

$\sum_{i\in V}s(i)=\sum_{i\in V}d(i).$ (3.32)

This can be understood as an equilibrium condition: prices have adjusted to
equalize initial supply and final demand in aggregate. We assume throughout
that the vectors $s$ and $d$ are nonnegative with at least one positive
element. Let $q(i,j)$ be the amount of the good shipped from node $i$ to node
$j$ for all $i,j\in V$. The minimum cost network flow problem is to minimize
total shipping cost

$\sum_{i\in V}\sum_{j\in V}c(i,j)q(i,j),$ (3.33)

subject to the restriction that $q\geqslant 0$ and

$s(i)+\sum_{v\in\mathscr{I}(i)}q(v,i)=d(i)+\sum_{j\in\mathscr{O}(i)}q(i,j)\quad\text{for
all }i\in V.$ (3.34)

The left hand side of (3.34) is total supply to node $i$ (initial supply plus
inflow from other nodes), while the right hand side is total demand (final
demand plus outflow to other nodes). Although we presented them separately,
the node-by-node restriction (3.34) implies the aggregate restriction (3.32).
Explain why this is the case. Since any outflow from some node $i$ is matched
by equal inflow into some node $j$, summing both sides of (3.34) across all
$i\in V$ yields (3.32).

######  3.3.4.2 Optimality

There are several ways to transform the network flow problem into a linear
program. We follow the presentation in Bertsimas and Tsitsiklis, (1997). We
take $m=|E|$ to be the total number of edges and enumerate them (in any
convenient way) as $e_{1},\ldots,e_{m}$. Let’s say that $e_{k}$ leaves node
$i$ if $e_{k}=(i,j)$ for some $j\in[n]$, and that $e_{k}$ enters node $i$ if
$e_{k}=(\ell,i)$ for some $\ell\in[n]$. Then we define the $n\times m$ node-
edge incidence matrix $A$ by

$A=(a_{ik})\quad\text{with}\quad a_{ik}:=\begin{cases}\;1&\text{if $e_{k}$
leaves $i$}\\\ -1&\text{if $e_{k}$ enters $i$}\\\
\;0&\text{otherwise}.\end{cases}$

######  Example 3.3.1.

Consider the very simple minimum cost flow problem in Figure 3.16. The
shipment costs $c(i,j)$ are listed next to each existing edge. Initial supply
is 10 at node 1 and zero elsewhere. Final demand is 10 at node 4 and zero
elsewhere. We enumerate the edges as

$E=\\{e_{1},\ldots,e_{4}\\}=\\{(1,2),(1,4),(2,3),(3,4)\\}.$ (3.35)

The node-edge incidence matrix is

$A=\begin{pmatrix}1&1&0&0\\\ -1&0&1&0\\\ 0&0&-1&1\\\ 0&-1&0&-1\end{pmatrix}.$

$1$$s_{1}=10$$2$$3$$4$$d_{4}=10$$4$$1$$1$$1$ Figure 3.16:  A simple network
flow problem

Now, returning to the general case, we rearrange $q$ and $c$ into $m\times 1$
vectors $(q_{k})$ and $(c_{k})$, where $q_{k}$ is the amount shipped along
edge $k$ and $c_{k}$ is the cost. For example, if $e_{k}=(i,j)$, then
$q_{k}=q(i,j)$ and $c_{k}=c(i,j)$. In addition, we set $b$ to be the vector in
$\mathbbm{R}^{n}$ with $i$-th element $s(i)-d(i)$, which is net exogenous
supply at node $i$. In this set up, show that (3.34) is equivalent to $Aq=b$
in the special case of Example 3.3.1. Let $(Aq)(i)$ be the $i$-th row of $Aq$.
Show that

$(Aq)(i)=\sum_{j\in\mathscr{O}(i)}q(i,j)-\sum_{v\in\mathscr{I}(i)}q(v,i).$
(3.36)

Verifying the claim is just a matter of working with the definition of $A$.
Fixing $i\in[n]$, we have

$(Aq)(i)=\sum_{k=1}^{m}a_{ik}q_{k}=\sum_{k=1}^{m}\mathbbm{1}\\{e_{k}\text{
leaves }i\\}q_{k}-\sum_{k=1}^{m}\mathbbm{1}\\{e_{k}\text{ points to
}i\\}q_{k},$

This is equal to
$\sum_{j\in\mathscr{O}(i)}q(i,j)-\sum_{v\in\mathscr{I}(i)}q(v,i)$, as was to
be shown. Equation (3.36) tells us that the $i$-th row of $Aq$ give us the net
outflow from node $i$ under the transport plan $q$. Now, with
$\left\langle{c,q}\right\rangle:=\sum_{k=1}^{m}c_{k}q_{k}$, the minimum cost
network flow problem can now be expressed as

$\min\left\langle{c,q}\right\rangle\quad\ \mathrm{s.t.}\ \quad q\geqslant
0\;\text{ and }\;Aq=b.$ (3.37)

This is a linear program in standard equality form, to which we can apply any
linear programming solver. For Example 3.3.1, we run the following:

⬇

A = (( 1, 1, 0, 0),

(-1, 0, 1, 0),

( 0, 0, -1, 1),

( 0, -1, 0, -1))

\parb = (10, 0, 0, -10)

c = (1, 4, 1, 1)

\parresult = linprog(c, A_eq=A, b_eq=b, method=’highs-ipm’)

print(result.x)

The output is [10. 0. 10. 10.]. Recalling the order of the paths in (3.35),
this means that the optimal transport plan is $q(1,4)=0$ and
$q(1,2)=q(2,3)=q(3,4)=10$, as our intuition suggests. Some network flow
problems have capacity constraints, which can be modeled as a map $g\colon
E\to[0,\infty]$, along with the restriction $q(e)\leqslant g(e)$ for all $e\in
E$. (If $g(e)=+\infty$, there is no capacity constraint over shipping on edge
$e$.) Formulate this as a linear program and modify the code above, which
solves Example 3.3.1, to include the capacity constraint $g(1,2)=5$. Solve for
the optimal plan. To the code that solves the original version of Example
3.3.1, we need to add

⬇

bounds = ((0, 5),

(0, None),

(0, None),

(0, None))

and then change the function call to

⬇

result = linprog(c, A_eq=A, b_eq=b, method=’highs-ipm’, bounds=bounds)

print(result.x)

The output is [5. 5. 5. 5.], which also agrees with our intuition. Explain how
the generic optimal transport problem treated in §3.3.1 is a special case of
the minimum cost network flow problem.

######  3.3.4.3 Reduction to Optimal Transport

In Exercise 3.3.4.2, we saw how every optimal transport problem is a special
kind of minimum cost network flow problem. There is a sense in which the
converse is also true. In particular, we can use optimal transport methods to
solve any network flow problem, provided that we first modify the network flow
problem via an application of shortest paths. To explain how this works, we
take the abstract network flow problem described in §3.3.4.1, on the weighted
digraph $(V,E,c)$, with $V=[n]$, initial supply vector
$s\in\mathbbm{R}^{n}_{+}$ and final demand vector $d\in\mathbbm{R}^{n}_{+}$.
For the purposes of this section, we agree to call a node $i$ with
$s(i)-d(i)>0$ a net supplier. A node $i$ with $d(i)-s(i)>0$ will be called a
net consumer. Nodes with $s(i)=d(i)$ will be called trading stations.

######  Example 3.3.2.

In the left hand side of Figure 3.17, nodes 1 and 2 are net suppliers, 3 is a
trading station and 4 and 5 are net consumers.

######  Example 3.3.3.

In the Silk Road application, Rome would be a net consumer, where final demand
is large and positive, while initial supply is zero. A city such as Yarkand
should probably be modeled as a trading station, with $s(i)=d(i)=0$.

The idea behind the reduction is to treat the net supplier nodes as source
locations and the net consumer nodes as target locations in an optimal
transport problem. The next step is to compute the shortest path (if there are
multiple, pick any one) from each net supplier $i$ to each net consumer $j$.
Let $\rho(i,j)$ denote this path, represented as a sequence of edges in $E$.
The cost of traversing $\rho(i,j)$ is

$\hat{c}(i,j):=\sum_{k=1}^{m}c_{k}\mathbbm{1}\\{e_{k}\in\rho(i,j)\\}.$

Now the trading stations are eliminated and we solve the optimal transport
problem with

  * •

$\mathsf{X}=$ the set of net suppliers,

  * •

$\mathsf{Y}=$ the set of net consumers,

  * •

$\varphi(i)=s(i)-d(i)$ on $\mathsf{X}$,

  * •

$\psi(j)=d(j)-s(j)$ on $\mathsf{Y}$, and

  * •

cost function $\hat{c}(i,j)$ as defined above.666If no path exists from $i$ to
$j$ then we set $\hat{c}(i,j)=\infty$. Such settings can be handled in linear
programming solvers by adding capacity constraints. See, for example, Peyré et
al., (2019), Section 10.3.

After we find the optimal transport plan $\pi$, the network minimum cost flow
$q_{k}$ along arbitrary edge $e_{k}\in E$ is recovered by setting

$q_{k}=\sum_{i\in\mathsf{X}}\sum_{j\in\mathsf{Y}}\pi(i,j)\mathbbm{1}\\{e_{k}\in\text{
the shortest path from $i$ to $j$}\\}.$

Recalling our assumptions on $s$ and $d$, prove that
$\sum_{i\in\mathsf{X}}\varphi(i)=\sum_{j\in\mathsf{Y}}\psi(j)$ and that this
sum is nonzero.

######  Remark 3.3.2.

We have not imposed $\sum_{i}\varphi(i)=\sum_{j}\psi(j)=1$, as required for
the standard formulation of the optimal transport problem. But this
normalization is only used for convenience in exposition and most solvers do
not require it.

Figure 3.17 illustrates the method. Trading station 3 is eliminated after the
shortest paths are computed.

$1$$2$$3$$4$$5$$1$$2$$4$$5$$1$$5$$1$$1$$1$$3$$3$$2$$2$ Figure 3.17:  Reducing
minimum cost optimal flow to optimal transport

####  3.4 Chapter Notes

Our treatment of shortest paths can be understood as a simplified version of
both the Bellman–Ford algorithm and Dijkstra’s algorithm, which are routinely
used to solve large shortest path planning problems. Our approach is intended
to emphasize recursive solution methods, which are valuable for analyzing a
vast range of economic problems, from intertemporal modeling (see, e.g., Lucas
and Stokey, (1989) or Ljungqvist and Sargent, (2018)) to production chains
(Kikuchi et al.,, 2021). For a more in-depth treatment of linear programming,
we recommend the excellent textbooks by Bertsimas and Tsitsiklis, (1997) and
Matousek and Gärtner, (2007). For bedtime reading, Cook, (2011) provides an
entertaining introduction to some of the main ideas and applications to
network problems, including a review of computation. While §3.2.1 provided a
rather whimsical introduction to matching and assignment problems, such
problems have great real world importance. Examples include assigning kidney
donors to recipients, mothers to maternity wards, doctors to hospitals,
students to schools, delivery drivers to orders and autonomous vehicles to
riders.777The latter problem occurs in the field of Autonomous Mobility on
Demand (AMoD). See, for example, Ascenzi and Palanza, (2021) or Simonetto et
al., (2019). A brief history of assignment, matching problems and
combinatorial optimization can be found in Schrijver, (2005). Greinecker and
Kah, (2021) study existence of stable matchings in a setting with many agents.
Villani, (2008) and Vershik, (2013) provide extensive historical background on
the optimal transport problem. Vershik, (2013) mentions some of the problems
that Kantorovich faced, as a Soviet mathematician working in the time of
Stalin and Khrushchev, given that his main duality theorem for optimal
transport can be seen as a proof that competitive market equilibria attain the
maximal transport plan. In §3.3.2, we mentioned that Kantorovich’s work
anticipated much of the later development of duality theory for arbitrary
linear programs. In fact, according to Vershik, (2013), Kantorovich
anticipated much of the general theory of linear programming itself, including
providing a version of the simplex algorithm later rediscovered and extended
by Dantzig. Optimal transport has a remarkably wide variety of applications,
spread across economics, econometrics, finance, statistics, artificial
intelligence, machine learning and other fields. Within economics, Galichon,
(2018) provides an excellent overview. Fajgelbaum and Schaal, (2020) consider
optimal transport in spacial equilibria. Beiglböck et al., (2022) review some
of the major milestones of modern finance theory and show their connections
via optimal transport. Connections to machine learning are surveyed in Kolouri
et al., (2017). The computational theory of optimal transport is now a major
field. A high quality exposition can be found in Peyré et al., (2019).
Blanchet et al., (2018) use computational optimal transport to solve for
Cournot–Nash equilibria in mean-field type games.

###  Chapter 4 Markov Chains and Networks

Markov chains evolving on finite sets are a foundational class of stochastic
processes, routinely employed in quantitative modeling within economics,
finance, operations research and social science. A Markov chain is most easily
understood as a weighted digraph, with graph-theoretic properties such as
connectedness and periodicity being key determinants of dynamics.

####  4.1 Markov Chains as Digraphs

We begin with fundamental definitions and then investigate dynamics.

#####  4.1.1 Markov Models

A finite Markov model is a weighted directed graph $\mathscr{M}=(S,E,p)$,
where $S$ is the (finite) set of vertices, $E$ is the set of edges and $p$ is
the weight function, with the additional restriction that

$\sum_{y\,\in\,\mathscr{O}(x)}p(x,y)=1\quad\text{for all }x\in S.$ (4.1)

Figure 1.15 on page 1.15 presented an example of such a digraph. The set of
vertices $S$ of a finite Markov model $\mathscr{M}=(S,E,p)$ is also called the
state space of the model, and vertices are called states. The two standard
interpretations are

  1. (i)

$S$ is a set of possible states for some random element (the state) and the
weight $p(x,y)$ represents the probability that the state moves from $x$ to
$y$ in one step.

  2. (ii)

$S$ is a set of possible values for some measurement over a large population
(e.g., hours worked per week measured across a large cross-section of
households) and $p(x,y)$ is the fraction of agents that transition from state
$x$ to state $y$ in one unit of time.

These two perspectives are related in ways that we explore below.

######  4.1.1.1 Transition Matrices

If $\mathscr{M}$ is a finite Markov model, then the restriction (4.1) is
equivalent to the statement that the adjacency matrix associated with
$\mathscr{M}$ is stochastic (see §1.3.1.3 for the definition). Identifying
$S=\\{\text{poor},\text{middle},\text{rich}\\}$ with $\\{1,2,3\\}$, the
adjacency matrix for this weighted digraph is

$P_{a}=\begin{pmatrix}0.9&0.1&0.0\\\ 0.4&0.4&0.2\\\ 0.1&0.1&0.8\end{pmatrix}$
(4.2)

Since $P_{a}\geqslant 0$ and rows sum to unity, $P_{a}$ is stochastic as
required. In the context of finite Markov models, the adjacency matrix of
$\mathscr{M}$ is also called the transition matrix. Regarding notation, when
$S$ has typical elements $x,y$, it turns out to be convenient to write
elements of the transition matrix $P$ as $P(x,y)$ rather than $P_{ij}$ or
similar. We can think of $P$ as extending the weight function $p$ from $E$ to
the set of all $(x,y)$ pairs in $S\times S$, assigning zero whenever
$(x,y)\notin E$. As such, for every possible choice of $(x,y)$, the value
$P(x,y)$ represents the probability of transitioning from $x$ to $y$ in one
step. The requirement that $P$ is stochastic can now be written as $P\geqslant
0$ and

$\sum_{y\in S}P(x,y)=1\text{ for all }x\in S.$ (4.3)

The restriction in (4.3) just says that the state space is “complete:” after
arriving at $x\in S$, the state must now move to some $y\in S$. Using notation
from §1.3.1, to say that $P$ is stochastic is the same as requiring that each
row of $P$ is in $\mathscr{D}(S)$.

######  Example 4.1.1.

A Markov model is estimated in the international growth dynamics study of
Quah, (1993). The state is real GDP per capita in a given country relative to
the world average. Quah discretizes the possible values to 0–1/4, 1/4–1/2,
1/2–1, 1–2 and 2–$\infty$, calling these states 1 to 5 respectively. The
transitions are over a one year period. Estimated one step transition
probabilities are represented as a weighted digraph in Figure 4.1, where

  * •

$S=\\{1,\ldots,5\\}$ is the state space

  * •

the set of edges $E$ is represented by arrows and

  * •

transition probabilities are identified with weights attached to these edges.

12345$0.03$$0.05$$0.03$$0.04$$0.04$$0.04$$0.02$$0.01$$0.97$$0.92$$0.92$$0.94$$0.99$

Figure 4.1:  Cross-country GDP dynamics as a digraph

The transition matrix for the Markov model in Example 4.1.1 is

$P_{Q}=\left(\begin{array}[]{ccccc}0.97&0.03&0.00&0.00&0.00\\\
0.05&0.92&0.03&0.00&0.00\\\ 0.00&0.04&0.92&0.04&0.00\\\
0.00&0.00&0.04&0.94&0.02\\\ 0.00&0.00&0.00&0.01&0.99\end{array}\right)$ (4.4)

Note the large values on the principal diagonal of $P_{Q}$. These indicate
strong persistence: the state stays constant from period to period with high
probability. Quah, (1993) estimated $P_{Q}$ by maximum likelihood, pooling
transitions over the years 1960–1984. (In this case maximum likelihood
estimation amounts to recording the relative frequency of transitions between
states.) Figure 4.2 shows how the numbers change if we repeat the exercise
using World Bank GDP data from 1985–2019. The numbers are quite stable
relative to the earlier estimate. Below we will examine how long run
predictions are affected.

12345$0.03$$0.05$$0.06$$0.05$$0.05$$0.06$$0.04$$0.01$$0.97$$0.89$$0.90$$0.90$$0.99$

Figure 4.2:  Cross-country GDP dynamics as a digraph, updated data

As another example, Benhabib et al., (2019) estimate the following transition
matrix for intergenerational social mobility:

$P_{B}:=\begin{pmatrix}0.222&0.222&0.215&0.187&0.081&0.038&0.029&0.006\\\
0.221&0.22&0.215&0.188&0.082&0.039&0.029&0.006\\\
0.207&0.209&0.21&0.194&0.09&0.046&0.036&0.008\\\
0.198&0.201&0.207&0.198&0.095&0.052&0.04&0.009\\\
0.175&0.178&0.197&0.207&0.11&0.067&0.054&0.012\\\
0.182&0.184&0.2&0.205&0.106&0.062&0.05&0.011\\\
0.123&0.125&0.166&0.216&0.141&0.114&0.094&0.021\\\
0.084&0.084&0.142&0.228&0.17&0.143&0.121&0.028\end{pmatrix}$ (4.5)

Here the states are percentiles of the wealth distribution. In particular,
with the states represented by $1,2,\ldots,8$, the corresponding percentiles
are

0–20%, 20–40%, 40–60%, 60–80%, 80–90%, 90–95%, 95–99%, 99–100%

Transition probabilities are estimated from US 2007–2009 Survey of Consumer
Finances data. Relative to the highly persistent matrix $P_{Q}$, less weight
on the principle diagonal suggests more mixing—the influence of initial
conditions is relatively short-lived. Additional insight about the dynamics
can be obtained from a contour plot of the matrix $P_{B}$, as in Figure 4.3.
Here $P_{B}$ is plotted as a heat map after rotating it 90 degrees
anticlockwise. The rotation is so that the dynamics are comparable to the 45
degree diagrams often used to understand discrete time dynamic systems. A
vertical line from state $x$ corresponds to the next period conditional
distribution $P(x,\cdot)$. In this case, we see that, for example, lower
states are quite persistent, whereas households in the highest state tend to
fall back towards the middle.

Figure 4.3:  Contour plot of transition matrix $P_{B}$

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Show that $U\subset S$ is absorbing (see §1.4.1.3) for the digraph
$\mathscr{M}$ if and only if

$\sum_{y\in U}P(x,y)=1\quad\text{for all}\quad x\in U.$ (4.6)

######  4.1.1.2 Markov Chains

Consider a finite Markov model $\mathscr{M}$ with state space $S$ and
transition matrix $P$. As before, $P(x,y)$ indicates the probability of
transitioning from $x$ to $y$ in one step. Another way to say this is that,
when in state $x$, we update to a new state by choosing it from $S$ via the
distribution $P(x,\cdot)$. The resulting stochastic process is called a Markov
chain. We can state this more formally as follows. Let
$(X_{t})_{t\in\mathbbm{Z}_{+}}$ be a sequence of random variables taking
values in $S$. We say that $(X_{t})$ is a Markov chain on $S$ if there exists
a stochastic matrix $P$ on $S$ such that

$\mathbbm{P}\\{X_{t+1}=y\mid
X_{0},X_{1},\ldots,X_{t}\\}=P(X_{t},y)\quad\text{for all}\quad t\geqslant
0,\;y\in S.$ (4.7)

To simplify terminology, we also call $(X_{t})$ $P$-Markov when it satisfies
(4.7) . We call either $X_{0}$ or its distribution $\psi_{0}$ the initial
condition of $(X_{t})$ depending on context. The definition of a Markov chain
says two things:

  1. (i)

When updating to $X_{t+1}$ from $X_{t}$, earlier states are not required.

  2. (ii)

The matrix $P$ encodes all of the information required to perform the update,
given $X_{t}$.

One way to think about Markov chains is algorithmically: Let $P$ be a
stochastic matrix and let $\psi_{0}$ be an element of $\mathscr{D}(S)$. Now
generate $(X_{t})$ via Algorithm 3. The resulting sequence is $P$-Markov with
initial condition $\psi_{0}$.

1set $t=0$ and draw $X_{t}$ from $\psi_{0}$

2 while  _$t <\infty$_ do draw $X_{t+1}$ from the distribution
$P(X_{t},\cdot)$

3 let $t=t+1$

4 end while

Algorithm 3 Generation of $P$-Markov $(X_{t})$ with initial condition
$\psi_{0}$

######  4.1.1.3 Simulation

For both simulation and theory, it is useful to be able to translate Algorithm
3 into a stochastic difference equation governing the evolution of
$(X_{t})_{t\geqslant 0}$. We now outline the procedure, which uses inverse
transform sampling (see §1.3.1.2). For simplicity, we assume that $S=[n]$,
with typical elements $i,j$. The basic idea is to apply the inverse transform
method to each row of $P$ and then sample by drawing a uniform random variable
at each update. To this end, we set

$F(i,u):=\sum_{j=1}^{n}j\mathbbm{1}\\{q(i,j-1)<u\leqslant q(i,j)\\}\qquad(i\in
S,\;u\in(0,1)),$

where, for each $i,j\in S$, the value $q(i,j)$ is defined recursively by

$q(i,j):=q(i,j-1)+P(i,j)\quad\text{with }q(i,0)=0.$

Let $U(0,1)$ represent the uniform distribution on $(0,1)$ and take

$X_{t+1}=F(X_{t},U_{t+1})\quad\text{where
}\;(U_{t})\stackrel{{\scriptstyle\textrm{ {\sc iid }}}}{{\sim}}U(0,1).$ (4.8)

If $X_{0}$ is an independently drawn random variable with distribution
$\psi_{0}$ on $S$, then $(X_{t})$ is $P$-Markov on $S$ with initial condition
$\psi_{0}$, as Exercise 4.1.1.3 asks you to show. Conditional on $X_{t}=i$,
show that, for given $j\in S$,

  1. (i)

$X_{t+1}=j$ if and only if $U_{t+1}$ lies in the interval $(q(i,j-1),q(i,j)]$.

  2. (ii)

This event has probability $P(i,j)$.

Conclude that $X_{t+1}$ in (4.8) is a draw from $P(i,\cdot)$. Point (i) is
immediate from the definition of $F$. Regarding (ii), from
$q(i,j)=q(i,j-1)+P(i,j)$ we have $P(i,j)=q(i,j)-q(i,j-1)$, which is the length
of the interval $(q(i,j-1),q(i,j)]$. The probability that $U_{t+1}$ falls in
this interval is its length, which we just agreed is $P(i,j)$. The proof is
now complete. Each subfigure in Figure 4.4 shows realizations of two Markov
chains, both generated using the stochastic difference equation (4.8). The
sequences are generated each with its own independent sequence of draws
$(U_{t})$. The underlying transition matrices are $P_{B}$ from (4.5) in the
top panel and $P_{Q}$ from (4.4) in the bottom panel. In both panels, one
chain starts from the lowest state and the other from the highest. Notice that
time series generated by $P_{B}$ mix faster than those generated by $P_{Q}$:
the difference in initial states is not a strong predictor of outcomes after
an initial “burn in” period. We discuss mixing and its connection to stability
below.

Figure 4.4:  Wealth percentile over time

######  4.1.1.4 Higher Order Transition Matrices

Given a finite Markov model $\mathscr{M}$ with state space $S$ and transition
matrix $P$, define $(P^{k})_{k\in\mathbbm{N}}$ by $P^{k+1}=PP^{k}$ for all
$k$, with the understanding that $P^{0}=I=$ the identity matrix. In other
words, for each $k$, the matrix $P^{k}$ is the $k$-th power of $P$. If we
spell out the matrix product $P^{k+1}=PP^{k}$ element-by-element, we get

$P^{k+1}(x,y):=\sum_{z}P(x,z)P^{k}(z,y)\qquad(x,y\in S,\;k\in\mathbbm{N}).$
(4.9)

Prove that $P^{k}$ is a stochastic matrix on $S$ for all $k\in\mathbbm{N}$.
$P^{k}$ is stochastic for all $k\in\mathbbm{N}$ by induction and the fact that
the set of stochastic matrices is closed under multiplication (see 1.3.1.3).
In this context, $P^{k}$ is called the $k$-step transition matrix
corresponding to $P$. The $k$-step transition matrix has the following
interpretation: If $(X_{t})$ is $P$-Markov, then, for any $t,k\in\mathbbm{N}$
and $x,y\in S$,

$P^{k}(x,y)=\mathbbm{P}\\{X_{t+k}=y\,|\,X_{t}=x\\}.$ (4.10)

In other words, $P^{k}$ provides the $k$-step transition probabilities for the
$P$-Markov chain $(X_{t})$, as suggested by its name. This claim can be
verified by induction. Fix $t\in\mathbbm{N}$ and $x,y\in S$. The claim is true
by definition when $k=1$. Suppose the claim is also true at $k$ and now
consider the case $k+1$. By the law of total probability, we have

$\mathbbm{P}\\{X_{t+k+1}=y\,|\,X_{t}=x\\}=\sum_{z}\mathbbm{P}\\{X_{t+k+1}=y\,|\,X_{t+k}=z\\}\mathbbm{P}\\{X_{t+k}=z\,|\,X_{t}=x\\}.$

The induction hypothesis allows us to use (4.10), so the last equation becomes

$\mathbbm{P}\\{X_{t+k+1}=y\,|\,X_{t}=x\\}=\sum_{z}P(z,y)P^{k}(x,z)=P^{k+1}(x,y).$

This completes our proof by induction. A useful identity for the higher order
Markov matrices is

$P^{j+k}(x,y)=\sum_{z}P^{k}(x,z)P^{j}(z,y)\qquad((x,y)\in S\times S)$ (4.11)

which holds for any $j,k$ in $\mathbbm{N}$. This is called the
Chapman–Kolmogorov equation. Note that

  * •

(4.9) is a special case of (4.11) and

  * •

(4.11) is a special case of (1.26) on page 1.26, written with different
notation.

To provide probabilistic intuition for the validity of the Chapman–Kolmogorov
equation, let $X_{0}=x$ and let $y\in S$ be given. Using the law of total
probability again, we have

$\mathbbm{P}\\{X_{j+k}=y\,|\,X_{0}=x\\}=\sum_{z}\mathbbm{P}\\{X_{j+k}=y\,|\,X_{0}=x,\,X_{k}=z\\}\mathbbm{P}\\{X_{k}=z\,|\,X_{0}=x\\}$

By Markov property (4.7), the future and past are independent given the
present, so

$\sum_{z}\mathbbm{P}\\{X_{j+k}=y\,|\,X_{0}=x,\,X_{k}=z\\}=\sum_{z}\mathbbm{P}\\{X_{j+k}=y\,|\,X_{k}=z\\}.$

As a result of this fact and (4.10), the equation before this one can be
rewritten as (4.11).

#####  4.1.2 Distribution Dynamics

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Let $(X_{t})$ be $P$-Markov and, for each $t\geqslant 0$, let
$\psi_{t}\in\mathscr{D}(S)$ be defined by

$\psi_{t}:=\mathbbm{P}\\{X_{t}=\cdot\\}=\text{ the distribution of $X_{t}$}.$

The vector $\psi_{t}$ is called the marginal distribution of $X_{t}$. While
$(X_{t})$ is random, the sequence $(\psi_{t})$ is deterministic. In this
section we investigate its dynamics.

######  4.1.2.1 Updating Marginal Distributions

The key idea for this section is that there is a simple link between
successive marginal distributions: by the law of total probability, we have

$\mathbbm{P}\\{X_{t+1}=y\\}=\sum_{x}\mathbbm{P}\\{X_{t+1}=y\,|\,X_{t}=x\\}\cdot\mathbbm{P}\\{X_{t}=x\\},$

which can be rewritten as

$\psi_{t+1}(y)=\sum_{x}P(x,y)\psi_{t}(x)\qquad(y\in S).$ (4.12)

This fundamental expression tells how to update marginal distributions given
the transition matrix $P$. When each $\psi_{t}$ is interpreted as a row
vector, we can write (4.12) as

$\psi_{t+1}=\psi_{t}P.$ (4.13)

(Henceforth, in expressions involving matrix algebra, distributions are row
vectors unless otherwise stated). Some authors refer to (4.13) as the forward
equation associated with $P$, while $\psi\mapsto\psi P$ is called the forward
operator, by analogy with the Kolmogorov forward equation from continuous
time. Think of (4.13) as a difference equation in distribution space.
Iterating backwards,

$\psi_{t}=\psi_{t-1}P=(\psi_{t-2}P)P=\psi_{t-2}P^{2}=\cdots=\psi_{0}P^{t}.$

Given any $\psi_{0}\in\mathscr{D}(S)$ and $t\in\mathbbm{N}$, we have

$\psi_{0}P^{t}=\text{ the distribution of $X_{t}$ given
$X_{0}\stackrel{{\scriptstyle d}}{{=}}\psi_{0}$}.$

######  Example 4.1.2.

As an illustration, let’s take the matrix $P_{Q}$ estimated by Quah, (1993)
using 1960–1984 data and use $P_{Q}^{t}$ to update the 1985 distribution
$t=2019-1985=34$ times, in order to obtain a prediction for the cross-country
income distribution in 2019. Figure 4.5 shows how this prediction fares
compared to the realized 2019 distribution, calculated from World Bank GDP
data.111Example 4.1.2 is intended as an illustration of the mechanics of
updating distributions. While the methodology in Quah, (1993) is thought
provoking, it struggles to make plausible long run predictions about something
as complex as the cross-country income distribution. Indeed, the
Kullback–Leibler deviation between the predicted and realized 2019
distributions is actually larger that the Kullback–Leibler deviation between
the 1985 and 2019 distributions. Evidently, a naive estimate “nothing will
change” model predicts better than Quah’s. Most of the focus in this text is
on short term forecasts and scenarios, where the system is approximately
stationary after suitable transformations, rather than heroic long term
predictions where nonstationary change is hard to quantify.

Figure 4.5:  Predicted vs realized cross-country income distributions for 2019

######  4.1.2.2 Trajectories in the Long Run

One of the main sub-fields of Markov chain analysis is asymptotics of
distribution sequences. This topic turns out to be important for network
theory as well. In §4.2.2 we will investigate asymptotics in depth. In this
section we build some intuition via simulations. Figure 4.6 shows the
trajectory $(\psi P_{a}^{t})$ when $S=\\{1,2,3\\}$, $\psi=(0,0,1)$, and
$P_{a}$ is the transition matrix displayed in (4.2). The blue triangle is the
unit simplex in $\mathbbm{R}^{3}_{+}$, consisting of all row vectors
$\psi\in\mathbbm{R}^{3}$ such that $\psi\geqslant 0$ and $\psi\mathbbm{1}=1$.
The unit simplex can be identified with $\mathscr{D}(S)$ when $S=\\{1,2,3\\}$.
Red dots in the figure are distributions in the sequence $(\psi P_{a}^{t})$
for $t=0,\ldots,20$. Figure 4.7 shows distribution dynamics for $P_{a}$ that
start from initial condition $\psi=(0,1/2,1/2)$. It seems that both of the
sequences are converging. This turns out to be the case—the black dot in the
figures is the limit of both sequences and, moreover, this point is a
stationary distribution for $P_{a}$, as defined in §1.3.1.3. In fact we
can—and will—also show that no other stationary distribution exists, and that
$\psi P_{a}^{t}$ converges to the stationary distribution regardless of the
choice of $\psi\in\mathscr{D}(S)$. This is due to certain properties of
$P_{a}$, related to connectivity and aperiodicity.

Figure 4.6:  A trajectory from $\psi_{0}=(0,0,1)$

Figure 4.7:  A trajectory from $\psi_{0}=(0,1/2,1/2)$

Figure 4.8 provides another view of a distribution sequence, this time
generated from the matrix $P_{B}$ in (4.5). The initial condition $\psi_{0}$
was uniform. Each distribution shown was calculated as $\psi P_{Q}^{t}$ for
different values of $t$. The distribution of across wealth classes converges
rapidly to what appears to be a long run limit. Below we confirm that this is
so, and that the limit is independent of the initial condition. The rapid rate
of convergence is due to the high level of mixing in the transition dynamics.

Figure 4.8:  Distribution projections from $P_{B}$

#####  4.1.3 Stationarity

In this section we examine stationary distributions and their properties. As
we will see, stationary distributions can be regarded as steady states for the
evolution of the sequence of marginal distributions. (Later, when we study
ergodicity, stationary distributions will aquire another important
interpretation.)

######  4.1.3.1 Stationary Distributions

Recall from §1.3.1.3 that if $P$ is a stochastic matrix and
$\psi\in\mathbbm{R}^{n}_{+}$ is a row vector such that $\psi\mathbbm{1}=1$ and
$\psi P=\psi$, then $\psi$ is called stationary for $P$. So now let
$\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Translating to the notation of Markov chains, a distribution
$\psi^{*}\in\mathscr{D}(S)$ is stationary for $\mathscr{M}$ if

$\psi^{*}(y)=\sum_{x\in S}P(x,y)\psi^{*}(x)\quad\text{for all }y\in S.$

We also write this expression as $\psi^{*}=\psi^{*}P$ and understand
$\psi^{*}$ as a fixed point of the map $\psi\mapsto\psi P$ that updates
distributions (cf., Equation (4.13)). If $X_{t}\stackrel{{\scriptstyle
d}}{{=}}\psi^{*}$, then, for any $j\in\mathbbm{N}$, the fixed point property
implies that $X_{t+j}\stackrel{{\scriptstyle d}}{{=}}\psi^{*}P^{j}=\psi^{*}$.
Hence

$X_{t}\stackrel{{\scriptstyle d}}{{=}}\psi^{*}\implies
X_{t+j}\stackrel{{\scriptstyle d}}{{=}}\psi^{*}\;\text{ for all
}j\in\mathbbm{N}.$

######  Example 4.1.3.

Suppose workers are hired at rate $\alpha$ and fired at rate $\beta$,
transitioning between unemployment and employment according to

$P_{w}=\begin{pmatrix}1-\alpha&\alpha\\\ \beta&1-\beta\end{pmatrix}.$ (4.14)

In §1.2.3.3 we showed $\alpha+\beta>0$ implies
$\psi=(\beta/(\alpha+\beta),\alpha/(\alpha+\beta))$ is a dominant left
eigenvector for $P_{w}$. Since $r(P_{w})=1$, the Perron–Frobenius Theorem
tells us that $\psi$ is stationary for $P_{w}$. One implication is that, if
the distribution of workers across the two employment states is given by
$\psi$ and updates obey the dynamics encoded in $P_{w}$, then no further
change is observed in the unemployment rate.

######  Example 4.1.4.

The black dot in each of Figures 4.6–4.7, indicating the limit of the marginal
sequences, is a stationary distribution for the Markov matrix $P_{a}$
displayed in (4.2), page 4.2. We discuss its calculation in §4.1.3.4.

Let $\mathscr{M}$ be a finite Markov model with state space $|S|=n$ and
transition matrix $P$. Let $\psi\equiv 1/n$ be the uniform distribution on
$S$. Prove that $\psi$ is stationary for $P$ if and only if $P$ is doubly
stochastic (i.e., has unit column sums as well as unit row sums). In Figures
4.6–4.7 all trajectories converge towards the stationary distribution. Not all
stationary distributions have this “attractor” property, and in general there
can be many stationary distributions. The next example illustrates. Let
$\mathscr{M}=(S,E,p)$ be a finite Markov model with $(x,y)\in E$ if and only
if $x=y$. Describe the implied weight function and corresponding transition
matrix. Show that every distribution in $\mathscr{D}(S)$ is stationary for
$\mathscr{M}$. Let $\mathscr{M}$ be as described. In this setting, the
requirement that the transition matrix $P$ is stochastic implies that
$p(x,x)=1$ and $P(x,y)=\mathbbm{1}\\{x=y\\}$. Thus, $P=I$, the $n\times n$
identity matrix. Every distribution is stationary because $\psi I=\psi$ for
all $\psi\in\mathscr{D}(S)$.

######  4.1.3.2 Existence and Uniqueness

From the Perron–Frobenius Theorem we easily obtain the following fundamental
result.

######  Theorem 4.1.1 (Existence and Uniqueness of Stationary Distributions).

Every finite Markov model $\mathscr{M}=(S,E,p)$ has at least one stationary
distribution $\psi^{*}$ in $\mathscr{D}(S)$. If the digraph $(S,E)$ is
strongly connected, then $\psi^{*}$ is unique and everywhere positive on $S$.

###### Proof.

Let $\mathscr{M}$ be a finite Markov model. Since the corresponding adjacency
matrix $P$ is stochastic, existence follows from Exercise 1.3.1.3 in §1.3.1.3.
By Theorem 1.4.3 on page 1.4.3, strong connectedness of $\mathscr{M}$ implies
irreducibility of $P$. When $P$ is irreducible, uniqueness and everywhere
positivity of the stationary distribution follows from the Perron–Frobenius
Theorem (page 1.2.6). ∎

The basic idea behind the uniqueness part of Theorem 4.1.1 is as follows:
Suppose to the contrary that $\mathscr{M}$ is strongly connected and yet two
distinct stationary distributions $\psi^{*}$ and $\psi^{**}$ exist in
$\mathscr{P}(S)$. Since a $P$-Markov chain started at $\psi^{*}$ always has
marginal distribution $\psi^{*}$ and likewise for $\psi^{**}$, different
initial conditions lead to different long run outcomes. This contradicts
strong connectedness in the following way: strong connectedness implies that
both chains traverse the whole state space. Moreover, being Markov chains,
they forget the past once they arrive at any state, so the starting draws
should eventually be irrelevant. (Actually, the story in the last paragraph is
incomplete. Initial conditions _can_ determine long run outcomes under strong
connectedness in one sense: a “periodic” Markov model can traverse the whole
space but _only at specific times_ that depend on the starting location. If we
rule out such periodicity, we get results that are even stronger than Theorem
4.1.1. This topic is treated in §4.2.2.)

######  4.1.3.3 Brouwer’s Theorem

Another way to prove the existence claim in Theorem 4.1.1 is via the famous
fixed point theorem of L. E. J. Brouwer (1881–1966).

######  Theorem 4.1.2 (Brouwer).

If $C$ is a convex compact subset of $\mathbbm{R}^{n}$ and $G$ is a continuous
self-map on $C$, then $G$ has at least one fixed point in $C$.

The proof of Theorem 4.1.2 in one dimension is not difficult, while the proof
in higher dimensions is challenging. See, for example, Aliprantis and Border,
(1999). Prove Brouwer’s fixed point theorem for the case $C=[0,1]$ by applying
the intermediate value theorem. Let $G$ be a continuous function from $[0,1]$
to itself and set $f(x):=Gx-x$. Since $G$ is a self-map on $[0,1]$ we have
$f(0)=G0\geqslant 0$ and $f(1)=Gx-1\leqslant 0$, so $f$ is a continuous
function on $[0,1]$ satisfying $f(0)\geqslant 0$ and $f(1)\leqslant 1$.
Existence of an $x$ satisfying $f(x)=0$ follows from the intermediate value
theorem. The same $x$ is a fixed point of $G$. One advantage of Brouwer’s
fixed point theorem is that its conditions are quite weak. One disadvantage is
that it provides only existence, without uniqueness or stability. Figure 6.4
provides an example of how multiple fixed points can coexist under the
conditions of the theorem.222There are many useful extensions to Brouwer’s
theorem, such as one for set-valued functions due to Kakutani, (1941). These
results have many applications in economics (see, for example, Nash, (1950)).
Another paper, due to Schauder, (1930), extends Brouwer’s result to infinite
dimensional spaces. Prove the first part of Theorem 4.1.1 (existence of a
stationary distribution) using Brouwer’s fixed point theorem. Let $P$ be a
Markov matrix on finite set $S$. Suppose in particular that $S$ has $d$
elements, so we can identify functions in $\mathbbm{R}^{S}$ with vectors in
$\mathbbm{R}^{d}$ and $\mathscr{D}(S)$ with the unit simplex in
$\mathbbm{R}^{d}$. The set $\mathscr{D}(S)$ is a closed, bounded and convex
subset of $\mathbbm{R}^{d}$ that $P$ maps into itself. As a linear matrix
operation, the map $\psi\mapsto\psi P$ is continuous. Existence of a fixed
point now follows from Brouwer’s fixed point theorem (p. 4.1.2).

######  4.1.3.4 Computation

Let’s consider how to compute stationary distributions. While
$\psi^{*}P=\psi^{*}$ is a finite set of linear equations that we might hope to
solve directly for $\psi^{*}$, there are problems with this idea. For example,
$\psi^{*}=0$ is a solution that fails to lie in $\mathscr{D}(S)$. To restrict
the solution to $\mathscr{D}(S)$ we proceed as follows: Suppose $|S|=n$ and
note that row vector $\psi\in\mathscr{D}(S)$ is stationary if and only if
$\psi(I-P)=0$, where $I$ is the $n\times n$ identity matrix. Let
$\mathbbm{1}_{n}$ be the $1\times n$ row vector $(1,\ldots,1)$. Let
$\mathbbm{1}_{n\times n}$ be the $n\times n$ matrix of ones. Consider the
linear system

$\mathbbm{1}_{n}=\psi(I-P+\mathbbm{1}_{n\times n})$ (4.15)

where $\psi$ is a row vector. Show that

  1. (i)

every solution $\psi$ of (4.15) lies in $\mathscr{D}(S)$ and

  2. (ii)

$\psi$ is stationary for $P$ if and only if (4.15) holds.

Solving the linear system (4.15) produces a stationary distribution in any
setting where the stationary distribution is unique. When this is not the
case, problems can arise. Give a counterexample to the claim that
$(I-P+\mathbbm{1}_{n\times n})$ is always nonsingular when $P$ is a stochastic
matrix. There are also graph-theoretic algorithms for computing all of the
stationary distributions of an arbitrary stochastic matrix. (The Python and
Julia libraries QuantEcon.py and QuantEcon.jl have efficient implementations
of this type.)

####  4.2 Asymptotics

In this section we turn to long run properties of Markov chains, including
ergodicity. With these ideas in hand, we will then investigate the evolution
of information over social networks.

#####  4.2.1 Ergodicity

Let’s begin with the fascinating and important topic of ergodicity. The
simplest way to understand ergodicity is as an extension of the law of large
numbers from iid sequences to more general settings, where the iid property
holds only in a limiting sense.

######  4.2.1.1 Harmonic Functions

Let $\mathscr{M}$ be a be a finite Markov model with state space $S$ and
transition matrix $P$. Fix $h\in\mathbbm{R}^{S}$ and let

$Ph(x):=\sum_{y\in S}P(x,y)h(y)\qquad(x\in S).$ (4.16)

If $h$ is understood as a column vector in $\mathbbm{R}^{|S|}$, then $Ph(x)$
is just element $x$ of the vector $Ph$. The map $h\mapsto Ph$ is sometimes
called the “conditional expectation operator,” since, given that
$P(x,y)=\mathbbm{P}\\{X_{t+1}=y\,|\,X_{t}=x\\}$, we have

$\sum_{y\in S}P(x,y)h(y)=\mathbbm{E}[h(X_{t+1})\,|\,X_{t}=x]$

A function $h\in\mathbbm{R}^{S}$ is called $P$-harmonic if $Ph=h$ pointwise on
$S$. Thus, $P$-harmonic functions are fixed points of the conditional
expectation operator $h\mapsto Ph$. If $h$ is $P$-harmonic and $(X_{t})$ is
$P$-Markov, then

$\mathbbm{E}\,[h(X_{t+1})\,|\,X_{t}]=(Ph)(X_{t})=h(X_{t}).$ (4.17)

(A stochastic process with this property—i.e., that the current value is the
best predictor of next period’s value—is called a _martingale_. Martingales
are one of the foundational concepts in probability, statistics and finance.)

######  Example 4.2.1.

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. If $A\subset S$ and $A^{c}:=S\setminus A$ are both absorbing for
$\mathscr{M}$, then $\mathbbm{1}_{A}$ and $\mathbbm{1}_{A^{c}}$ are both
$P$-harmonic. To see this, we apply Exercise 4.1.1.1 to obtain

$(P\mathbbm{1}_{A})(x)=\sum_{y\in S}P(x,y)\mathbbm{1}_{A}(y)=\sum_{y\in
A}P(x,y)=\begin{cases}1&\text{ if }x\in A\\\ 0&\text{ if }x\in
A^{c}.\end{cases}$

In other words, $P\mathbbm{1}_{A}=\mathbbm{1}_{A}$. A similar argument gives
$P\mathbbm{1}_{A^{c}}=\mathbbm{1}_{A^{c}}$.

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Show that every constant function in $\mathbbm{R}^{S}$ is
$P$-harmonic.

######  4.2.1.2 Definition and Implications

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. We know that every constant function in $\mathbbm{R}^{S}$ is
$P$-harmonic. We call $\mathscr{M}$ ergodic if the only $P$-harmonic functions
in $\mathbbm{R}^{S}$ are the constant functions.

######  Example 4.2.2.

If $P(x,y)=\varphi(y)$ for all $x,y$ in $S$, where $\varphi$ is some fixed
distribution on $S$, then $P$ generates an iid Markov chain, since the next
period draw has no dependence on the current state. Any model $\mathscr{M}$
with an adjacency matrix of this type is ergodic. Indeed, if $P$ has this
property and $h\in\mathbbm{R}^{S}$ is $P$-harmonic, then, for any $x\in S$, we
have

$h(x)=(Ph)(x)=\sum_{y}P(x,y)h(y)=\sum_{y}\varphi(y)h(y)$

Hence, $h$ is constant. This proves that any $P$-harmonic function is a
constant function.

######  Example 4.2.3.

A laborer is either unemployed (state $1$) or employed (state $2$). In state
$1$ he is hired with probability $\alpha$. In state $2$ he is fired with
probability $\beta$. The corresponding Markov model $\mathscr{M}$ has state
and transition matrix given by

$S=\\{1,2\\}\quad\text{and}\quad P_{w}=\begin{pmatrix}1-\alpha&\alpha\\\
\beta&1-\beta\end{pmatrix}.$

The statement $P_{w}h=h$ becomes

$\begin{pmatrix}1-\alpha&\alpha\\\
\beta&1-\beta\end{pmatrix}\begin{pmatrix}h(1)\\\
h(2)\end{pmatrix}=\begin{pmatrix}h(1)\\\ h(2)\end{pmatrix}$

The first row is $(1-\alpha)h(1)+\alpha h(2)=h(1)$, or $\alpha h(1)=\alpha
h(2)$. Thus, $\mathscr{M}$ is ergodic whenever $\alpha>0$. By a similar
argument, $\mathscr{M}$ is ergodic whenever $\beta>0$.

######  Example 4.2.4.

Let $\mathscr{M}$ be any finite Markov model with state space $S$ and
transition matrix $P$. It is immediate from Example 4.2.1 that if $S$ can be
partitioned into two nonempty absorbing sets, then $\mathscr{M}$ is not
ergodic. Hence, the poverty trap model in Figure 4.9 is not ergodic.
Similarly, if $\alpha=\beta=0$ in the matrix $P_{w}$ discussed in Example
4.2.3, then $P_{w}=I$ and each state is a disjoint absorbing set. Ergodicity
fails.

middle classrichpoor Figure 4.9:  A poverty trap

The examples discussed above suggest that, for a finite Markov model
$\mathscr{M}$, ergodicity depends on the connectivity properties of the
digraph. The next result confirms this.

######  Proposition 4.2.1.

If $\mathscr{M}$ is strongly connected, then $\mathscr{M}$ is ergodic.

###### Proof.

Let $\mathscr{M}$ be a strongly connected Markov model with state space $S$
and transition matrix $P$. Let $h$ be $P$-harmonic and let $x\in S$ be the
maximizer of $h$ on $S$. Let $m=h(x)$. Suppose there exists a $y$ in $S$ with
$h(y)<m$. Since $\mathscr{M}$ is strongly connected, we can choose a
$k\in\mathbbm{N}$ such that $P^{k}(x,y)>0$. Since $h$ is $P$-harmonic, we have
$h=P^{k}h$, and hence

$m=h(x)=\sum_{z}P^{k}(x,z)h(z)=P^{k}(x,y)h(y)+\sum_{z\not=y}P^{k}(x,z)h(z)<m.$

This contradiction shows that $h$ is constant. Hence $\mathscr{M}$ is ergodic.
∎

The implication in Proposition 4.2.1 is one way, as the next exercise asks you
to confirm. Provide an example of a finite Markov model that is ergodic but
not strongly connected. Let $\mathscr{M}$ be a finite Markov model with states
$S=\\{1,2\\}$ and edge set $E=\\{(1,2)\\}$. Thus, the chain immediately moves
to state 2 and stays there forever. The corresponding transition matrix is

$P=\begin{pmatrix}0&1\\\ 0&1\end{pmatrix}$

If $h$ is $P$-harmonic, then $h(x)=Ph(x)=h(2)$. Hence $h$ is constant. This
shows that $P$ is ergodic. At the same time $(S,E)$ is not strongly connected.

######  4.2.1.3 Implications of Ergodicity

One of the most important results in probability theory is the law of large
numbers (LLN). In the finite state setting, the classical version of this
theorem states that

$\mathbbm{P}\left\\{\lim_{k\to\infty}\frac{1}{k}\sum_{t=0}^{k-1}h(X_{t})=\sum_{x\in
S}h(x)\varphi(x)\right\\}=1$

when $(X_{t})_{t\geqslant 0}$ is an iid sequence of random variables with
common distribution $\varphi\in\mathscr{D}(S)$ and $h$ is an arbitrary element
of $\mathbbm{R}^{S}$. This version of the LLN is classical in the sense that
the iid assumption is imposed. It turns out that the iid assumption can be
weakened to allow for a degree of dependence between observations, which leads
us to ask whether or not the LLN holds for Markov chains as well. The answer
to this question is yes, provided that dependence between observations dies
out fast enough. One candidate for this condition is ergodicity. In fact, it
turns out that ergodicity is the exact necessary and sufficient condition
required to extend the law of large numbers to Markov chains. The next theorem
gives details.

######  Theorem 4.2.2.

Let $\mathscr{M}$ be any finite Markov model with state space $S$ and
transition matrix $P$. The following statements are equivalent:

  1. (i)

$\mathscr{M}$ is ergodic.

  2. (ii)

$\mathscr{M}$ has a unique stationary distribution $\psi^{*}$ and, for any
$P$-Markov chain $(X_{t})$ and any $h\in\mathbbm{R}^{S}$,

$\mathbbm{P}\left\\{\lim_{k\to\infty}\frac{1}{k}\sum_{t=0}^{k-1}h(X_{t})=\sum_{x\in
S}h(x)\psi^{*}(x)\right\\}=1.$ (4.18)

The proof of Theorem 4.2.2 is given in Chapter 17 of Meyn and Tweedie, (2009).
We skip the proof but provide intuition through the following examples.

######  Example 4.2.5.

Recall the iid case from Example 4.2.2. We showed that $\mathscr{M}$ is
ergodic. By Theorem 4.2.2, the convergence in (4.18) holds with
$\psi^{*}=\varphi$. This is consistent with the LLN for iid sequences.

######  Example 4.2.6.

Let $\mathscr{M}$ be a finite Markov model with $S=\\{1,2\\}$ and $P=I$, the
identity matrix. Markov chains generated by $P$ are constant. Since every
$h\in\mathbbm{R}^{2}$ satisfies $Ph=h$, we see that $\mathscr{M}$ is not
ergodic. This means that the LLN result in (4.18) fails. But this is exactly
what we would expect, since a constant chain $(X_{t})$ implies
$\frac{1}{k}\sum_{t=0}^{k-1}X_{t}=X_{0}$ for all $k$. In particular, if
$X_{0}$ is drawn from a nondegenerate distribution, then the sample mean does
not converge to any constant value.

######  Example 4.2.7.

Consider again the poverty trap model in Figure 4.9. Say that $h(x)$ is
earnings in state $x$, and that $h(\text{poor})=1$, $h(\text{middle})=2$ and
$h(\text{rich})=3$. Households that start with $X_{0}=\text{poor}$ will always
be poor, so $\frac{1}{k}\sum_{t=0}^{k-1}h(X_{t})=1$ for all $k$. Households
that start with $X_{0}$ in $\\{\text{middle, rich}\\}$ remain in this
absorbing set forever, so $\frac{1}{k}\sum_{t=0}^{k-1}h(X_{t})\geqslant 2$ for
all $k$. In particular, the limit of the sum depends on the initial condition.
This violates part (ii) of Theorem 4.2.2, which states that the limit is
independent of the distribution of $X_{0}$.

######  4.2.1.4 Reinterpreting the Stationary Distribution

Ergodicity has many useful implications. One is a new _interpretation_ for the
stationary distribution. To illustrate this, let $\mathscr{M}$ be an ergodic
Markov model with state space $S$ and transition matrix $P$. Let $(X_{t})$ be
a $P$-chain and let

$\hat{\psi}_{k}(y):=\frac{1}{k}\sum_{t=1}^{k}\mathbbm{1}\\{X_{t}=y\\}\qquad(y\in
S).$

The value $\hat{\psi}_{k}(y)$ measures the fraction of time that the $P$-chain
spends in state $y$ over the time interval $1,\ldots,k$. Under ergodicity, for
fixed $y\in S$, we can use (4.18) with $h(x)=\mathbbm{1}\\{x=y\\}$ to obtain

$\hat{\psi}_{k}(y)\to\sum_{x\in
S}\mathbbm{1}\\{x=y\\}\psi^{*}(x)=\psi^{*}(y).$ (4.19)

Turning (4.19) around, we see that

$\psi^{*}(y)\approx\text{ the fraction of time that any $P$-chain $(X_{t})$
spends in state $y$}.$ (4.20)

Figure 4.10 illustrates this idea for a simulated Markov chain $(X_{t})$
generated from the matrix $P_{B}$ introduced in (4.5). The figure compares
$\hat{\psi}_{k}$ and $\psi^{*}$ for different values of $k$. As $k\to\infty$,
the convergence claimed in (4.19) occurs.

Figure 4.10:  Convergence of the empirical distribution to $\psi^{*}$

Of course, we must remember that (4.20) is only valid under ergodicity. For
example, if $P=I$, the identity, then every distribution is stationary, every
$P$-Markov chain is constant, and (4.20) does not hold. Notice that, in view
of Theorem 4.2.2, the convergence in (4.18) occurs for _any_ initial
condition. This gives us, under ergodicity, a way of computing the stationary
distribution via simulation (and (4.19)). While this method is slower than
algebraic techniques (see, e.g., §4.1.3.4) for small problems, it can be the
only feasible option when $|S|$ is large.

#####  4.2.2 Aperiodicity and Stability

In §4.2.1 we discussed sample path properties of Markov chains, finding that
strong connectivity is sufficient for stability of time series averages. Next
we turn to marginal distributions of the chain and their long run properties,
which we examined informally in §4.1.2.2. We will see that, to guarantee
convergence of these sequences, we also need a condition that governs
periodicity.

######  4.2.2.1 Convergence of Marginals

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. We call $\mathscr{M}$ globally stable if there is only one
stationary distribution $\psi^{*}$ in $\mathscr{D}(S)$ and, moreover,

$\lim_{t\to\infty}\psi P^{t}=\psi^{*}\quad\text{for all
}\psi\in\mathscr{D}(S).$ (4.21)

In other words, marginal distributions of every $P$-Markov chain converge to
the unique stationary distribution of the model. There is a useful connection
between global stability and absorbing sets. Intuitively, if $\mathscr{M}$ is
globally stable and has an absorbing set $A$ that can be reached from any
$x\in S$, then any Markov chain generated by these dynamics will eventually
arrive in $A$ and never leave. Hence the stationary distribution must put all
of its mass on $A$. The next exercise asks you to confirm this. Let
$\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Let $A$ be a nonempty absorbing subset of $S$. Assume that

  1. (a)

$\mathscr{M}$ is globally stable with unique stationary distribution
$\psi^{*}$ and

  2. (b)

for each $x\in A^{c}:=S\setminus A$, there exists an $a\in A$ such that $a$ is
accessible from $x$,

Show that, under the stated conditions, there exists an $\varepsilon>0$ such
that $\sum_{y\in A^{c}}P^{n}(x,y)\\}\leqslant 1-\varepsilon$ for all $x\in S$,
where $n=|S|$. Using this fact, show in addition that $\sum_{y\in
A}\psi^{*}(y)=1$. Let $\mathscr{M}$ and $A$ have the stated properties. Fix
$x\in S$. Consider a $P$-chain $(X_{t})$ that starts at $x$. By property (b),
there exists an $a\in A$ and a $k\leqslant n:=|S|$ such that
$\varepsilon_{x}:=\mathbbm{P}\\{X_{k}=a\\}$ is strictly positive. Once
$(X_{t})$ enters $A$ it never leaves, so $X_{k}=a$ implies $X_{n}\in A$. Hence
$\mathbbm{P}\\{X_{n}\in A\\}\geqslant\varepsilon_{x}$. With
$\varepsilon:=\min_{x\in S}\varepsilon_{x}>0$, we then have
$\mathbbm{P}\\{X_{n}\in A\\}\geqslant\varepsilon$ for any initial condition
$x$. Another way to state this is $\sum_{y\in A^{c}}P^{n}(x,y)\\}\leqslant
1-\varepsilon$ for all $x\in S$, which proves the first claim. Regarding the
second claim, fix $\psi\in\mathscr{D}(S)$ and let $\psi_{t}=\psi P^{t}$ for
all $t$. Observe that, for fixed $m\in\mathbbm{N}$, we have
$\psi_{(m+1)n}=P^{n}\psi_{mn}$. As a result, for $y\in A^{c}$, we have

$\psi_{(m+1)n}(y)=\sum_{x\in S}P^{n}(x,y)\psi_{mn}(x)=\sum_{x\in
A^{c}}P^{n}(x,y)\psi_{mn}(x).$

Summing over $y$ gives

$\sum_{y\in A^{c}}\psi_{(m+1)n}(y)=\sum_{y\in A^{c}}\sum_{x\in
A^{c}}P^{n}(x,y)\psi_{mn}(x)=\sum_{x\in A^{c}}\left[\sum_{y\in
A^{c}}P^{n}(x,y)\right]\psi_{mn}(x).$

Let $\eta_{t}:=\sum_{x\in A^{c}}\psi_{t}(x)$ be the amount of probability mass
on $A^{c}$ at time $t$. Using the first claim and the definition of $\eta_{t}$
now gives $\eta_{(m+1)n}\leqslant(1-\varepsilon)\eta_{mn}$. Hence
$\eta_{mn}\to 0$ as $m\to\infty$. At the same time, $\psi_{mn}\to\psi^{*}$ as
$m\to\infty$, so

$\sum_{x\in A^{c}}\psi^{*}(x)=\lim_{m\to\infty}\sum_{x\in
A^{c}}\psi_{mn}(x)=\lim_{m\to\infty}\eta_{mn}=0.$

The second claim is now verified.

######  4.2.2.2 A Key Theorem

From Theorem 4.1.1 we saw that strong connectedness is sufficient for
uniqueness of the stationary distribution. One might hope that strong
connectedness is enough for global stability too, but this is not true. To see
why, suppose for example that $S=\\{0,1\\}$ and $E=\\{(0,1),(1,0)\\}$. The
transition matrix is

$P_{d}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}.$ (4.22)

While this model is strongly connected, global stability fails. Indeed, if
$(X_{t})$ is $P_{d}$-Markov and starts at $0$, then $(X_{t})$ will visit state
1 on odd dates and state 0 on even dates. That is,
$P^{t}_{d}\,\delta_{0}=\delta_{t\text{ mod }2}$. This sequence does not
converge. The issue with $P_{d}$ is that, even though the chain traverses the
whole state space, the distribution of $X_{t}$ will affect that of $X_{t+j}$
for all $j$ due to periodicity. This causes stability to fail. If, however, we
rule out periodicity, then we have enough for stability to hold. This line of
reasoning leads to the following famous theorem.

######  Theorem 4.2.3.

Let $\mathscr{M}$ be a finite Markov model. If $\mathscr{M}$ is strongly
connected and aperiodic, then $\mathscr{M}$ is globally stable.

###### Proof.

Let $\mathscr{M}$ be a finite Markov model with state $S$ and transition
matrix $P$. Suppose $\mathscr{M}$ is aperiodic and strongly connected. Let
$\psi^{*}$ be the unique stationary distribution of $\mathscr{M}$. By Theorem
1.4.4, $P$ is primitive. Hence we can apply the last part of the
Perron–Frobenius Theorem (see page 1.2.6). Using $r(P)=1$, this result tells
us that $P^{t}\to e_{r}e_{\ell}^{\top}$ as $t\to\infty$, where $e_{r}$ and
$e_{\ell}$ are the unique right and left eigenvectors of $P$ corresponding to
the eigenvalue $r(P)=1$ that also satisfy the normalization
$\left\langle{e_{\ell},e_{r}}\right\rangle=1$. Now observe that $\psi^{*}$
obeys $\psi^{*}=\psi^{*}P$ and, in addition, $P\mathbbm{1}=\mathbbm{1}$. Hence
$\psi^{*}$ and $\mathbbm{1}$ are left and right eigenvectors corresponding to
$r(P)=1$. Moreover, $\left\langle{\psi^{*},\mathbbm{1}}\right\rangle=1$. Hence
$e_{r}=\mathbbm{1}$ and $e_{\ell}^{\top}=\psi^{*}$. Combining these facts
leads to

$\lim_{t\to\infty}P^{t}=P^{*}\quad\text{where}\quad
P^{*}:=\mathbbm{1}\psi^{*}.$ (4.23)

If we pick any $\psi\in\mathscr{D}(S)$, then, by (4.23), we get $\psi
P^{t}\to\left\langle{\psi,\mathbbm{1}}\right\rangle\psi^{*}=\psi^{*}$. Hence
global stability holds. ∎

######  Example 4.2.8.

The Markov models $P_{Q}$ and $P_{B}$ in §4.1.1 are both aperiodic by the
results in §1.4.1.4. Being strongly connected, they are also globally stable.

The aperiodicity condition in Theorem 4.2.3 is, in general, not stringent. On
the other hand, the strong connectedness requirement is arguably quite strict.
Weaker conditions for global stability are available, as we shall see in
§4.2.3.

######  4.2.2.3 Rates of Convergence: Spectral Gap

While global stability is a very useful property, the implications are
qualitative rather than quantitative. In practice, we usually want to know
something about the _rate_ of convergence. There are several ways of looking
at this. We cover the two most common in this section and the next. As a first
step, fixing a Markov model $\mathscr{M}$ with state space $S$ and transition
matrix $P$, we use (1.4) on page 1.4 to write $P^{t}$ as

$P^{t}=\sum_{i=1}^{n-1}\lambda_{i}^{t}e_{i}\varepsilon_{i}^{\top}+\mathbbm{1}\psi^{*},$
(4.24)

where each $\lambda_{i}$ is an eigenvalue of $P$ and $e_{i}$ and
$\varepsilon_{i}$ are the right and left eigenvectors corresponding to
$\lambda_{i}$. We have also ordered the eigenvalues from smallest to largest,
and used the Perron–Frobenius Theorem to infer that $\lambda_{n}=r(P)=1$, as
well as the arguments in the proof of Theorem 4.2.3 that showed
$e_{n}=\mathbbm{1}$ and $\varepsilon_{n}^{\top}=\psi^{*}$. Premultiplying
$P^{t}$ by arbitrary $\psi\in\mathscr{D}(S)$ and rearranging now gives

$\psi P^{t}-\psi^{*}=\sum_{i=1}^{n-1}\lambda_{i}^{t}\psi
e_{i}\varepsilon_{i}^{\top}$ (4.25)

Recall that eigenvalues are ordered from smallest to largest and. Moreover, by
the Perron–Frobenius Theorem, $\lambda_{i}<1$ for all $i<n$ when $P$ is
primitive (i.e., $\mathscr{M}$ is strongly connected and aperiodic). Hence,
after taking the Euclidean norm deviation, we obtain

$\|\psi
P^{t}-\psi^{*}\|=O(\eta^{t})\quad\text{where}\quad\eta:=|\lambda_{n-1}|<1.$
(4.26)

Thus, the rate of convergence is governed by the modulus of the second largest
eigenvalue. The difference between the largest and second largest eigenvalue
of a nonnegative matrix is often called the spectral gap. For this reason, we
can also say that, for primitive stochastic matrices, the rate of convergence
is determined by the (nonzero) spectral gap.

######  Example 4.2.9.

When studying the worker model with hiring rate $\alpha$ and firing rate
$\beta$ in §1.2.1.3, we found that the eigenvalues of the transition matrix
$P_{w}$ are $\lambda_{1}=1-\alpha-\beta$ and $\lambda_{2}=1$. Hence the
spectral gap is $\alpha+\beta$ and the rate of convergence is
$O((1-\alpha-\beta)^{t})$. High hiring and firing rates both produce faster
convergence. In essence, this is because higher hiring and firing rates mean
workers do not stay in any state for long, so initial conditions die out
faster.

#####  4.2.3 The Markov–Dobrushin Coefficient

The rate of convergence of $\psi P^{t}$ to $\psi^{*}$ given in §4.2.2.3
restricts the Euclidean distance between these vectors as a function of $t$.
There are, however, other metrics we could use in studying rates of
convergence, and sometimes these other metrics give more convenient results.
In fact, as we show in this section, a good choice of metric leads us to a
more general stability result than the (better known) aperiodicity-based
result in §4.2.2.

######  4.2.3.1 An Alternative Metric

For the purposes of this section, For $\varphi,\psi\in\mathscr{D}(S)$, we set

$\rho(\varphi,\psi):=\|\varphi-\psi\|_{1}:=\sum_{x\in S}|\varphi(x)-\psi(x)|,$

which is just the $\ell_{1}$ deviation between $\varphi$ and $\psi$ (see
§2.3.1.1). Show that, for any $\varphi,\psi\in\mathscr{D}(S)$, we have

  1. (i)

$\rho(\varphi,\psi)\leqslant 2$.

  2. (ii)

$\rho(\varphi P,\psi P)\leqslant\rho(\varphi,\psi)$ for any stochastic matrix
$P$

Regarding part (i), the triangle inequality, combined with the assumption that
$\varphi,\psi\in\mathscr{D}(S)$, gives the bound

$\sum_{x}|\varphi(x)-\psi(x)|\leqslant\sum_{x}|\varphi(x)|+\sum_{x}|\psi(x)|\leqslant
2.$

Regarding part (ii), if $P$ is a stochastic matrix, then

$\rho(\varphi P,\psi
P)=\sum_{y}\left|\sum_{x}P(x,y)\varphi(x)-\sum_{x}P(x,y)\psi(x)\right|\leqslant\sum_{y}\sum_{x}P(x,y)\left|\varphi(x)-\psi(x)\right|.$

Swapping the order of summation and using $\sum_{y}P(x,y)=1$ proves the claim.
Property (ii) is called the nonexpansiveness property of stochastic matrices
under the $\ell_{1}$ deviation. As we will see, we can tighten this bound when
$P$ satisfies certain properties. As a first step we note that, for the
$\ell_{1}$ deviation, given any stochastic matrix $P$, we have

$\rho(\varphi P,\psi P)\leqslant(1-\alpha(P))\rho(\varphi,\psi),$ (4.27)

where

$\alpha(P):=\min\left\\{\sum_{y\in S}[P(x,y)\wedge
P(x^{\prime},y)]\;:\;(x,x^{\prime})\in S\times S\right\\}.$ (4.28)

Here $a\wedge b:=\min\\{a,b\\}$ when $a,b\in\mathbbm{R}$. We call $\alpha(P)$
the Markov–Dobrushin coefficient of $P$, although other names are also used in
the literature. A proof of the bound in (4.27) can be found in Stachurski,
2022a  or Seneta, 2006a .333Seneta, 2006a  also discusses the history of
Andrey Markov’s work, which originated the kinds of contraction based
arguments used in this section. Consider the stochastic matrices

$P_{w}=\begin{pmatrix}1-\alpha&\alpha\\\ \beta&1-\beta\end{pmatrix}.$

Show that $\alpha(P_{w})=0$ if and only if $\alpha=\beta=0$ or
$\alpha=\beta=1$. How should the Markov–Dobrushin coefficient be understood?
The coefficient is large when the rows of $P$ are relatively similar. For
example, if rows $P(x,\cdot)$ and $P(x^{\prime},\cdot)$ are identical, the
$\sum_{y\in S}[P(x,y)\wedge P(x^{\prime},y)]=1$. Similarity of rows is related
to stability. The next exercise helps to illustrate. Let $P$ be such that all
rows are identical and equal to $\varphi\in\mathscr{D}(S)$. Prove that global
stability holds in one step, in the sense that $\varphi$ is the unique
stationary distribution and $\psi P=\varphi$ for all $\psi\in\mathscr{D}(S)$.
Using (4.27), show that, for any $\varphi,\psi\in\mathscr{D}(S)$, we have

$\rho(\varphi P^{t},\psi
P^{t})\leqslant(1-\alpha(P))^{t}\rho(\varphi,\psi)\quad\text{for all
}\;t\in\mathbbm{N}.$ (4.29)

Fix $\varphi,\psi\in\mathscr{D}(S)$. From (4.27) we know that (4.29) is true
when $t=1$. Now suppose it is true at $t$. Then, using the fact that (4.27)
holds for any pair of distributions,

$\rho(\varphi P^{t+1},\psi P^{t+1})\leqslant(1-\alpha(P))\rho(\varphi
P^{t},\psi P^{t})\leqslant(1-\alpha(P))^{t+1}\rho(\varphi,\psi)$

where the last step uses the induction hypothesis. Hence (4.29) also holds at
$t+1$, and, by induction, at all $t\in\mathbbm{N}$. Since powers of stochastic
matrices are again stochastic, and since (4.27) is valid for any stochastic
matrix, the bound in (4.29) can be generalized by replacing $P$ with $P^{k}$
for any given $k\in\mathbbm{N}$, which gives (with $t$ replaced by $\tau$)

$\rho(\varphi P^{\tau k},\psi P^{\tau
k})\leqslant(1-\alpha(P^{k}))^{\tau}\rho(\varphi,\psi)\quad\text{for all
}\;\tau\in\mathbbm{N}.$

Now observe that, for any $t\in\mathbbm{N}$, we can write $t=\tau k+j$, where
$\tau\in\mathbbm{Z}_{+}$ and $j\in\\{0,\ldots,k-1\\}$. Fixing $t$ and choosing
$j$ to make this equality hold, we get

$\rho(\varphi P^{t},\psi P^{t})=\rho(\varphi P^{\tau k+j},\psi P^{\tau
k+j})\leqslant\rho(\varphi P^{\tau k},\psi P^{\tau
k})\leqslant(1-\alpha(P^{k}))^{\tau}\rho(\varphi,\psi)$

where the second inequality is due to the nonexpansive property of stochastic
matrices (Exercise 4.2.3.1). Since $\tau$ is an integer satisfying
$\tau=\lfloor t/k\rfloor$, where $\lfloor\cdot\rfloor$ is the floor function,
we can now state the following.

######  Theorem 4.2.4.

Let $\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. For all $\varphi,\psi\in\mathscr{D}(S)$ and all
$k,t\in\mathbbm{N}$, we have

$\rho(\varphi P^{t},\psi P^{t})\leqslant(1-\alpha(P^{k}))^{\lfloor
t/k\rfloor}\rho(\varphi,\psi).$ (4.30)

In particular, if there exists a $k\in\mathbbm{N}$ such that
$\alpha(P^{k})>0$, then $\mathscr{M}$ is globally stable.

To see why the global stability implication stated in Theorem 4.2.4 holds,
suppose there exists a $k\in\mathbbm{N}$ such that $\alpha(P^{k})>0$. Now
substitute arbitrary $\psi\in\mathscr{D}(S)$ and any stationary distribution
$\psi^{*}$ for $P$ into (4.30) to obtain

$\rho(\psi P^{t},\psi^{*})\leqslant(1-\alpha(P^{k}))^{\lfloor
t/k\rfloor}\rho(\psi,\psi^{*})\leqslant 2(1-\alpha(P^{k}))^{\lfloor
t/k\rfloor}$ (4.31)

for all $t\in\mathbbm{N}$, where the second bound is due to Exercise 4.2.3.1.
In the preceding discussion, the distribution $\psi^{*}$ was taken to be any
stationary distribution of $P$. Using (4.31), prove that $P$ has only one
stationary distribution when $\alpha(P^{k})>0$. One of the major advantages of
Theorem 4.2.4 is that strong connectivity of $\mathscr{M}$ is not required. In
the next section we will see an example of a finite Markov model $\mathscr{M}$
where (i) strong connectivity fails but (ii) the conditions of Theorem 4.2.4
are satisfied.444It can be shown that the condition $\alpha(P^{k})>0$ for some
$k\in\mathbbm{N}$ is necessary as well as sufficient for global stability of a
finite Markov model. Hence the conditions of Theorem 4.2.4 are strictly weaker
than strong connectedness plus aperiodicity, as used in Theorem 4.1.1. Let
$\mathscr{M}$ be a finite Markov model with state space $S$ and transition
matrix $P$. Prove the existence of a $k\in\mathbbm{N}$ with $\alpha(P^{k})>0$
whenever $\mathscr{M}$ is strongly connected and aperiodic. If $\mathscr{M}$
is strongly connected and aperiodic, then $P$ is primitive, in which case
there exists a $k\in\mathbbm{N}$ with $P^{k}\gg 0$. Clearly $\alpha(P^{k})>0$.

######  4.2.3.2 Sufficient Conditions

While the Markov–Dobrushin coefficient $\alpha(P^{k})$ can be calculated for
any given $k$ on a computer by stepping through all pairs of rows in $P^{k}$,
this calculation is computationally intensive when $S$ is large. When
applicable, the following lemma simplifies the problem by providing a
sufficient condition for $\alpha(P^{k})>0$.

######  Lemma 4.2.5.

Let $k$ be a natural number and let $\mathscr{M}$ be a finite Markov model
with state space $S$ and transition matrix $P$. If there is a state $z\in S$
such that, for every $x\in S$, there exists a directed walk from $x$ to $z$ of
length $k$, then $\alpha(P^{k})>0$.

###### Proof.

Let $k\in\mathbbm{N}$ and $z\in S$ be such that, for every $x\in S$, there
exists a directed walk from $x$ to $z$ of length $k$. By Proposition 1.4.2, we
then have $r:=\min_{x\in S}P^{k}(x,z)>0$. Since, for any $x,x^{\prime}\in S$,

$\sum_{y\in S}[P^{k}(x,y)\wedge P^{k}(x^{\prime},y)]\geqslant P^{k}(x,z)\wedge
P^{k}(x^{\prime},z)\geqslant r>0,$

strict positivity of $\alpha(P^{k})$ now follows from the definition of the
Markov–Dobrushin coefficient. ∎

######  Example 4.2.10.

Consider the digraph in Figure 4.11. This digraph is not strongly connected
because 4 is not accessible from anywhere. However, there exists a directed
walk from any vertex to vertex 1 in $k=2$ steps. For example, from 2 we can
choose $2\to 1$ and then $1\to 1$, from 1 we can choose $1\to 1$ and then
$1\to 1$, etc. Hence, if Figure 4.11 is the digraph of a finite Markov model
with transition matrix $P$, then $\alpha(P^{2})>0$.

1234 Figure 4.11:  A digraph with walk of length 2 from any node to 1

######  Example 4.2.11.

Consider the Markov dynamics suggested in Figure 1.10 on page 1.10. Although
there are no weights, we can see that poor is accessible from every state in
one step, so $\mathscr{M}$ must be globally stable. In addition, poor is, by
itself, an absorbing set. Hence, by Exericse 4.2.2.1, for any choice of
weights compatible with these edges, the stationary distribution will
concentrate all its mass on poor.

######  Remark 4.2.1.

As was pointed out in proof of Lemma 4.2.5, under the conditions of that lemma
we have $P^{k}(x,z)>0$ for all $x\in S$. This means that $\alpha(P^{k})>0$
whenever $P^{k}$ has a strictly positive column.

######  4.2.3.3 Application: PageRank

In §1.4.3 we discussed centrality measures for networks. Centrality measures
provide a ranking of vertices in the network according to their “centrality”
or “importance.” One of the most important applications of ranking of vertices
in a network is ranking the importance of web pages on the internet.
Historically, the most prominent example of a ranking mechanism for the
internet is PageRank, which transformed Google from a minor start up to a
technology behemoth. In this section we provide a simple introduction to the
original form of PageRank and connect it to previously discussed measures of
centrality. Consider a finite collection of web pages $W$ and let $L$ be the
hyperlinks between them. We understand $(W,L)$ as a digraph $\mathscr{G}$,
where $W$ is the vertices and $L$ is the edges. Let $A$ be the associated
adjacency matrix, so that $A(u,v)=1$ if there is a link from $u$ to $v$ and
zero otherwise. We set $n=|W|$, so that $A$ is $n\times n$. To start our
analysis, we consider the case where $\mathscr{G}$ is strongly connected, such
as the small network in Figure 4.12. Furthermore, we introduce a second matrix
$P$ in which each row of $A$ has been normalized so that it sums to one. For
the network in Figure 4.12, this means that

$A=\begin{pmatrix}0&1&1&1\\\ 0&0&1&0\\\ 0&0&0&1\\\
1&0&0&0\end{pmatrix}\quad\text{and}\quad P=\begin{pmatrix}0&1/3&1/3&1/3\\\
0&0&1&0\\\ 0&0&0&1\\\ 1&0&0&0\end{pmatrix}.$

Now consider an internet surfer who, once per minute, randomly clicks on one
of the $k$ outgoing links on a page, each link selected with uniform
probability $1/k$. The idea of PageRank is to assign to each page $u\in W$ a
value $g(u)$ equal to the fraction of time that this surfer spends on page $u$
over the long run. Intuitively, a high value for $g(u)$ indicates a heavily
visited and hence important site.

1234 Figure 4.12:  A digraph with walk of length 2 from any node to 1

The vector $g$ is easy to compute, given our knowledge of Markov chains. Let
$\mathscr{M}$ be the finite Markov model associated with the random surfer,
with state space $W$ and adjacency matrix given by $P$. Since $\mathscr{M}$ is
strongly connected, the ergodicity theorem (page 4.2.2) tells us that $P$ has
a unique stationary distribution $\psi^{*}$, and that the fraction of time the
surfer spends at page $u$ is equal to the probability assigned to $u$ under
the stationary distribution (see, in particular, (4.20)). Hence $g=\psi^{*}$.
As $g$ is stationary and $r(P)=1$, we can write $g=(1/r(P))gP$. Taking
transposes gives $g^{\top}=(1/r(P))P^{\top}g^{\top}$. Comparing with (1.32) on
page 1.32, we see that, for this simple case, the PageRank vector $g$ is just
the authority-based eigenvector centrality measure of $\mathscr{M}$. Thus,
PageRank gives high ranking to pages with many inbound links, attaching high
value to inbound links from pages that are themselves highly ranked. There are
two problems with the preceding analysis. First, we assumed that the internet
is strongly connected, which is clearly violated in practice (we need only one
page with no outbound links). Second, internet users sometimes select pages
without using hyperlinks, by manually entering the URL. The PageRank solution
to this problem is to replace $P$ with the so-called Google matrix

$G:=\delta P+(1-\delta)\frac{1}{n}\mathbbm{1},$

where $\mathbbm{1}$ is the $n\times n$ matrix of ones. The value
$\delta\in(0,1)$ is called the damping factor. Prove that $G$ is a stochastic
matrix for all $\delta\in(0,1)$. The Markov dynamics embedded in the
stochastic matrix $G$ can be understood as follows: The surfer begins by
flipping a coin with heads probability $\delta$. (See Exercise 1.3.1.2 and its
solution for the connection between convex combinations and coin flips.) If
the coin is heads then the surfer randomly selects and follows one of the
links on the current page. If not then the surfer randomly selects and moves
to any page on the internet. For given $\delta$, the PageRank vector for this
setting is adjusted to be the stationary distribution of the Google matrix
$G$. Verify that the digraph associated with the transition probabilities in
$G$ is always strongly connected (assuming, as above, that $\delta\in(0,1)$).
All elements of $G$ are strictly positive, so a directed edge exists between
every pair of pages $u,v\in W$. This clearly implies strong connectedness. As
a result of Exercise 4.2.3.3, we can always interpret the stationary of $G$ as
telling us the fraction of time that the surfer spends on each page in the
long run. Use (4.31) to obtain a rate of convergence of $\psi G^{t}$ to the
adjusted PageRank vector $g^{*}$ (i.e., the unique stationary distribution
$g^{*}$ of $G$), where $\psi$ is an arbitrary initial distribution on $W$.
(Set $k=1$.) For any $u,u^{\prime},v\in W$, we have $G(u,v)\wedge
G(u^{\prime},v)\geqslant 1-\delta$. Hence $\alpha(G)\geqslant 1-\delta$.
Therefore, by (4.31), we have

$\rho(\psi G^{t},g^{*})\leqslant 2\delta^{t}.$

#####  4.2.4 Information and Social Networks

In recent years, the way that opinions spread across social networks has
become a major topic of concern in many countries around the globe. A well-
known mathematical model of this phenomenon is De Groot learning, which was
originally proposed in DeGroot, (1974). This mechanism has linear properties
that make it relatively easy to analyze (although large and complex underlying
networks can cause significant challenges). In De Groot learning, a group of
agents, labeled from $1$ to $n$, is connected by a social or information
network of some type. Connections are indicated by a trust matrix $T$, where,
informally,

$T(i,j)=\text{amount that $i$ trusts the opinion of $j$}.$

In other words, $T(i,j)$ is large if agent $i$ puts a large positive weight on
the opinion of agent $j$. The matrix $T$ is assumed to be stochastic. We can
view the trust matrix as an adjacency matrix for a weighed digraph
$\mathscr{S}$ with vertex set $V:=[n]$ and edges

$E=\\{(i,j)\in V\times V\,:\,T(i,j)>0\\}.$

1234$0.5$$0.4$$0.4$$0.5$$0.5$$0.5$$0.1$$0.1$$1.0$ (a) Social network
$\mathscr{S}_{a}$

1234$0.2$$0.5$$0.4$$0.4$$0.5$$0.5$$0.5$$0.1$$0.1$$0.8$ (b) Social network
$\mathscr{S}_{b}$

Figure 4.13: Two social networks

Figure 4.13 shows two social networks $\mathscr{S}_{a}$ and $\mathscr{S}_{b}$
with trust matrices given by

$T_{a}=\begin{pmatrix}1&0&0&0\\\ 0.5&0.1&0.4&0\\\ 0.5&0.4&0.1&0\\\
0&0.5&0&0.5\end{pmatrix}\quad\text{and}\quad
T_{b}=\begin{pmatrix}0.8&0&0&0.2\\\ 0.5&0.1&0.4&0\\\ 0.5&0.4&0.1&0\\\
0&0.5&0&0.5\end{pmatrix}$

respectively. In network A, agent 1 places no trust in anyone’s opinion but
his own. In network B, he places at least some trust in the opinion of agent
4\. Below we show how these differences matter for the dynamics of beliefs.

######  4.2.4.1 Learning

At time zero, all agents have an initial subjective belief concerning the
validity of a given statement. Belief takes values in $[0,1]$, with 1
indicating complete (subjective) certainty that the statement is true. Let
$b_{0}(i)$ be the belief of agent $i$ at time zero. An agent updates beliefs
sequentially based on the beliefs of others, weighted by the amount of trust
placed in their opinion. Specifically, agent $i$ updates her belief after one
unit of time to $\sum_{j=1}^{n}T(i,j)b_{0}(j)$. More generally, at time $t+1$,
beliefs update to

$b_{t+1}(i)=\sum_{j=1}^{n}T(i,j)b_{t}(j)\qquad(i\in V).$ (4.32)

In matrix notation this is $b_{t+1}=Tb_{t}$, where each $b_{t}$ is treated as
a column vector. (Notice that this update rule is similar but not identical to
the marginal distribution updating rule for Markov chains (the forward
equation) discussed on page 4.1.2.1. Here we are postmultiplying by a column
vector rather than premultiplying be a row vector.) If some subgroup of agents
$U\subset V$ is an absorbing set for the digraph $\mathscr{S}$, indicating
that members of this group place no trust in outsiders, then the initial
beliefs $\\{b_{0}(i)\\}_{i\in U^{c}}$ of the outsiders (members of
$U^{c}=V\setminus U$) have no influence on the beliefs of the insiders
(members of $U$) at any point in time. Prove that this is true.

######  4.2.4.2 Consensus

A social network $\mathscr{S}$ is said to lead to consensus if
$|b_{t}(i)-b_{t}(j)|\to 0$ as $t\to\infty$ for all $i,j\in V$. Consensus
implies that all agents eventually share the same belief. An important
question is, what conditions on the network lead to a consensus outcome?

######  Proposition 4.2.6.

If there exists a $k\in\mathbbm{N}$ such that $\alpha(T^{k})>0$, then
consensus is obtained. In particular,

$|b_{t}(i)-b_{t}(j)|\leqslant 2(1-\alpha(T^{k}))^{\lfloor
t/k\rfloor}\quad\text{for all $t\in\mathbbm{N}$ and $i,j\in V$}.$ (4.33)

###### Proof.

Fix $i,j\in V$ and $t\in\mathbbm{N}$. Let $b=b_{0}$. For any
$\varphi,\psi\in\mathscr{D}(V)$, an application of the triangle inequality
gives

$|\varphi T^{t}b-\psi T^{t}b|=\left|\sum_{j}(\varphi
T^{t})(j)b(j)-\sum_{j}(\psi
T^{t})(j)b(j)\right|\leqslant\sum_{j}\left|(\varphi T^{t})(j)-(\psi
T^{t})(j)\right|,$

where we have used the fact that $|b(j)|=b(j)\leqslant 1$. Applying the
definition of the $\ell_{1}$ deviation and (4.30), we obtain the bound

$|\varphi T^{t}b-\psi T^{t}b|\leqslant(1-\alpha(T^{k}))^{\lfloor
t/k\rfloor}\rho(\varphi,\psi)\leqslant 2(1-\alpha(T^{k}))^{\lfloor
t/k\rfloor}.$ (4.34)

Since this bound is valid for any choice of $\varphi,\psi\in\mathscr{D}(V)$,
we can specialize to $\varphi=\delta_{i}$ and $\psi=\delta_{j}$ and, applying
$T^{t}b=b_{t}$, get

$|b_{t}(i)-b_{t}(j)|=|\delta_{i}T^{t}b-\delta_{j}T^{t}b|\leqslant
2(1-\alpha(T^{k}))^{\lfloor t/k\rfloor}.\qed$

Proposition 4.2.6 can be applied to the two social networks $\mathscr{S}_{a}$
and $\mathscr{S}_{b}$ in Figure 4.13. For example, in network
$\mathscr{S}_{a}$, for every node $i$, there exists a walk of length 2 from
$i$ to node $1$. Hence, by Lemma 4.2.5, we have $\alpha(T^{2})>0$. In network
$\mathscr{S}_{b}$ the same is true. Let $\mathscr{S}$ be a social network with
trust (and adjacency) matrix $T$. Use Proposition 4.2.6 to show that
$\mathscr{S}$ leads to consensus whenever $\mathscr{S}$ is strongly connected
and aperiodic. If $\mathscr{S}$ is strongly connected and aperiodic, then the
adjacency matrix $T$ is primitive, so there exists a $k\in\mathbbm{N}$ such
that $T^{k}\gg 0$. Hence $\alpha(T^{k})>0$, and Proposition 4.2.6 applies.

######  4.2.4.3 Influence of Authorities

Now let’s consider what beliefs converge to when consensus emerges. In
particular, we are interested in discovering whose opinions are most
influential under De Groot learning, for a given trust matrix. To answer the
question, let $\mathscr{S}$ be a given social network with trust matrix $T$.
Suppose that $\alpha(T^{k})>0$ for some $k\in\mathbbm{N}$. By Theorem 4.2.4,
the network $\mathscr{S}$ is globally stable when viewed as a finite Markov
model. Let $\psi^{*}$ be the unique stationary distribution, so that
$\psi^{*}=\psi^{*}T$. Applying (4.34) with $\varphi=\delta_{i}$ and
$\psi=\psi^{*}$ yields

$\left|b_{t}(i)-b^{*}\right|\leqslant(1-\alpha(T^{k}))^{\lfloor
t/k\rfloor}\rho(\varphi,\psi)\leqslant 2(1-\alpha(T^{k}))^{\lfloor
t/k\rfloor}$

where

$b^{*}:=\sum_{j\in V}\psi^{*}(j)b_{0}(j).$

We conclude that the belief of every agent converges geometrically to $b^{*}$,
which is a weighted average of the initial beliefs of all agents. In
particular, those agents with high weighting under the stationary distribution
have a large influence on these equilibrium beliefs. We can interpret this
through notions of centrality. Since $r(T)=1$, we have
$(\psi^{*})^{\top}=(1/r(T))T^{\top}(\psi^{*})^{\top}$, so $\psi^{*}$ is the
authority-based eigenvector centrality measure on $\mathscr{S}$. Thus, the
influence of each agent on long run beliefs is proportional to their
authority-based eigenvector centrality. This makes sense because such agents
are highly trusted by many agents who are themselves highly trusted. For
network $\mathscr{S}_{a}$ in Figure 4.13, show that $b^{*}=b_{0}(1)$. That is,
all agents’ beliefs converge to the belief of agent 1. Using a computer, show
that the stationary distribution for $\mathscr{S}_{b}$, rounded to two decimal
places, is $\psi^{*}=(0.56,0.15,0.07,0.22)$. (Notice how the relatively slight
change in network structure from $\mathscr{S}_{a}$ to $\mathscr{S}_{b}$
substatially reduces the influence of agent 1.)

####  4.3 Chapter Notes

High quality treatments of finite-state Markov dynamics include Norris,
(1998), Privault, (2013) and Häggström et al., (2002). For the general state
case see Meyn and Tweedie, (2009). A review of De Groot learning is available
in Jackson, (2010). Some interesting extensions related to the “wisdom of
crowds” phenomenon are provided in Golub and Jackson, (2010). Acemoglu et al.,
2021a  study misinformation and echo chambers in information networks. Board
and Meyer-ter Vehn, (2021) analyze learning dynamics in continuous time on
large social networks. Shiller, (2020) provides an interesting discussion of
how ideas spread across social networks and shape economic outcomes.

###  Chapter 5 Nonlinear Interactions

Much of what makes network analysis interesting is how ramifications of
choices flow across networks. In general, decisions made at a given node $i$
affect responses of neighboring nodes and, through them, neighbors of
neighboring nodes, and so on. As these consequences flow across the network,
they in turn affect choices at $i$. This is a tail-chasing scenario, which can
be unraveled through fixed point theory. In some network settings, such as the
input-output model in §2.1, interactions are linear and fixed point problems
reduce a system of linear equations. In other settings, interactions are
inherently nonlinear and, as a result, we need more sophisticated fixed point
theory. This chapter is dedicated to the study of networks with nonlinear
interactions. We begin with relevant fixed point theory and then apply it to a
sequence of problems that arise in analysis of economic networks, including
production models with supply constraints and financial networks.

####  5.1 Fixed Point Theory

Let $S$ be any set. Recall from §6.1.1.6 that, given a self-map $G$ on $S$, a
point $x\in S$ is called a fixed point of $G$ if $Gx=x$. (A self-map on $S$ is
a function $G$ from $S$ to itself. When working with self-maps it is common to
abbreviate $G(x)$ to $Gx$.) In this chapter, we will say that $G$ is globally
stable on $S$ if $G$ has a unique fixed point $x^{*}\in S$ and $G^{k}x\to
x^{*}$ as $k\to\infty$ for all $x\in S$. In other words, under this property,
the fixed point is not only unique but also globally attracting under
iteration of $G$. We have already discussed fixed points, indirectly or
directly, in multiple contexts:

  * •

In Chapter 2 we studied the equation $x=Ax+d$, where $x$ is an output vector,
$A$ is a matrix of coefficients and $d$ is a demand vector. A solution $x$ to
this equation can also be thought of as a fixed point of the affine map
$Fx=Ax+d$.

  * •

In Chapter 4 we learned that a stationary distribution of a finite Markov
model with state space $S$ and adjacency matrix $P$ is a
$\psi\in\mathscr{D}(S)$ with $\psi=\psi P$. In other words, $\psi$ is a fixed
point of $\psi\mapsto\psi P$ in $\mathscr{D}(S)$.

  * •

In Chapter 3 we studied the Bellman equation
$q(x)=\min_{y\in\mathscr{O}(x)}\\{c(x,y)+q(y)\\}$ and introduced an operator
the Bellman operator $T$, with the property that its fixed points exactly
coincide with solutions to the Bellman equation.

In each case, when we introduced these fixed point problems, we immediately
needed to consider questions of existence and uniqueness of fixed points. Now
we address these same issues more systematically in an abstract setting.

#####  5.1.1 Contraction Mappings

In Chapter 2 we studied solutions of the system $x=Ax+b$ that are fixed points
of the affine map $Fx=Ax+b$ on $\mathbbm{R}^{n}$ studied in Example 5.1.1. The
Neumann series lemma on page 1.2.5 is, in essence, a statement about existence
and uniqueness of fixed points of this map. Here we investigate another fixed
point theorem, due to Stefan Banach (1892–1945), that can be thought of as
extending the Neumann series lemma to nonlinear systems.

######  5.1.1.1 Contractions

Let $S$ be a nonempty subset of $\mathbbm{R}^{n}$. A self-map $F$ on $S$ is
called contracting or a contraction of modulus $\lambda$ on $S$ if there
exists a $\lambda<1$ and a norm $\|\cdot\|$ on $\mathbbm{R}^{n}$ such that

$\|Fu-Fv\|\leqslant\lambda\|u-v\|\quad\text{for all}\quad u,v\in S.$ (5.1)

Let $F$ be a contraction of modulus $\lambda$ on $S$. Show that

  1. (i)

$F$ is continuous on $S$ and

  2. (ii)

$F$ has at most one fixed point on $S$.

######  Example 5.1.1.

Let $S=\mathbbm{R}^{n}$, paired with the Euclidean norm $\|\cdot\|$. Let
$Fx=Ax+b$, where $A\in\mathbbm{M}^{n\times n}$ and $b\in\mathbbm{R}^{n}$. If
$\|A\|<1$, where $\|\cdot\|$ is the operator norm on $\mathbbm{M}^{n\times
n}$, then $F$ is a contraction of modulus $\|A\|$, since, for any $x,y\in S$,

$\|Ax+b-Ay-b\|=\|A(x-y)\|\leqslant\|A\|\|x-y\|.$

The next example uses a similar idea but is based on a different norm.

######  Example 5.1.2.

In (2.22) we studied the system $\rho=A^{\top}\rho-\varepsilon$, where
$A=(a_{ij})$ is a matrix of coefficients satisfying $\sum_{i}a_{ij}=1-\alpha$
for some $\alpha\in(0,1)$, the vector $\varepsilon$ is given and $\rho$ is
unknown. Solutions can be viewed as a fixed points of the map
$F\colon\mathbbm{R}^{n}\to\mathbbm{R}^{n}$ defined by
$Fp=A^{\top}p-\varepsilon$. Under the $\|\cdot\|_{\infty}$ norm, $F$ is a
contraction of modulus $1-\alpha$ on $\mathbbm{R}^{n}$. Indeed, for any
$p,q\in\mathbbm{R}^{n}$,

$\|Fp-
Fq\|_{\infty}=\max_{j}\left|\sum_{i=1}^{n}a_{ij}(p_{i}-q_{i})\right|\leqslant\max_{j}\sum_{i=1}^{n}a_{ij}\left|p_{i}-q_{i}\right|.$

Since $|p_{i}-q_{i}|\leqslant\|p-q\|_{\infty}$, we obtain

$\|Fp-
Fq\|_{\infty}\leqslant\max_{j}\sum_{i=1}^{n}a_{ij}\|p-q\|_{\infty}=(1-\alpha)\|p-q\|_{\infty}.$

Consider again Example 5.1.1. For the affine map $Fx=Ax+b$, the condition
$\|A\|<1$ used to obtain contraction is stronger than the condition $r(A)<1$
used to obtain a unique fixed point in the Neumann series lemma (see Exercise
2.3.2.3 on page 2.3.2.3). Furthermore, the Neumann series lemma provides a
geometric series representation of the fixed point. On the other hand, as we
now show, the contraction property can be used to obtain unique fixed points
when the map in question is not affine.

######  5.1.1.2 Banach’s Theorem

The fundamental importance of contractions stems from the following theorem.

######  Theorem 5.1.1 (Banach’s contraction mapping theorem).

If $S$ is closed in $\mathbbm{R}^{n}$ and $F$ is a contraction of modulus
$\lambda$ on $S$, then $F$ has a unique fixed point $u^{*}$ in $S$ and

$\|F^{n}u-u^{*}\|\leqslant\lambda^{n}\|u-u^{*}\|\quad\text{for all
}n\in\mathbbm{N}\text{ and }u\in S.$ (5.2)

In particular, $F$ is globally stable on $S$.

We complete a proof of Theorem 5.1.1 in stages. Let $S$ and $F$ have the
properties stated in Theorem 5.1.1. Fix $u_{0}\in S$ and let
$u_{m}:=F^{m}u_{0}$. Show that

$\|u_{m}-u_{k}\|\leqslant\sum_{i=m}^{k-1}\lambda^{i}\|u_{0}-u_{1}\|$

holds for all $m,k\in\mathbbm{N}$ with $m<k$. Using the results in Exercise
5.1.1.2, prove that $(u_{m})$ is a Cauchy sequence in $\mathbbm{R}^{n}$ (see
§6.1.3.2 for notes on the Cauchy property). From the bound in Exercise
5.1.1.2, we obtain

$\|u_{m}-u_{k}\|\leqslant\frac{\lambda^{m}}{1-\lambda}\lambda^{i}\|u_{0}-u_{1}\|\qquad(m,k\in\mathbbm{N}\text{
with }m<k).$

Hence $(u_{m})$ is Cauchy, as claimed. Using Exercise 5.1.1.2, argue that
$(u_{m})$ hence has a limit $u^{*}\in\mathbbm{R}^{n}$. Prove that $u^{*}\in
S$.

###### Proof of Theorem 5.1.1.

In the exercises we proved existence of a point $u^{*}\in S$ such that
$F^{m}u\to u^{*}$. The fact that $u^{*}$ is a fixed point of $F$ now follows
from Lemma 6.1.6 on page 6.1.6 and Exercise 5.1.1.1. Uniqueness is implied by
Exercise 5.1.1.1. The bound (5.2) follows from iteration on the contraction
inequality (5.1) while setting $v=u^{*}$. ∎

######  5.1.1.3 Eventual Contractions

Let $S$ be a nonempty subset of $\mathbbm{R}^{n}$. A self-map $F$ on $S$ is
called eventually contracting if there exists a $k\in\mathbbm{N}$ such that
$F^{k}$ is a contraction on $S$. Significantly, most of the conclusions of
Banach’s theorem carry over to the case where $F$ is eventually contracting.

######  Theorem 5.1.2.

Let $F$ be a self-map on $S\subset\mathbbm{R}^{n}$. If $S$ is closed and $F$
is eventually contracting, then $F$ is globally stable on $S$.

Prove Theorem 5.1.2.111Hint: Theorem 5.1.1 is self-improving: it implies this
seemingly stronger result. The proof is not trivial but see if you can get it
started. You might like to note that $F^{k}$ has unique fixed point $u^{*}$ in
$S$. (Why?) Now consider the fact that
$\|Fu^{*}-u^{*}\|=\|FF^{nk}u^{*}-u^{*}\|$ for all $n\in\mathbbm{N}$. Let $S$
be complete, let $F$ be a self-map on $S$ and let $F^{k}$ be a uniform
contraction. Let $u^{*}$ be the unique fixed point of $F^{k}$. Fix
$\varepsilon>0$. We can choose $n$ such that
$\|F^{nk}Fu^{*}-u^{*}\|<\varepsilon$. But then

$\|Fu^{*}-u^{*}\|=\|FF^{nk}u^{*}-u^{*}\|=\|F^{nk}Fu^{*}-u^{*}\|<\varepsilon.$

Since $\varepsilon$ was arbitrary we have $\|Fu^{*}-u^{*}\|=0$, implying that
$u^{*}$ is a fixed point of $F$. Regarding convergence, fix $u\in S$. Given
$n\in\mathbbm{N}$, there exist integers $j(n)$ and $i(n)$ such that
$n=j(n)k+i(n)$, and $j(n)\to\infty$ as $n\to\infty$. Hence

$\|F^{n}u-u^{*}\|=\|F^{j(n)k+i(n)}u-u^{*}\|=\|F^{j(n)k}F^{i(n)}u-u^{*}\|\to
0\qquad(n\to\infty),$

by the assumptions on $F^{k}$. Convergence implies uniqueness of the fixed
point (why?). There is a close connection between Theorem 5.1.2 and the
Neumann series lemma (NSL). If $S=\mathbbm{R}^{n}$ and $Fx=Ax+b$ with
$r(A)<1$, then the NSL implies a unique fixed point. We can also obtain this
result from Theorem 5.1.2, since, for any $k\in\mathbbm{N}$,

$\|F^{k}x-F^{k}y\|=\|A^{k}x-A^{k}y\|=\|A^{k}(x-y)\|\leqslant\|A^{k}\|\|x-y\|.$

As $r(A)<1$, we can choose $k$ such that $\|A^{k}\|<1$ (see §2.3.3). Hence $F$
is eventually contracting and Theorem 5.1.2 applies. As mentioned above,
contractions and eventual contractions have much wider scope than the NSL,
since they can also be applied in nonlinear settings. At the same time, the
NSL is preferred when its conditions hold, since it also gives inverse and
power series representations of the fixed point.

######  5.1.1.4 A Condition for Eventual Contractions

The result below provides a useful test for the eventual contraction property.
(In the statement, the absolute value of a vector is defined pointwise, as in
§6.1.2.2.)

######  Proposition 5.1.3.

Let $F$ be a self-map on $S\subset\mathbbm{R}^{n}$ such that, for some
$n\times n$ matrix $A$,

$|Fx-Fy|\leqslant A|x-y|\quad\text{for all }\,x,y\in S.$

If, in addition, $A\geqslant 0$ and $r(A)<1$, then $F$ is eventually
contracting on $S$ with respect to the Euclidean norm.

###### Proof.

Our first claim is that, under the conditions of the proposition,

$|F^{k}x-F^{k}y|\leqslant A^{k}|x-y|\quad\text{for all $k\in\mathbbm{N}$ and
$x,y\in S$}.$ (5.3)

This is true at $k=1$ by assumption. If it is true at $k-1$, then

$|F^{k}x-F^{k}y|\leqslant A|F^{k-1}x-F^{k-1}y|\leqslant AA^{k-1}|x-y|,$ (5.4)

where the second inequality uses the induction hypothesis and $A\geqslant 0$
(so that $u\leqslant v$ implies $Au\leqslant Av$). Hence (5.3) holds. It
follows from the definition of the Euclidean norm that $\||u|\|=\|u\|$ for any
vector $u$. Also, for the same norm, $|u|\leqslant|v|$ implies
$\|u\|\leqslant\|v\|$. Hence, for all $k\in\mathbbm{N}$ and $x,y\in S$,

$\|F^{k}x-F^{k}y\|\leqslant\|A^{k}|x-y|\|\leqslant\|A^{k}\|_{o}\|x-y\|.$

In the second inequality, we used $\|\cdot\|_{o}$ for the operator norm,
combined with the fact that $\|Au\|\leqslant\|A\|_{o}\|u\|$ always holds, as
discussed in §2.3.2.3. By Gelfand’s lemma (see in particular Exercise 2.3.3.1
on page 2.3.3.1), we obtain existence of a $\lambda\in(0,1)$ and a
$k\in\mathbbm{N}$ with $\|A^{k}\|_{o}\leqslant\lambda<1$. Hence, for this $k$,

$\|F^{k}x-F^{k}y\|\leqslant\lambda\|x-y\|.$

Since $\lambda$ does not depend on $x$ or $y$, we have shown that $F$ is an
eventual contraction on $S$ with respect to the Euclidean norm. ∎

#####  5.1.2 Shortest Paths Revisited

Consider again the shortest path problem introduced in §3.1. One modification
that sometimes appears in applications is the addition of discounting during
travel between vertices. For example, if the vertices are international ports
and travel takes place by sea, then port-to-port travel time is measured in
weeks or even months. It is natural to apply time discounting to future costs
associated with that travel, to implement the idea that paying a given dollar
amount in the future is preferable to paying it now. Suppose it is possible to
borrow and lend risk-free at a positive interest rate $r$. Explain why it is
always preferable to have $100 now than $100 in a year’s time in this setting.
If the risk-free real interest rate $r$ is positive, then $100 received now
can be converted with probability one into $(1+r)100>100$ dollars in one year.
Recall that, without discounting, the Bellman equation for the shortest path
problem takes the form $q(x)=\min_{y\in\mathscr{O}(x)}\\{c(x,y)+q(y)\\}$ for
all $x\in V$, where $c$ is the cost function $V$ is the set of vertices and
$q$ is a candidate for the cost-to-go function. We showed that the minimum
cost-to-go function $q^{*}$ satisfies the Bellman equation and is the unique
fixed point of the Bellman operator. The Bellman equation neatly divides the
problem into current costs, embedded in the term $c(x,y)$, and future costs
embedded in $q(y)$. To add discounting, we need only discount $q(y)$. We do
this by multiplying it by a discount factor $\beta\in(0,1)$. The Bellman
equation is then $q(x)=\min_{y\in\mathscr{O}(x)}\\{c(x,y)+\beta q(y)\\}$ for
all $x\in V$ and the Bellman operator is

$Tq(x)=\min_{y\in\mathscr{O}(x)}\\{c(x,y)+\beta q(y)\\}\qquad(x\in V).$ (5.5)

In §3.1.2, without discounting, we had to work hard to show that the Bellman
operator has a unique fixed point in $U$, the set of all $q\colon
V\to\mathbbm{R}_{+}$ with $q(d)=0$. With discounting the proof is easier,
since we can leverage the Banach contraction mapping theorem. In what follows,
we identify the vertices in $V$ with integers $1,\ldots,n$, where $d$ is
identified with $n$. We then understand $U$ as all nonnegative vectors $q$ in
$\mathbbm{R}^{n}$ with $q(n)=0$. (We continue to write $q(x)$ for the $x$-th
element of the vector $q$, but now $x$ is in $[n]$.) Prove that $U$ is a
closed subset of $\mathbbm{R}^{n}$. Take $q_{k}\to q$ where $(q_{k})$ is a
sequence of $n$-vectors contained in $U$. By Exercise 6.1.2.1 on page 6.1.2.1,
since $q_{k}\geqslant 0$ for all $k$, we must have $q\geqslant 0$. It remains
only to show that $q(n)=0$. As $q_{k}\in U$ for all $k$, we have $q_{k}(n)=0$
for all $k$. By Lemma 2.3.1, we also have $q_{k}(n)\to q(n)$. Hence $q(n)=0$.
Prove that $T$ is order-preserving on $U$ with respect to the pointwise order.
We need to show that if $p,q\in U$ and $p\leqslant q$, then $Tp\leqslant Tq$.
This follows easily from the definition of $T$ in (5.5). Prove that, for any
$q\in U$ and $\alpha\in\mathbbm{R}_{+}$, we have
$T(q+\alpha\mathbbm{1})=Tq+\beta\alpha\mathbbm{1}$. Fix $q\in U$
$\alpha\in\mathbbm{R}_{+}$ and $x\in V$. By definition,

$T(q+\alpha\mathbbm{1})(x)=\min_{y\in\mathscr{O}(x)}\\{c(x,y)+\beta
q(y)+\alpha\beta\\}=Tq(x)+\alpha\beta.$

Hence $T(q+\alpha\mathbbm{1})=Tq+\beta\alpha\mathbbm{1}$ as claimed. Now let
$\|\cdot\|_{\infty}$ be the supremum norm on $\mathbbm{R}^{n}$ (see §2.3.1).
We claim that $T$ is a contraction on $U$ of modulus $\beta$. To see that this
is so, fix $p,q\in U$ and observe that, pointwise,

$Tq=T(p+q-p)\leqslant T(p+\|q-p\|_{\infty}\mathbbm{1})\leqslant
Tp+\beta\|q-p\|_{\infty}\mathbbm{1},$

where the first inequality is by the order-preserving of $T$ and the second
follows from Exercise 5.1.2. Hence

$Tq-Tp\leqslant\beta\|q-p\|_{\infty}\mathbbm{1}.$

Reversing the roles of $p$ and $q$ gives the reverse inequality. Hence

$|Tq(x)-Tp(x)|\leqslant\beta\|q-p\|_{\infty}$

for all $x\in[n]$. Taking the maximum on the left hand side yields $\|Tq-
Tp\|_{\infty}\leqslant\beta\|q-p\|_{\infty}$, which shows that $T$ is a
contraction of modulus $\beta$. Hence Banach’s theorem applies and a unique
fixed point exists.

#####  5.1.3 Supply Constraints

While the input-output model from §2.1.2 has many useful applications, its
linear structure can be a liability. One natural objection to linearity is
supply constraints: if sector $j$ doubles its orders from sector $i$, we
cannot always expect that sector $i$ will quickly meet this jump in demand. In
this sector we investigate the impact of supply constraints on equilibrium.
These constraints introduce nonlinear relationships between nodes that affect
equilibria and make analysis more challenging.

######  5.1.3.1 Production with Constraints

We recall from §2.1.2 that $d_{i}$ is final demand for good $i$, $x_{i}$ is
total sales of sector $i$, $z_{ij}$ is inter-industry sales from sector $i$ to
sector $j$, and $a_{ij}=z_{ij}/x_{j}$ is dollar value of inputs from $i$ per
dollar output from $j$. Departing from our previous formulation of equilibrium
in the input-output model, suppose that, in the short run, the total output
value of sector $i$ is constrained by positive constant $\bar{x}(i)$. Holding
prices fixed (in the short run), this means that sector $i$ has a capacity
constraint in terms of unit output. For the purposes of our model, the vector
of capacity constraints $\bar{x}:=(\bar{x}(i))_{i=1}^{n}$ can be any vector in
$\mathbbm{R}^{n}_{+}$. For each sector $i$, we modify the accounting identity
(2.1) from page 2.1 to

$x_{i}=\min\left\\{\sum_{j=1}^{n}o_{ij}+d_{i},\;\bar{x}(i)\right\\},$ (5.6)

where $o_{ij}$ is the value of orders from sector $i$ made by sector $j$.
Thus, if the capacity constraint in sector $i$ is not binding, then output is
the sum of orders from other sectors and orders from final consumers. If
$\bar{x}(i)$ is less than this sum, however, then sector $i$ produces to
capacity $\bar{x}(i)$. An equilibrium for this model is one where all orders