@@ -11,7 +11,7 @@ kernelspec:
11
11
name : python3
12
12
---
13
13
14
- +++ {"user_expressions": [ ] }
14
+
15
15
16
16
# Markov Chains: Basic Concepts and Stationarity
17
17
@@ -36,10 +36,10 @@ In addition to what's in Anaconda, this lecture will need the following librarie
36
36
:class: warning
37
37
If you are running this lecture locally it requires [graphviz](https://www.graphviz.org)
38
38
to be installed on your computer. Installation instructions for graphviz can be found
39
- [here](https://www.graphviz.org/download/)
39
+ [here](https://www.graphviz.org/download/)
40
40
```
41
41
42
- +++ {"user_expressions": [ ] }
42
+
43
43
44
44
## Overview
45
45
@@ -74,7 +74,7 @@ import matplotlib as mpl
74
74
from itertools import cycle
75
75
```
76
76
77
- +++ {"user_expressions": [ ] }
77
+
78
78
79
79
## Definitions and examples
80
80
@@ -137,7 +137,7 @@ dot.edge("sr", "sr", label="0.492")
137
137
dot
138
138
```
139
139
140
- +++ {"user_expressions": [ ] }
140
+
141
141
142
142
Here there are three ** states**
143
143
@@ -180,11 +180,11 @@ We can collect all of these conditional probabilities into a matrix, as follows
180
180
181
181
$$
182
182
P =
183
- \begin{bmatrix}
183
+ \begin{bmatrix}
184
184
0.971 & 0.029 & 0 \\
185
185
0.145 & 0.778 & 0.077 \\
186
186
0 & 0.508 & 0.492
187
- \end{bmatrix}
187
+ \end{bmatrix}
188
188
$$
189
189
190
190
Notice that $P$ is a stochastic matrix.
@@ -220,11 +220,11 @@ Given the above information, we can write out the transition probabilities in ma
220
220
``` {math}
221
221
:label: p_unempemp
222
222
223
- P =
224
- \begin{bmatrix}
223
+ P =
224
+ \begin{bmatrix}
225
225
1 - \alpha & \alpha \\
226
226
\beta & 1 - \beta
227
- \end{bmatrix}
227
+ \end{bmatrix}
228
228
```
229
229
230
230
For example,
@@ -255,26 +255,26 @@ We'll cover some of these applications below.
255
255
(mc_eg3)=
256
256
#### Example 3
257
257
258
- Imam and Temple {cite}` imampolitical ` categorize political institutions into three types: democracy (D), autocracy (A), and an intermediate state called anocracy (N).
258
+ Imam and Temple {cite}` imampolitical ` categorize political institutions into three types: democracy (D), autocracy (A), and an intermediate state called anocracy (N).
259
259
260
- Each institution can have two potential development regimes: collapse (C) and growth (G). This results in six possible states: DG, DC, NG, NC, AG, and AC.
260
+ Each institution can have two potential development regimes: collapse (C) and growth (G). This results in six possible states: DG, DC, NG, NC, AG, and AC.
261
261
262
- The lower probability of transitioning from NC to itself indicates that collapses in anocracies quickly evolve into changes in the political institution.
262
+ The lower probability of transitioning from NC to itself indicates that collapses in anocracies quickly evolve into changes in the political institution.
263
263
264
264
Democracies tend to have longer-lasting growth regimes compared to autocracies as indicated by the lower probability of transitioning from growth to growth in autocracies.
265
265
266
266
We can also find a higher probability from collapse to growth in democratic regimes
267
267
268
268
$$
269
269
P :=
270
- \begin{bmatrix}
270
+ \begin{bmatrix}
271
271
0.86 & 0.11 & 0.03 & 0.00 & 0.00 & 0.00 \\
272
272
0.52 & 0.33 & 0.13 & 0.02 & 0.00 & 0.00 \\
273
273
0.12 & 0.03 & 0.70 & 0.11 & 0.03 & 0.01 \\
274
274
0.13 & 0.02 & 0.35 & 0.36 & 0.10 & 0.04 \\
275
275
0.00 & 0.00 & 0.09 & 0.11 & 0.55 & 0.25 \\
276
276
0.00 & 0.00 & 0.09 & 0.15 & 0.26 & 0.50
277
- \end{bmatrix}
277
+ \end{bmatrix}
278
278
$$
279
279
280
280
``` {code-cell} ipython3
@@ -297,7 +297,7 @@ for start_idx, node_start in enumerate(nodes):
297
297
value = P[start_idx][end_idx]
298
298
if value != 0:
299
299
G.add_edge(node_start,node_end, weight=value, len=100)
300
-
300
+
301
301
pos = nx.spring_layout(G, seed=10)
302
302
fig, ax = plt.subplots()
303
303
nx.draw_networkx_nodes(G, pos, node_size=600, edgecolors='black', node_color='white')
@@ -327,7 +327,7 @@ The set $S$ is called the **state space** and $x_1, \ldots, x_n$ are the **state
327
327
328
328
A ** distribution** $\psi$ on $S$ is a probability mass function of length $n$, where $\psi(i)$ is the amount of probability allocated to state $x_i$.
329
329
330
- A ** Markov chain** $\{ X_t\} $ on $S$ is a sequence of random variables taking values in $S$
330
+ A ** Markov chain** $\{ X_t\} $ on $S$ is a sequence of random variables taking values in $S$
331
331
that have the ** Markov property** .
332
332
333
333
This means that, for any date $t$ and any state $y \in S$,
@@ -390,9 +390,9 @@ In these exercises, we'll take the state space to be $S = 0,\ldots, n-1$.
390
390
391
391
### Writing our own simulation code
392
392
393
- To simulate a Markov chain, we need
393
+ To simulate a Markov chain, we need
394
394
395
- 1 . a stochastic matrix $P$ and
395
+ 1 . a stochastic matrix $P$ and
396
396
1 . a probability mass function $\psi_0$ of length $n$ from which to draw a initial realization of $X_0$.
397
397
398
398
The Markov chain is then constructed as follows:
@@ -416,7 +416,7 @@ cdf = np.cumsum(ψ_0) # convert into cumulative distribution
416
416
qe.random.draw(cdf, 5) # generate 5 independent draws from ψ
417
417
```
418
418
419
- +++ {"user_expressions": [ ] }
419
+
420
420
421
421
We'll write our code as a function that accepts the following three arguments
422
422
@@ -449,7 +449,7 @@ def mc_sample_path(P, ψ_0=None, ts_length=1_000):
449
449
return X
450
450
```
451
451
452
- +++ {"user_expressions": [ ] }
452
+
453
453
454
454
Let's see how it works using the small matrix
455
455
@@ -458,15 +458,15 @@ P = [[0.4, 0.6],
458
458
[0.2, 0.8]]
459
459
```
460
460
461
- +++ {"user_expressions": [ ] }
461
+
462
462
463
463
Here's a short time series.
464
464
465
465
``` {code-cell} ipython3
466
466
mc_sample_path(P, ψ_0=[1.0, 0.0], ts_length=10)
467
467
```
468
468
469
- +++ {"user_expressions": [ ] }
469
+
470
470
471
471
It can be shown that for a long series drawn from ` P ` , the fraction of the
472
472
sample that takes value 0 will be about 0.25.
@@ -483,7 +483,7 @@ X = mc_sample_path(P, ψ_0=[0.1, 0.9], ts_length=1_000_000)
483
483
np.mean(X == 0)
484
484
```
485
485
486
- +++ {"user_expressions": [ ] }
486
+
487
487
488
488
You can try changing the initial distribution to confirm that the output is
489
489
always close to 0.25 (for the ` P ` matrix above).
@@ -501,7 +501,7 @@ X = mc.simulate(ts_length=1_000_000)
501
501
np.mean(X == 0)
502
502
```
503
503
504
- +++ {"user_expressions": [ ] }
504
+
505
505
506
506
The ` simulate ` routine is faster (because it is [ JIT compiled] ( https://python-programming.quantecon.org/numba.html#numba-link ) ).
507
507
@@ -513,7 +513,7 @@ The `simulate` routine is faster (because it is [JIT compiled](https://python-pr
513
513
%time mc.simulate(ts_length=1_000_000) # qe code version
514
514
```
515
515
516
- +++ {"user_expressions": [ ] }
516
+
517
517
518
518
#### Adding state values and initial conditions
519
519
@@ -536,15 +536,15 @@ mc.simulate(ts_length=4, init='unemployed')
536
536
mc.simulate(ts_length=4) # Start at randomly chosen initial state
537
537
```
538
538
539
- +++ {"user_expressions": [ ] }
539
+
540
540
541
541
If we want to see indices rather than state values as outputs as we can use
542
542
543
543
``` {code-cell} ipython3
544
544
mc.simulate_indices(ts_length=4)
545
545
```
546
546
547
- +++ {"user_expressions": [ ] }
547
+
548
548
549
549
(mc_md)=
550
550
## Distributions over time
@@ -621,7 +621,7 @@ Hence the following is also valid.
621
621
X_t \sim \psi_t \quad \implies \quad X_{t+m} \sim \psi_t P^m
622
622
```
623
623
624
- +++ {"user_expressions": [ ] }
624
+
625
625
626
626
(finite_mc_mstp)=
627
627
### Multiple step transition probabilities
@@ -657,20 +657,20 @@ Suppose that the current state is unknown --- perhaps statistics are available o
657
657
658
658
We guess that the probability that the economy is in state $x$ is $\psi_t(x)$ at time t.
659
659
660
- The probability of being in recession (either mild or severe) in 6 months time is given by
660
+ The probability of being in recession (either mild or severe) in 6 months time is given by
661
661
662
662
$$
663
663
(\psi_t P^6)(1) + (\psi_t P^6)(2)
664
664
$$
665
665
666
- +++ {"user_expressions": [ ] }
666
+
667
667
668
668
(mc_eg1-1)=
669
669
### Example 2: Cross-sectional distributions
670
670
671
- The distributions we have been studying can be viewed either
671
+ The distributions we have been studying can be viewed either
672
672
673
- 1 . as probabilities or
673
+ 1 . as probabilities or
674
674
1 . as cross-sectional frequencies that a Law of Large Numbers leads us to anticipate for large samples.
675
675
676
676
To illustrate, recall our model of employment/unemployment dynamics for a given worker {ref}` discussed above <mc_eg1> ` .
@@ -720,11 +720,11 @@ P = np.array([[0.4, 0.6],
720
720
ψ @ P
721
721
```
722
722
723
- +++ {"user_expressions": [ ] }
723
+
724
724
725
725
Notice that ` ψ @ P ` is the same as ` ψ `
726
726
727
- +++ {"user_expressions": [ ] }
727
+
728
728
729
729
Such distributions are called ** stationary** or ** invariant** .
730
730
@@ -765,7 +765,7 @@ distribution.
765
765
766
766
We will come back to this when we introduce irreducibility in the next lecture
767
767
768
- +++ {"user_expressions": [ ] }
768
+
769
769
770
770
### Example
771
771
822
822
823
823
See, for example, {cite}` sargent2023economic ` Chapter 4.
824
824
825
- +++ {"user_expressions": [ ] }
825
+
826
826
827
827
(hamilton)=
828
828
#### Example: Hamilton's chain
@@ -836,7 +836,7 @@ P = np.array([[0.971, 0.029, 0.000],
836
836
P @ P
837
837
```
838
838
839
- +++ {"user_expressions": [ ] }
839
+
840
840
841
841
Let's pick an initial distribution $\psi_0$ and trace out the sequence of distributions $\psi_0 P^t$ for $t = 0, 1, 2, \ldots$
842
842
@@ -900,13 +900,13 @@ First, we write a function to draw initial distributions $\psi_0$ of size `num_d
900
900
def generate_initial_values(num_distributions):
901
901
n = len(P)
902
902
ψ_0s = np.empty((num_distributions, n))
903
-
903
+
904
904
for i in range(num_distributions):
905
905
draws = np.random.randint(1, 10_000_000, size=n)
906
906
907
907
# Scale them so that they add up into 1
908
908
ψ_0s[i,:] = np.array(draws/sum(draws))
909
-
909
+
910
910
return ψ_0s
911
911
```
912
912
@@ -929,14 +929,14 @@ def plot_distribution(P, ts_length, num_distributions):
929
929
# Get the path for each starting value
930
930
for ψ_0 in ψ_0s:
931
931
ψ_t = iterate_ψ(ψ_0, P, ts_length)
932
-
932
+
933
933
# Obtain and plot distributions at each state
934
934
for i in range(n):
935
935
axes[i].plot(range(0, ts_length), ψ_t[:,i], alpha=0.3)
936
936
937
937
# Add labels
938
938
for i in range(n):
939
- axes[i].axhline(ψ_star[i], linestyle='dashed', lw=2, color = 'black',
939
+ axes[i].axhline(ψ_star[i], linestyle='dashed', lw=2, color = 'black',
940
940
label = fr'$\psi^*({i})$')
941
941
axes[i].set_xlabel('t')
942
942
axes[i].set_ylabel(fr'$\psi_t({i})$')
@@ -948,7 +948,7 @@ def plot_distribution(P, ts_length, num_distributions):
948
948
The following figure shows
949
949
950
950
``` {code-cell} ipython3
951
- # Define the number of iterations
951
+ # Define the number of iterations
952
952
# and initial distributions
953
953
ts_length = 50
954
954
num_distributions = 25
@@ -962,7 +962,7 @@ plot_distribution(P, ts_length, num_distributions)
962
962
963
963
The convergence to $\psi^* $ holds for different initial distributions.
964
964
965
- +++ {"user_expressions": [ ] }
965
+
966
966
967
967
#### Example: Failure of convergence
968
968
@@ -979,7 +979,7 @@ num_distributions = 30
979
979
plot_distribution(P, ts_length, num_distributions)
980
980
```
981
981
982
- +++ {"user_expressions": [ ] }
982
+
983
983
984
984
(finite_mc_expec)=
985
985
## Computing expectations
@@ -1010,8 +1010,8 @@ where
1010
1010
algebra, we'll think of as the column vector
1011
1011
1012
1012
$$
1013
- h =
1014
- \begin{bmatrix}
1013
+ h =
1014
+ \begin{bmatrix}
1015
1015
h(x_1) \\
1016
1016
\vdots \\
1017
1017
h(x_n)
@@ -1058,15 +1058,15 @@ $\sum_t \beta^t h(X_t)$.
1058
1058
In view of the preceding discussion, this is
1059
1059
1060
1060
$$
1061
- \mathbb{E}
1061
+ \mathbb{E}
1062
1062
\left[
1063
- \sum_{j=0}^\infty \beta^j h(X_{t+j}) \mid X_t
1063
+ \sum_{j=0}^\infty \beta^j h(X_{t+j}) \mid X_t
1064
1064
= x
1065
1065
\right]
1066
1066
= x + \beta (Ph)(x) + \beta^2 (P^2 h)(x) + \cdots
1067
1067
$$
1068
1068
1069
- By the {ref}` Neumann series lemma <la_neumann> ` , this sum can be calculated using
1069
+ By the {ref}` Neumann series lemma <la_neumann> ` , this sum can be calculated using
1070
1070
1071
1071
$$
1072
1072
I + \beta P + \beta^2 P^2 + \cdots = (I - \beta P)^{-1}
@@ -1082,11 +1082,11 @@ Imam and Temple {cite}`imampolitical` used a three-state transition matrix to de
1082
1082
1083
1083
$$
1084
1084
P :=
1085
- \begin{bmatrix}
1085
+ \begin{bmatrix}
1086
1086
0.68 & 0.12 & 0.20 \\
1087
1087
0.50 & 0.24 & 0.26 \\
1088
1088
0.36 & 0.18 & 0.46
1089
- \end{bmatrix}
1089
+ \end{bmatrix}
1090
1090
$$
1091
1091
1092
1092
where rows, from top to down, correspond to growth, stagnation, and collapse.
0 commit comments