Skip to content

Commit

Permalink
[markov_chain_I] latest edits (#510)
Browse files Browse the repository at this point in the history 
* latest edits to mc I lecture

* misc

* misc

* update the label

there are two labels for this theorem.

* Revert "update the label"

This reverts commit be885de.

* update label

update the ref to the theorem

---------

Co-authored-by: Longye Tian <[email protected]>
  • Loading branch information
@jstac @longye-tian
jstac and longye-tian authored Jul 8, 2024
1 parent fbe4c82 commit 4c08c22
Showing 1 changed file with 87 additions and 111 deletions.
198 changes: 87 additions & 111 deletions lectures/markov_chains_I.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,8 @@ import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
```

## Definitions and examples
Expand Down Expand Up @@ -743,6 +745,11 @@ This is, in some sense, a steady state probability of unemployment.

Not surprisingly it tends to zero as $\beta \to 0$, and to one as $\alpha \to 0$.






### Calculating stationary distributions

A stable algorithm for computing stationary distributions is implemented in [QuantEcon.py](http://quantecon.org/quantecon-py).
Expand All @@ -757,6 +764,11 @@ mc = qe.MarkovChain(P)
mc.stationary_distributions # Show all stationary distributions
```






### Asymptotic stationarity

Consider an everywhere positive stochastic matrix with unique stationary distribution $\psi^*$.
Expand All @@ -767,17 +779,24 @@ For example, we have the following result

(strict_stationary)=
```{prf:theorem}
:label: mc_gs_thm
If there exists an integer $m$ such that all entries of $P^m$ are
strictly positive, with unique stationary distribution $\psi^*$, then
strictly positive, then
$$
\psi_0 P^t \to \psi^*
\quad \text{ as } t \to \infty
$$
where $\psi^*$ is the unique stationary distribution.
```

This situation is often referred to as **asymptotic stationarity** or **global stability**.

A proof of the theorem can be found in Chapter 4 of {cite}`sargent2023economic`, as well as many other sources.


See, for example, {cite}`sargent2023economic` Chapter 4.



Expand Down Expand Up @@ -848,103 +867,96 @@ Here
You might like to try experimenting with different initial conditions.


#### An alternative illustration

We can show this in a slightly different way by focusing on the probability that $\psi_t$ puts on each state.

First, we write a function to draw initial distributions $\psi_0$ of size `num_distributions`

```{code-cell} ipython3
def generate_initial_values(num_distributions):
n = len(P)
draws = np.random.randint(1, 10_000_000, size=(num_distributions,n))
ψ_0s = draws/draws.sum(axis=1)[:, None]
return ψ_0s
```
#### Example: failure of convergence

We then write a function to plot the dynamics of $(\psi_0 P^t)(i)$ as $t$ gets large, for each state $i$ with different initial distributions

```{code-cell} ipython3
def plot_distribution(P, ts_length, num_distributions):
Consider the periodic chain with stochastic matrix

# Get parameters of transition matrix
n = len(P)
mc = qe.MarkovChain(P)
ψ_star = mc.stationary_distributions[0]
$$
P =
\begin{bmatrix}
0 & 1 \\
1 & 0 \\
\end{bmatrix}
$$

## Draw the plot
fig, axes = plt.subplots(nrows=1, ncols=n, figsize=[11, 5])
plt.subplots_adjust(wspace=0.35)
This matrix does not satisfy the conditions of
{ref}`strict_stationary` because, as you can readily check,

ψ_0s = generate_initial_values(num_distributions)
* $P^m = P$ when $m$ is odd and
* $P^m = I$, the identity matrix, when $m$ is even.

# Get the path for each starting value
for ψ_0 in ψ_0s:
ψ_t = iterate_ψ(ψ_0, P, ts_length)
Hence there is no $m$ such that all elements of $P^m$ are strictly positive.

# Obtain and plot distributions at each state
for i in range(n):
axes[i].plot(range(0, ts_length), ψ_t[:,i], alpha=0.3)
Moreover, we can see that global stability does not hold.

# Add labels
for i in range(n):
axes[i].axhline(ψ_star[i], linestyle='dashed', lw=2, color = 'black',
label = fr'$\psi^*({i})$')
axes[i].set_xlabel('t')
axes[i].set_ylabel(fr'$\psi_t({i})$')
axes[i].legend()
For instance, if we start at $\psi_0 = (1,0)$, then $\psi_m = \psi_0 P^m$ is $(1, 0)$ when $m$ is even and $(0,1)$ when $m$ is odd.

plt.show()
```
We can see similar phenomena in higher dimensions.

The following figure shows
The next figure illustrates this for a periodic Markov chain with three states.

```{code-cell} ipython3
# Define the number of iterations
# and initial distributions
ts_length = 50
num_distributions = 25
P = np.array([[0.971, 0.029, 0.000],
[0.145, 0.778, 0.077],
[0.000, 0.508, 0.492]])
plot_distribution(P, ts_length, num_distributions)
```

The convergence to $\psi^*$ holds for different initial distributions.
ψ_1 = (0.0, 0.0, 1.0)
ψ_2 = (0.5, 0.5, 0.0)
ψ_3 = (0.25, 0.25, 0.5)
ψ_4 = (1/3, 1/3, 1/3)
P = np.array([[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
[1.0, 0.0, 0.0]])
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
colors = ['red','yellow', 'green', 'blue'] # Different colors for each initial point
#### Example: failure of convergence

# Define the vertices of the unit simplex
v = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]])
In the case of a periodic chain, with
# Define the faces of the unit simplex
faces = [
[v[0], v[1], v[2]],
[v[0], v[1], v[3]],
[v[0], v[2], v[3]],
[v[1], v[2], v[3]]
]
$$
P =
\begin{bmatrix}
0 & 1 \\
1 & 0 \\
\end{bmatrix}
$$
def update(n):
ax.clear()
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])
ax.set_zlim([0, 1])
ax.view_init(45, 45)
# Plot the 3D unit simplex as planes
simplex = Poly3DCollection(faces,alpha=0.05)
ax.add_collection3d(simplex)
for idx, ψ_0 in enumerate([ψ_1, ψ_2, ψ_3, ψ_4]):
ψ_t = iterate_ψ(ψ_0, P, n+1)
point = ψ_t[-1]
ax.scatter(point[0], point[1], point[2], color=colors[idx], s=60)
points = np.array(ψ_t)
ax.plot(points[:, 0], points[:, 1], points[:, 2], color=colors[idx],linewidth=0.75)
return fig,
we find the distribution oscillates
anim = FuncAnimation(fig, update, frames=range(20), blit=False, repeat=False)
plt.close()
HTML(anim.to_jshtml())
```
This animation demonstrates the behavior of an irreducible and periodic stochastic matrix.

```{code-cell} ipython3
P = np.array([[0, 1],
[1, 0]])
The red, yellow, and green dots represent different initial probability distributions.

ts_length = 20
num_distributions = 30
The blue dot represents the unique stationary distribution.

plot_distribution(P, ts_length, num_distributions)
```
Unlike Hamilton’s Markov chain, these initial distributions do not converge to the unique stationary distribution.

Indeed, this $P$ fails our asymptotic stationarity condition, since, as you can
verify, $P^t$ is not everywhere positive for any $t$.
Instead, they cycle periodically around the probability simplex, illustrating that asymptotic stability fails.


(finite_mc_expec)=
Expand Down Expand Up @@ -1105,42 +1117,6 @@ mc = qe.MarkovChain(P)
ψ_star
```

Solution 3:

We find the distribution $\psi$ converges to the stationary distribution more quickly compared to the {ref}`hamilton's chain <hamilton>`.

```{code-cell} ipython3
ts_length = 10
num_distributions = 25
plot_distribution(P, ts_length, num_distributions)
```

In fact, the rate of convergence is governed by {ref}`eigenvalues<eigen>` {cite}`sargent2023economic`.

```{code-cell} ipython3
P_eigenvals = np.linalg.eigvals(P)
P_eigenvals
```

```{code-cell} ipython3
P_hamilton = np.array([[0.971, 0.029, 0.000],
[0.145, 0.778, 0.077],
[0.000, 0.508, 0.492]])
hamilton_eigenvals = np.linalg.eigvals(P_hamilton)
hamilton_eigenvals
```

More specifically, it is governed by the spectral gap, the difference between the largest and the second largest eigenvalue.

```{code-cell} ipython3
sp_gap_P = P_eigenvals[0] - np.diff(P_eigenvals)[0]
sp_gap_hamilton = hamilton_eigenvals[0] - np.diff(hamilton_eigenvals)[0]
sp_gap_P > sp_gap_hamilton
```

We will come back to this when we discuss {ref}`spectral theory<spec_markov>`.

```{solution-end}
```
Expand Down

1 comment on commit 4c08c22

@github-actions
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🚀 Deployed on https://668bb425101d2c16a9ef2641--taupe-gaufre-c4e660.netlify.app

Please sign in to comment.