diff --git a/lectures/_static/networks/mc.png b/lectures/_static/networks/mc.png
new file mode 100644
index 0000000..caeccf2
Binary files /dev/null and b/lectures/_static/networks/mc.png differ
diff --git a/lectures/_static/networks/poverty_trap_1.png b/lectures/_static/networks/poverty_trap_1.png
new file mode 100644
index 0000000..e3bda2a
Binary files /dev/null and b/lectures/_static/networks/poverty_trap_1.png differ
diff --git a/lectures/_static/networks/poverty_trap_2.png b/lectures/_static/networks/poverty_trap_2.png
new file mode 100644
index 0000000..71f2d31
Binary files /dev/null and b/lectures/_static/networks/poverty_trap_2.png differ
diff --git a/lectures/_static/networks/properties.png b/lectures/_static/networks/properties.png
new file mode 100644
index 0000000..7098688
Binary files /dev/null and b/lectures/_static/networks/properties.png differ
diff --git a/lectures/_static/networks/weighted.png b/lectures/_static/networks/weighted.png
new file mode 100644
index 0000000..5510b4d
Binary files /dev/null and b/lectures/_static/networks/weighted.png differ
diff --git a/lectures/information_market_equilibrium.md b/lectures/information_market_equilibrium.md
new file mode 100644
index 0000000..2fb1f92
--- /dev/null
+++ b/lectures/information_market_equilibrium.md
@@ -0,0 +1,1574 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+ jupytext_version: 1.17.1
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+---
+
+(information_market_equilibrium)=
+```{raw} jupyter
+
+```
+
+# Information and Market Equilibrium
+
+```{contents} Contents
+:depth: 2
+```
+
+## Overview
+
+This lecture studies two questions about the **informational role of prices**
+posed and
+answered by {cite:t}`kihlstrom_mirman1975`.
+
+1. *When do prices transmit inside information?*
+ - An informed insider observes a private
+ signal correlated with an unknown state of the world and adjusts demand
+ accordingly.
+ - Equilibrium prices shift.
+ - Under what conditions can an outside observer *infer* the insider's
+ posterior distribution from the equilibrium price?
+
+2. *Do Bayesian price expectations converge?*
+ - In a stationary stochastic exchange
+ economy, an uninformed observer uses the history of market prices and
+ Bayes' Law to form
+ beliefs about the economy's structure and hence about its induced price
+ distribution.
+ - Do those expectations eventually
+ agree with those of a fully informed observer?
+
+Kihlstrom and Mirman's answers rely on two classical ideas from statistics:
+
+- **Blackwell sufficiency**: a random variable $\tilde{y}$ is said to be
+ *sufficient* for a random variable
+ $\tilde{y}'$ with respect to an unknown state if knowing $\tilde{y}$ gives
+ all the
+ information about the state that $\tilde{y}'$ contains.
+- **Bayesian consistency**: as the sample grows, posterior beliefs eliminate
+ models that imply the wrong **price distribution**, so even when structure is
+ not identified from prices the posterior mass on the true **reduced form**
+ still converges to one.
+
+Important findings of {cite:t}`kihlstrom_mirman1975` are:
+
+- Equilibrium prices transmit inside information *if and only if* the map from
+ the
+ insider's posterior distribution to the equilibrium price is one-to-one on
+ the set of
+ posteriors that can actually arise from the signal.
+ - For the two-state case ($S = 2$), invertibility holds when the informed
+ agent's utility is homothetic and the elasticity of substitution is everywhere
+ either below one or above one.
+- In the dynamic economy, as information accumulates, Bayesian price
+ expectations converge to **rational expectations**, even when the deep
+ structure is not identified from prices alone.
+
+```{note}
+{cite:t}`kihlstrom_mirman1975` use the terms "reduced form" and "structural"
+models in a
+way that careful econometricians do.
+
+Reduced-form and structural models come in pairs.
+
+To each structure or structural model
+there is a reduced form, or collection of reduced forms, underlying different
+possible regressions.
+```
+
+The lecture is organized as follows.
+
+1. Set up the static two-commodity model and define equilibrium.
+2. State the price-revelation theorem and the invertibility conditions.
+3. Illustrate invertibility and its failure with numerical examples using CES
+ and
+ Cobb-Douglas preferences.
+4. Introduce the dynamic stochastic economy and derive the Bayesian convergence
+ result.
+5. Simulate Bayesian learning from price observations.
+
+This lecture builds on ideas in {doc}`blackwell_kihlstrom` and
+{doc}`likelihood_bayes`.
+
+We start by importing some Python packages.
+
+```{code-cell} ipython3
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import brentq
+from scipy.stats import norm
+```
+
+
+## Setup
+
+### Preferences, endowments, and the unknown state
+
+The economy has two goods.
+
+Good 2 is the numeraire (price normalized to 1).
+
+Good 1 trades at price $p > 0$.
+
+An unknown parameter $\bar{a}$ affects the value of good 1.
+
+Agent $i$'s expected utility
+from a bundle $(x_1^i, x_2^i)$ is
+
+$$
+U^i(x_1^i, x_2^i)
+ = \sum_{s=1}^{S} u^i(a_s x_1^i,\, x_2^i)\, P^i(\bar{a} = a_s),
+$$
+
+where $P^i$ is agent $i$'s subjective probability distribution over the finite
+state space
+$A = \{a_1, \ldots, a_S\}$.
+
+Each agent starts with an endowment $w^i$ of good 2 and a share $\theta^i$ of
+the
+representative firm.
+
+In the paper's formal model, a single firm transforms good 2 into good 1
+according to
+$y_1 = f(y_2)$ with $f' < 0$ and chooses production to maximize
+
+$$
+\pi(p) = \max_{y_2 \leq 0} \{p f(y_2) + y_2\}.
+$$
+
+The firm's profit $\pi$ is then distributed to households according to the
+shares
+$\theta^i$.
+
+Agent
+$i$'s budget constraint is
+
+$$
+p x_1^i + x_2^i = w^i + \theta^i \pi.
+$$
+
+Agents maximize expected utility subject to their budget constraints.
+
+A **competitive
+equilibrium** is a price $\hat{p}$ that clears both markets simultaneously.
+
+Under the maintained convexity assumptions equilibrium exists, and following
+{cite:t}`kihlstrom_mirman1975` we assume the equilibrium price is unique, so
+that we can write $\hat p = p(\mu)$ as a well-defined function of the informed
+agent's posterior.
+
+For most of what follows, the production side matters only through the induced
+equilibrium price map, so when we turn to numerical illustrations we will
+suppress production and use a pure-exchange / portfolio interpretation to keep
+the calculations transparent.
+
+### The informed agent's problem
+
+Suppose **agent 1** (the insider) observes a private signal $\tilde{y}$
+correlated with
+$\bar{a}$ before trading, where $\tilde{y}$ takes values in a finite set $Y$.
+
+Before the signal arrives, agent 1 has prior beliefs
+$\mu_0 = P^1$.
+
+Upon observing $\tilde{y} = y$, agent 1 updates to the
+**posterior** $\mu_y = (\mu_{y1}, \ldots, \mu_{yS})$ via Bayes' rule:
+
+$$
+\mu_{ys} = P(\bar{a} = a_s \mid \tilde{y} = y).
+$$
+
+Because agent 1's demand depends on $\mu_y$, the new equilibrium price satisfies
+
+$$
+\hat{p} = p(\mu_y).
+$$
+
+Outside observers who see $\hat{p}$ but not $\tilde{y}$ can try to *back out*
+the
+insider's posterior from the price.
+
+Define the set of realized posteriors
+
+$$
+M = \{\mu_y : y \in Y,\; P(\tilde y = y) > 0\}.
+$$
+
+The key question is whether the map $\mu \mapsto p(\mu)$ is one-to-one on $M$.
+
+To answer that question, we now translate "information in prices" into
+Blackwell's language of sufficiency.
+
+(price_revelation_theorem)=
+## Price revelation
+
+### Blackwell sufficiency
+
+The price variable $p(\mu_{\tilde{y}})$ *accurately transmits* the insider's
+private
+information if observing the equilibrium price is just as informative about
+$\bar{a}$ as
+observing the signal $\tilde{y}$ directly.
+
+In Blackwell's language ({cite:t}`blackwell1951` and {cite:t}`blackwell1953`),
+this means
+$p(\mu_{\tilde{y}})$ is **sufficient** for $\tilde{y}$.
+
+```{prf:definition} Sufficiency
+:label: ime_def_sufficiency
+
+A random variable $\tilde{y}$ is *sufficient* for $\tilde{y}'$ with
+respect to $\bar{a}$ if there exists a conditional distribution $P(y' \mid y)$,
+**independent of** $\bar{a}$, such that
+
+$$
+\phi'_a(y') = \sum_{y \in Y} P(y' \mid y)\, \phi_a(y)
+\quad \text{for all } a \text{ and all } y',
+$$
+
+where $\phi_a(y) = P(\tilde{y} = y \mid \bar{a} = a)$.
+
+Thus, once $\tilde{y}$ is known, $\tilde{y}'$ provides no additional information
+about $\bar{a}$.
+```
+
+{cite:t}`kihlstrom_mirman1975` show that
+
+```{prf:lemma} Posterior Sufficiency
+:label: ime_lemma_posterior_sufficiency
+
+The posterior distribution $\mu_{\tilde{y}}$ is a sufficient statistic for
+$\tilde{y}$.
+```
+
+```{prf:proof} (Sketch)
+The posterior $\mu_{\tilde{y}}$ satisfies
+
+$$
+P(\bar{a} = a_s \mid \mu_{\tilde{y}} = \mu_y,\; \tilde{y} = y) = \mu_{ys}
+ = P(\bar{a} = a_s \mid \mu_{\tilde{y}} = \mu_y).
+$$
+
+This identity says that once the posterior is known, conditioning on the
+original signal
+$\tilde y$ does not change beliefs about $\bar a$.
+
+Equivalently, the conditional law of $\tilde y$ given $\mu_{\tilde y}$ is
+independent of
+$\bar a$, so $\mu_{\tilde y}$ is sufficient for $\tilde y$ in Blackwell's sense.
+```
+
+Now let's think about the mapping from
+belief to price.
+
+```{prf:theorem} Price Revelation
+:label: ime_theorem_price_revelation
+
+In the model outlined above, the price random variable $p(\mu_{\tilde{y}})$ is
+sufficient for the random variable $\tilde{y}$ if and only if the function
+$p(P^1)$ is invertible on the set of prices
+
+$$
+\mathcal{P} = \Bigl\{\, p(\mu_y) : y \in Y,\;
+ P(\tilde{y} = y) = \sum_{a \in A} \phi_a(y)\,\mu_0(a) > 0 \Bigr\}.
+$$
+```
+
+The logic is
+
+$$
+\tilde y \quad \longrightarrow \quad \mu_{\tilde y} \quad \longrightarrow \quad
+p(\mu_{\tilde y}).
+$$
+
+The first arrow loses no information about $\bar a$ by
+{prf:ref}`ime_lemma_posterior_sufficiency`, and the theorem asks when the second
+arrow also loses no information.
+
+The proof has two parts.
+
+If $p(\cdot)$ is one-to-one on $M$, then observing the price is equivalent to
+observing the
+posterior itself because
+
+$$
+P(\mu_{\tilde y} = \mu \mid p(\mu_{\tilde y}) = p)
+= \begin{cases}
+1 & \text{if } \mu = p^{-1}(p), \\
+0 & \text{otherwise.}
+\end{cases}
+$$
+
+This conditional distribution is independent of the state, so price is
+sufficient for the
+posterior; together with {prf:ref}`ime_lemma_posterior_sufficiency`, price is
+therefore
+sufficient for the signal.
+
+Conversely, if two different posteriors in $M$ generated the same price, an
+observer of the price could not tell which posterior had occurred, and the paper
+shows formally that in this case the conditional distribution of the posterior
+given price would depend on the state, so price could not be sufficient.
+
+Before turning to invertibility itself, it helps to keep in mind the two
+economic interpretations emphasized in the paper.
+
+### Two interpretations
+
+#### Insider trading in a stock market
+
+Good 1 is a risky asset with random return $\bar{a}$; good 2 is "money".
+
+An insider's demand reveals private information about the return.
+
+If the invertibility condition holds, outside observers can read the insider's
+posterior distribution -- the useful information the insider's signal carries
+about $\bar a$ -- from the equilibrium stock price.
+
+#### Price as a quality signal
+
+Good 1 has uncertain quality $\bar{a}$.
+
+Experienced consumers (who have sampled the good) observe a signal correlated
+with quality
+and buy accordingly.
+
+Uninformed consumers can infer quality from the market price, provided
+invertibility holds.
+
+(invertibility_conditions)=
+## Invertibility and the elasticity of substitution
+
+When does the belief-to-price map fail to be invertible?
+
+{prf:ref}`ime_theorem_invertibility_conditions`
+shows that for a two-state economy ($S = 2$), the answer depends on the
+**elasticity of
+substitution** $\sigma$ of agent 1's utility function.
+
+Before stating the theorem, it helps to see the two intermediate steps in the
+paper's
+argument.
+
+```{prf:lemma} Same Price Implies Same Allocation
+:label: ime_lemma_same_price_same_allocation
+
+Assume that $u^i$ has continuous first partial derivatives and that $u^i$ is
+quasi-concave.
+
+Let $p \in \mathcal{P}$.
+
+If there exist two measures $\mu^*$ and $\mu'$ in $M$ such that
+$p(\mu^*, P^2, \ldots, P^n) = p(\mu', P^2, \ldots, P^n) = p$, then
+
+$$
+x^i(\mu^*, P^2, \ldots, P^n) = x^i(\mu', P^2, \ldots, P^n), \quad
+i = 1, \ldots, n.
+$$
+```
+
+Fix the beliefs of all agents except agent 1.
+
+The lemma says that if two posterior beliefs $\mu^*$ and $\mu'$ for agent 1
+both support the same equilibrium price $p$, then they support the same
+equilibrium allocation for every trader.
+
+The intuition is that when the price is unchanged, the demands of the
+uninformed traders are unchanged too, so market clearing forces the informed
+agent's bundle to be unchanged as well.
+
+This lemma lets us define the informed agent's equilibrium bundle as a function
+of price alone:
+
+$$
+x(p) = (x_1(p), x_2(p)).
+$$
+
+Throughout, $u^i_j$ denotes the partial derivative of $u^i$ with respect to its
+$j$-th argument.
+
+Whenever the informed agent consumes positive amounts of both goods, optimality
+of $x(p)$
+under posterior $\mu$ gives the interior first-order condition
+
+$$
+p = \frac{\sum_{s=1}^S a_s u_1^1(a_s x_1(p), x_2(p))\, \mu(a_s)}
+ {\sum_{s=1}^S u_2^1(a_s x_1(p), x_2(p))\, \mu(a_s)}.
+$$
+
+For a fixed price $p$, the bundle $x(p)$ is fixed too, so invertibility boils
+down to
+whether this equation admits a unique posterior $\mu$.
+
+```{prf:lemma} Unique Posterior at a Given Price
+:label: ime_lemma_unique_posterior
+
+Assume that the first partial derivatives of $u^1$ exist and that $u^1$ is
+quasi-concave.
+
+Also assume that agent 1 always consumes positive quantities of both goods.
+
+Then $p(P^1)$ is invertible on $\mathcal{P}$ if for each $p \in \mathcal{P}$
+there exists a unique probability measure $\mu \in M$ such that
+
+$$
+\frac{\sum_{s=1}^S a_s\, u^1_1(a_s x_1(p), x_2(p))\, \mu(a_s)}
+ {\sum_{s=1}^S u^1_2(a_s x_1(p), x_2(p))\, \mu(a_s)} = p.
+$$
+```
+
+If two different posteriors gave the same price, then by
+{prf:ref}`ime_lemma_same_price_same_allocation` they would share the same bundle
+$x(p)$, contradicting uniqueness of the posterior that solves the first-order
+condition at that price.
+
+### The two-state first-order condition
+
+With $S = 2$ and $\mu = (q,\, 1-q)$, define
+
+$$
+\alpha_s(p) = a_s\, u^1_1(a_s x_1(p),\, x_2(p)), \qquad
+\beta_s(p) = u^1_2(a_s x_1(p),\, x_2(p)), \qquad s = 1, 2.
+$$
+
+Then the first-order condition becomes
+
+$$
+p = \frac{\alpha_1(p)\, q + \alpha_2(p)\, (1-q)}
+ {\beta_1(p)\, q + \beta_2(p)\, (1-q)}.
+$$
+
+At a fixed price $p$, the quantities $\alpha_s(p)$ and $\beta_s(p)$ are
+constants, so
+uniqueness of the posterior is the same as uniqueness of the scalar $q$ solving
+this
+equation.
+
+```{prf:theorem} Invertibility Conditions
+:label: ime_theorem_invertibility_conditions
+
+Assume that the first partial derivatives of $u^1$ exist and that $u^1$ is
+quasi-concave and homothetic.
+
+Also suppose that the informed agent always consumes positive quantities of
+both goods in all equilibrium allocations.
+
+If $S = 2$ and the elasticity of substitution of $u^1$ is either always less
+than one or always greater than one, then $p(P^1)$ is invertible on
+$\mathcal{P}$.
+
+If $u^1$ is Cobb-Douglas (elasticity of substitution constant and equal to
+one), then $p(P^1)$ is constant on $\mathcal{P}$.
+```
+
+When $\sigma = 1$ the income and substitution effects exactly cancel, so
+agent 1's demand for good 1 does not respond to changes in beliefs about
+$\bar{a}$.
+
+Because the demand is unchanged, the market-clearing price is unchanged too,
+and the price reveals nothing about the insider's signal.
+
+### CES utility
+
+For concreteness we work with a simplified example with the **constant-elasticity-of-substitution** (CES)
+utility
+function
+
+$$
+u(c_1, c_2) = \bigl(c_1^{\rho} + c_2^{\rho}\bigr)^{1/\rho}, \qquad \rho \in
+(-\infty,0) \cup (0,1),
+$$
+
+whose elasticity of substitution is $\sigma = 1/(1-\rho)$.
+
+- $\rho \to 0$: Cobb-Douglas ($\sigma = 1$).
+- $\rho < 0$: $\sigma < 1$ (complements).
+- $0 < \rho < 1$: $\sigma > 1$ (substitutes).
+
+Pertinent partial derivatives are
+
+$$
+u_1(c_1,c_2) = \bigl(c_1^\rho + c_2^\rho\bigr)^{1/\rho - 1}\, c_1^{\rho-1},
+\qquad
+u_2(c_1,c_2) = \bigl(c_1^\rho + c_2^\rho\bigr)^{1/\rho - 1}\, c_2^{\rho-1}.
+$$
+
+This CES example is only an illustration, because the theorem itself covers any
+homothetic utility with elasticity everywhere above one or everywhere below one.
+
+With that example in hand, we can compute the equilibrium price directly as a
+function of the posterior.
+
+### Equilibrium price as a function of the posterior
+
+We focus on agent 1 as the *only* informed trader who absorbs one unit of good 1
+at
+equilibrium (i.e., $x_1 = 1$).
+
+Let $W_1 = w^1 + \theta^1 \pi$ denote agent 1's total wealth (endowment plus
+profit share).
+
+Agent 1's budget constraint then reduces to
+$x_2 = W_1 - p$, and the equilibrium price is the unique $p \in (0, W_1)$
+satisfying
+the first-order condition
+
+$$
+p \bigl[q\, u_2(a_1,\, W_1-p) + (1-q)\, u_2(a_2,\, W_1-p)\bigr]
+= q\, a_1\, u_1(a_1,\, W_1-p) + (1-q)\, a_2\, u_1(a_2,\, W_1-p).
+$$
+
+For Cobb-Douglas utility ($\sigma = 1$), the first-order condition becomes $p =
+W_1 - p$,
+giving $p^* = W_1/2$ regardless of the posterior $q$, confirming that no
+information
+is transmitted through the price in the Cobb-Douglas case.
+
+We compute first-order conditions numerically below.
+
+```{code-cell} ipython3
+def ces_derivatives(c1, c2, ρ):
+ """
+ Return CES marginal utilities.
+
+ Use the Cobb-Douglas limit near rho = 0.
+ """
+ if abs(ρ) < 1e-4:
+ u1 = 0.5 * np.sqrt(c2 / c1)
+ u2 = 0.5 * np.sqrt(c1 / c2)
+ else:
+ common = (c1**ρ + c2**ρ)**(1 / ρ - 1)
+ u1 = common * c1**(ρ - 1)
+ u2 = common * c2**(ρ - 1)
+ return u1, u2
+
+
+def eq_price(q, a1, a2, W1, ρ):
+ """Return the equilibrium price for posterior q."""
+ def residual(p):
+ x2 = W1 - p
+ u1_s1, u2_s1 = ces_derivatives(a1, x2, ρ)
+ u1_s2, u2_s2 = ces_derivatives(a2, x2, ρ)
+ lhs = p * (q * u2_s1 + (1 - q) * u2_s2)
+ rhs = q * a1 * u1_s1 + (1 - q) * a2 * u1_s2
+ return lhs - rhs
+
+ try:
+ return brentq(residual, 1e-6, W1 - 1e-6, xtol=1e-10)
+ except ValueError:
+ return np.nan
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: equilibrium price vs posterior
+ name: fig-eq-price-posterior
+---
+a1, a2 = 2.0, 0.5 # state values (a1 > a2)
+W1 = 4.0
+
+q_grid = np.linspace(0.05, 0.95, 200)
+
+ρ_values = [-0.5, 0.0, 0.5]
+ρ_labels = [
+ r"$\rho = -0.5$ ($\sigma = 0.67$, complements)",
+ r"$\rho = 0$ ($\sigma = 1$, Cobb-Douglas)",
+ r"$\rho = 0.5$ ($\sigma = 2$, substitutes)",
+]
+
+fig, ax = plt.subplots(figsize=(8, 5))
+
+for ρ, label in zip(ρ_values, ρ_labels):
+ prices = [eq_price(q, a1, a2, W1, ρ) for q in q_grid]
+ ax.plot(q_grid, prices, label=label, lw=2)
+
+ax.set_xlabel(r"posterior probability $q = \Pr(\bar{a} = a_1)$", fontsize=12)
+ax.set_ylabel("equilibrium price $p^*(q)$", fontsize=12)
+ax.legend(fontsize=10)
+plt.tight_layout()
+plt.show()
+```
+
+The plot confirms {prf:ref}`ime_theorem_invertibility_conditions`.
+
+For CES with $\sigma \neq 1$, the equilibrium price is strictly monotone in $q$.
+
+An outside observer who knows the equilibrium map $p^*(\cdot)$ can therefore
+invert the price uniquely to recover $q$, so the inside information is fully
+transmitted.
+
+For Cobb-Douglas ($\sigma = 1$), the price is flat in $q$, so information is
+never transmitted through the market.
+
+```{code-cell} ipython3
+p_cd = [eq_price(q, a1, a2, W1, ρ=0.0) for q in q_grid]
+
+print(f"Cobb-Douglas (rho=0): min p* = {min(p_cd):.6f}, "
+ f"max p* = {max(p_cd):.6f}, "
+ f"range = {max(p_cd)-min(p_cd):.2e}")
+print(f"Analytical CD price = W1/2 = {W1/2:.6f}")
+```
+
+Every entry equals $W_1/2 = 2.0$ exactly, confirming analytically that the
+Cobb-Douglas
+equilibrium price is independent of $q$ and of the state values $a_1, a_2$.
+
+The numerical plot shows monotonicity, and the next subsection connects that
+pattern back to the proof of {prf:ref}`ime_theorem_invertibility_conditions`.
+
+(price_monotonicity)=
+### Why monotonicity depends on $\sigma$
+
+Fix a price $p$ and treat $\alpha_s(p)$ and $\beta_s(p)$ as constants.
+
+The right-hand side of the two-state first-order condition
+
+$$
+\frac{\alpha_1(p)\, q + \alpha_2(p)\, (1-q)}
+ {\beta_1(p)\, q + \beta_2(p)\, (1-q)}
+$$
+
+is then a function of $q$ alone, with derivative
+
+$$
+\frac{\partial}{\partial q}
+\frac{\alpha_1 q + \alpha_2 (1-q)}
+ {\beta_1 q + \beta_2 (1-q)}
+= \frac{\alpha_1 \beta_2 - \alpha_2 \beta_1}
+ {\bigl[\beta_1 q + \beta_2 (1-q)\bigr]^2}.
+$$
+
+So the sign is determined by $\alpha_1 \beta_2 - \alpha_2 \beta_1$, and if that
+sign is constant then for each fixed price there is at most one posterior weight
+$q$ consistent with the first-order condition, which is exactly what
+{prf:ref}`ime_theorem_invertibility_conditions` requires.
+
+Using
+
+$$
+\frac{\alpha_s}{\beta_s}
+ = \frac{a_s\, u_1(a_s x_1, x_2)}{u_2(a_s x_1, x_2)}
+ = a_s^{(\sigma-1)/\sigma}\,\Bigl(\frac{x_2}{x_1}\Bigr)^{1/\sigma},
+$$
+
+one can show
+
+$$
+\frac{\partial}{\partial a}\,\frac{\alpha}{\beta}
+ = \frac{(\sigma - 1)}{\sigma}\, a^{-1/\sigma}\,
+ \Bigl(\frac{x_2}{x_1}\Bigr)^{1/\sigma}.
+$$
+
+For the CES specification, this derivative is positive when $\sigma > 1$,
+negative when
+$\sigma < 1$, and *zero when $\sigma = 1$*.
+
+In other words, for CES utility the ratio $\alpha_s / \beta_s$ moves
+monotonically with the state value $a_s$ unless $\sigma = 1$, which makes the
+fixed-price first-order-condition expression monotone in $q$ and in turn
+delivers invertibility.
+
+The vanishing derivative in the Cobb-Douglas case means the marginal rate of
+substitution is
+independent of $a_s$, so the informed agent's demand, and hence the equilibrium
+price, does
+not respond to changes in beliefs.
+
+Let us visualize the ratio $\alpha_s / \beta_s$ as a function of $a_s$ for
+different
+values of $\sigma$:
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: marginal rate of substitution
+ name: fig-mrs-alpha-beta
+---
+a_vals = np.linspace(0.3, 3.0, 300)
+x1_fix, x2_fix = 1.0, 1.0
+
+fig, ax = plt.subplots(figsize=(7, 4))
+for ρ in [-0.5, -1e-6, 0.5]:
+ σ = 1 / (1 - ρ) if abs(ρ) > 1e-8 else 1.0
+ ratios = []
+ for a in a_vals:
+ u1, u2 = ces_derivatives(a * x1_fix, x2_fix, ρ)
+ ratios.append(a * u1 / u2)
+ ax.plot(a_vals, ratios, label=rf"$\sigma = {σ:.2f}$", lw=2)
+
+ax.set_xlabel(r"state value $a_s$", fontsize=12)
+ax.set_ylabel(r"$\alpha_s / \beta_s = a_s u_1 / u_2$", fontsize=12)
+ax.axhline(y=1.0, color="black", lw=0.8, ls="--")
+ax.legend(fontsize=10)
+plt.tight_layout()
+plt.show()
+```
+
+When $\sigma = 1$ the ratio is constant across all $a_s$ values, so
+information about the state has no effect on the marginal rate of substitution.
+
+For $\sigma < 1$ the ratio is decreasing in $a_s$, and for $\sigma > 1$ it is
+increasing, making the equilibrium price strictly monotone in the posterior $q$
+in both cases.
+
+The static analysis asks whether a current price reveals current private
+information, whereas the next section asks what a whole history of prices
+reveals over time.
+
+(bayesian_price_expectations)=
+## Bayesian price expectations in a dynamic economy
+
+We now turn to a question addressed in Section 3 of
+{cite:t}`kihlstrom_mirman1975`.
+
+### A stochastic exchange economy
+
+Time is discrete: $t = 1, 2, \ldots$
+
+In each period $t$:
+
+1. Consumer $i$ receives a random endowment $\omega_i^t$.
+2. Markets open; competitive prices $p^t = p(\omega^t)$ clear all markets.
+3. Consumers trade and consume.
+
+The endowment vectors $\{\tilde{\omega}^t\}$ are **i.i.d.** with density
+$f(\omega^t \mid \lambda)$, where $\lambda = (\lambda_1, \ldots, \lambda_K)$ is
+a **structural parameter vector** (of dimension $K$) that is *fixed but
+unknown*.
+
+The equilibrium price at time $t$ is a deterministic function of $\omega^t$, so
+$\{p^t\}$ is also i.i.d.
+
+For any measurable price set $P$, let
+
+$$
+W(P) = \{\omega^t : p(\omega^t) \in P\}.
+$$
+
+Then
+
+$$
+P_\lambda(p^t \in P) = P_\lambda(\omega^t \in W(P))
+= \int_{W(P)} f(\omega^t \mid \lambda)\, d\omega^t.
+$$
+
+The induced price density is denoted by $g(p^t \mid \lambda)$.
+
+For a given structure $\lambda$, this density is the observable implication of
+the model, and when several structures imply the same density we group them
+into a single reduced-form class.
+
+The next issue is therefore what an observer can and cannot infer about the
+structure from price data alone.
+
+### The identification problem
+
+Because price observations identify only the induced price density
+$g(\cdot \mid \lambda)$, and because the structural-to-reduced-form map
+$\lambda \mapsto g(\cdot \mid \lambda)$ may be many-to-one, price data may
+identify only a reduced-form class rather than the exact structure.
+
+In particular, it may be impossible to recover $\lambda$ from
+$g(p \mid \lambda)$ even with infinite price data.
+
+To handle this, partition $\Lambda$ into equivalence classes $\mu$ such that
+$\lambda \in \mu$ and $\lambda' \in \mu$ whenever $g(p \mid \lambda) = g(p \mid
+\lambda')$
+for all $p$.
+
+The equivalence class $\mu$ containing the true $\lambda$ is the **reduced
+form** relevant for price data.
+
+An observer who knows the infinite price history learns
+$\mu$ but not necessarily $\lambda$.
+
+Once that distinction is clear, Bayesian updating can be written down directly.
+
+### Bayesian updating
+
+An uninformed observer begins with a prior $h(\lambda)$ over $\lambda \in
+\Lambda$.
+
+If the observer could see endowments directly, the posterior would be
+
+$$
+h(\lambda \mid \omega^1, \ldots, \omega^t)
+ = \frac{h(\lambda)\, \prod_{\tau=1}^{t} f(\omega^\tau \mid \lambda)}
+ {\displaystyle\sum_{\lambda' \in \Lambda}
+ h(\lambda')\, \prod_{\tau=1}^{t} f(\omega^\tau \mid \lambda')},
+$$
+
+and the paper appeals to a Bayesian consistency result to conclude that this
+posterior concentrates on the true structure $\bar \lambda$.
+
+After observing the price sequence $(p^1, \ldots, p^t)$, the observer's Bayesian
+posterior is
+
+$$
+h(\lambda \mid p^1, \ldots, p^t)
+ = \frac{h(\lambda)\, \prod_{\tau=1}^{t} g(p^\tau \mid \lambda)}
+ {\displaystyle\sum_{\lambda' \in \Lambda}
+ h(\lambda')\, \prod_{\tau=1}^{t} g(p^\tau \mid \lambda')}.
+$$
+
+Price data cannot distinguish structures inside the same reduced-form class.
+
+Indeed, if
+$\lambda$ and $\lambda'$ belong to the same class $\mu$, then
+$g(\cdot \mid \lambda) = g(\cdot \mid \lambda')$, so
+
+$$
+\frac{h(\lambda \mid p^1, \ldots, p^t)}
+ {h(\lambda' \mid p^1, \ldots, p^t)}
+= \frac{h(\lambda)}{h(\lambda')}
+$$
+
+for every sample history, so the relative odds within an observationally
+equivalent class never change.
+
+At time $t$, the observer's price expectations for the next period are
+
+$$
+g(p^{t+1} \mid p^1, \ldots, p^t)
+ = \sum_{\lambda \in \Lambda} g(p^{t+1} \mid \lambda)\,
+ h(\lambda \mid p^1, \ldots, p^t).
+$$
+
+### The convergence theorem
+
+```{prf:theorem} Bayesian Convergence
+:label: ime_theorem_bayesian_convergence
+
+Let $\bar\lambda$ be the true
+structural parameter and $\bar\mu$ the reduced form that contains $\bar\lambda$.
+
+Assume the prior assigns positive probability to the reduced-form class $\bar\mu$.
+
+Define the posterior mass on a reduced-form class by
+
+$$
+H_t(\mu) = \sum_{\lambda \in \mu} h(\lambda \mid p^1, \ldots, p^t).
+$$
+
+Because all structures inside a class imply the same $g(\cdot \mid \lambda)$,
+the
+predictive density can equivalently be written as
+
+$$
+g(p^{t+1} \mid p^1, \ldots, p^t)
+ = \sum_{\mu} g(p^{t+1} \mid \mu)\, H_t(\mu).
+$$
+
+Then
+
+$$
+\lim_{t \to \infty} H_t(\mu)
+ = \begin{cases} 1 & \text{if } \mu = \bar\mu, \\ 0 & \text{otherwise,}
+ \end{cases}
+$$
+
+with probability one.
+
+Consequently,
+
+$$
+\lim_{t \to \infty} g(p^{t+1} \mid p^1, \ldots, p^t) = g(p \mid \bar\mu),
+$$
+
+which equals the rational-expectations price distribution for a fully informed
+observer.
+```
+
+```{note}
+Note that the theorem only requires the prior to assign positive probability to the reduced-form class $\bar\mu$ that contains the true structure $\bar\lambda$.
+
+This is implied by, but weaker than, assigning positive probability to the true
+structural parameter $\bar\lambda$ itself.
+
+A prior could place zero mass on $\bar\lambda$
+while still placing positive mass on other structures inside $\bar\mu$.
+```
+
+The important distinction is that price observers need not learn $\bar \lambda$
+itself.
+
+They only learn which reduced-form class is correct.
+
+That is enough for forecasting because every $\lambda \in \bar \mu$ generates
+the same price density $g(\cdot \mid \bar \mu)$.
+
+Rational price expectations emerge from
+learning the
+reduced form, not from identifying every structural detail of the economy.
+
+Here "rational expectations" means that the observer's predictive distribution
+for next
+period's price matches the objective price distribution generated by the true
+reduced form.
+
+Let's now turn to a simple simulation.
+
+(bayesian_simulation)=
+## Simulating Bayesian learning from prices
+
+We illustrate the theorem with a two-state example.
+
+Two possible reduced forms $\mu_1$ and $\mu_2$ generate prices
+$p^t \sim N(\bar{p}_i, \sigma_p^2)$ for $i = 1, 2$ respectively.
+
+The observer knows the two possible price distributions (the reduced forms) but
+not which
+one governs the data.
+
+This is a **Bayesian model selection** problem we have seen in {doc}`likelihood_bayes`.
+
+With a prior $h_0$ on $\mu_1$ and the observed price $p^t$, the posterior weight
+on $\mu_1$
+after period $t$ is
+
+$$
+h_t = \frac{h_{t-1}\, g(p^t \mid \mu_1)}{h_{t-1}\, g(p^t \mid \mu_1)
+ + (1-h_{t-1})\, g(p^t \mid \mu_2)}.
+$$
+
+We consider a numerical example with two normal distributions with different means
+
+```{code-cell} ipython3
+def simulate_bayesian_learning(
+ p_bar_true, p_bar_alt, σ_p, T, h0, n_paths, seed=42
+):
+ """Simulate posterior learning between two Gaussian reduced forms."""
+ rng = np.random.default_rng(seed)
+ h_paths = np.zeros((n_paths, T + 1))
+ h_paths[:, 0] = h0
+
+ for path in range(n_paths):
+ h = h0
+ prices = rng.normal(p_bar_true, σ_p, size=T)
+ for t, p in enumerate(prices):
+ g_true = norm.pdf(p, loc=p_bar_true, scale=σ_p)
+ g_alt = norm.pdf(p, loc=p_bar_alt, scale=σ_p)
+ denom = h * g_true + (1 - h) * g_alt
+ h = h * g_true / denom
+ h_paths[path, t + 1] = h
+
+ return h_paths
+
+
+def plot_bayesian_learning(h_paths, p_bar_true, p_bar_alt, ax):
+ """Plot posterior beliefs over time."""
+ T = h_paths.shape[1] - 1
+ t_grid = np.arange(T + 1)
+
+ for path in h_paths:
+ ax.plot(t_grid, path, alpha=0.25, lw=0.8, color="steelblue")
+
+ median_path = np.median(h_paths, axis=0)
+ ax.plot(t_grid, median_path, color="navy", lw=2, label="median posterior")
+
+ ax.axhline(
+ y=1.0,
+ color="black",
+ ls="--",
+ lw=1.2,
+ label="true model weight = 1",
+ )
+ ax.set_xlabel("period $t$", fontsize=12)
+ ax.set_ylabel(r"$h_t$ = posterior weight on true model", fontsize=12)
+ ax.legend(fontsize=10)
+```
+
+We consider two cases, one that is easy to learn and another one that is harder to learn,
+using $T = 300$ periods, $n = 40$ simulated paths, a diffuse prior $h_0 = 0.5$, and
+common standard deviation $\sigma_p = 0.4$.
+
+- *Easy case*: true model $N(2.0,\, 0.4^2)$, alternative $N(1.2,\, 0.4^2)$.
+- *Hard case*: true model $N(2.0,\, 0.4^2)$, alternative $N(1.8,\, 0.4^2)$.
+
+Whether easy or hard to learn depends on "how close" the true distribution is compared to the
+alternative hypothesis.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: bayesian learning across paths
+ name: fig-bayesian-learning
+---
+T = 300
+h0 = 0.5 # diffuse prior
+n_paths = 40
+σ_p = 0.4
+
+fig, axes = plt.subplots(1, 2, figsize=(12, 5))
+
+# Distinct reduced forms
+p_bar_true, p_bar_alt = 2.0, 1.2
+h_paths = simulate_bayesian_learning(p_bar_true, p_bar_alt, σ_p, T, h0, n_paths)
+plot_bayesian_learning(h_paths, p_bar_true, p_bar_alt, axes[0])
+
+# Similar reduced forms
+p_bar_true, p_bar_alt = 2.0, 1.8
+h_paths_hard = simulate_bayesian_learning(
+ p_bar_true, p_bar_alt, σ_p, T, h0, n_paths
+)
+plot_bayesian_learning(h_paths_hard, p_bar_true, p_bar_alt, axes[1])
+
+plt.tight_layout()
+plt.show()
+```
+
+In both panels the posterior weight on the true model converges to 1 with
+probability one,
+though convergence is slower when the two price distributions are similar (right
+panel).
+
+### Price expectations vs. rational expectations
+
+We now verify that the observer's price expectations converge to the
+rational-expectations
+distribution $g(p \mid \bar\mu)$.
+
+We use the parameterization of the "hard-to-learn" example above
+($\bar{p}_{\text{true}} = 2.0$, $\bar{p}_{\text{alt}} = 1.8$, $\sigma_p = 0.4$),
+extending to $T = 1{,}000$ periods with a single simulated path and prior $h_0 = 0.5$
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: price distribution convergence
+ name: fig-price-convergence
+---
+def price_expectation(h_t, p_bar_true, p_bar_alt, σ_p, p_grid):
+ """Return the predictive price density at posterior weight h_t."""
+ return (
+ h_t * norm.pdf(p_grid, loc=p_bar_true, scale=σ_p)
+ + (1 - h_t) * norm.pdf(p_grid, loc=p_bar_alt, scale=σ_p)
+ )
+
+
+p_bar_true, p_bar_alt = 2.0, 1.8
+σ_p = 0.4
+n_paths = 1
+T_long = 1000
+
+h_paths_long = simulate_bayesian_learning(
+ p_bar_true, p_bar_alt, σ_p, T_long, h0=0.5, n_paths=n_paths, seed=7
+)
+
+p_grid = np.linspace(0.0, 3.5, 300)
+re_density = norm.pdf(p_grid, loc=p_bar_true, scale=σ_p)
+
+fig, ax = plt.subplots(figsize=(8, 5))
+snapshots = [0, 25, 100, 300, 1000]
+palette = plt.cm.Blues(np.linspace(0.3, 1.0, len(snapshots)))
+
+for t_snap, col in zip(snapshots, palette):
+ h_t = h_paths_long[0, t_snap]
+ dens = price_expectation(h_t, p_bar_true, p_bar_alt, σ_p, p_grid)
+ ax.plot(
+ p_grid,
+ dens,
+ color=col,
+ lw=2,
+ label=rf"$t = {t_snap}$, $h_t = {h_t:.3f}$",
+ )
+
+ax.plot(p_grid, re_density, "k--", lw=2,
+ label=r"rational expectations $g(p \mid \bar{\mu})$")
+ax.set_xlabel("price $p$", fontsize=12)
+ax.set_ylabel("density", fontsize=12)
+ax.legend(fontsize=9)
+plt.tight_layout()
+plt.show()
+```
+
+The sequence of predictive densities (shades of blue) converges to the
+rational-expectations
+density (dashed black line) as experience accumulates.
+
+This illustrates {prf:ref}`ime_theorem_bayesian_convergence`.
+
+We can now sharpen the point by looking at a case in which the reduced form is
+learned but the underlying structure is not.
+
+(km_extension_nonidentification)=
+### Learning the reduced form without identifying the structure
+
+The convergence result is particularly striking because the observer converges
+to
+*rational expectations* even when the underlying **structure** $\lambda$ is
+*not identified* by prices.
+
+To illustrate this, consider a case with *three* possible structures
+$\lambda^{(1)}, \lambda^{(2)}, \lambda^{(3)}$ but only *two* reduced forms
+$\mu_1 = \{\lambda^{(1)}, \lambda^{(2)}\}$ and $\mu_2 = \{\lambda^{(3)}\}$
+(because $\lambda^{(1)}$ and $\lambda^{(2)}$ generate the same price
+distribution).
+
+We continue with the hard-to-learn parameterization, so the three structures
+have price means $\bar{p}_1 = \bar{p}_2 = 2.0$ and $\bar{p}_3 = 1.8$, with
+common standard deviation $\sigma_p = 0.4$, a uniform prior
+$h_0 = (1/3, 1/3, 1/3)$, and $T = 400$ periods over $30$ paths.
+
+The true structure is $\lambda^{(1)}$.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: learning with non-identification
+ name: fig-nonidentification
+---
+def simulate_learning_3struct(
+ T, h0_vec, p_bar_vec, σ_p, true_idx, n_paths, seed=0
+):
+ """Simulate learning with three structures and two reduced forms."""
+ rng = np.random.default_rng(seed)
+ h_paths = np.zeros((n_paths, T + 1, 3))
+ h_paths[:, 0, :] = h0_vec
+
+ for path in range(n_paths):
+ h = np.array(h0_vec, dtype=float)
+ prices = rng.normal(p_bar_vec[true_idx], σ_p, size=T)
+ for t, p in enumerate(prices):
+ likelihoods = norm.pdf(p, loc=p_bar_vec, scale=σ_p)
+ h = h * likelihoods
+ h /= h.sum()
+ h_paths[path, t + 1, :] = h
+
+ return h_paths
+
+
+# Structures 0 and 1 share the same reduced form
+p_bar_vec = np.array([2.0, 2.0, 1.8])
+h0_vec = np.array([1 / 3, 1 / 3, 1 / 3])
+σ_p = 0.4
+T = 400
+true_idx = 0 # Structure 0 is observationally equivalent to 1
+
+h_paths_3 = simulate_learning_3struct(
+ T, h0_vec, p_bar_vec, σ_p, true_idx, n_paths=30
+)
+t_grid = np.arange(T + 1)
+
+fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True)
+struct_labels = [
+ r"$\lambda^{(1)}$",
+ r"$\lambda^{(2)}$",
+ r"$\lambda^{(3)}$",
+]
+
+for k, (ax, label) in enumerate(zip(axes, struct_labels)):
+ for path in h_paths_3:
+ ax.plot(t_grid, path[:, k], alpha=0.25, lw=0.8, color="steelblue")
+ ax.plot(t_grid, np.median(h_paths_3[:, :, k], axis=0),
+ color="navy", lw=2, label=f"median weight on {label}")
+ ax.set_xlabel("period $t$", fontsize=11)
+ ax.legend(fontsize=9)
+
+axes[0].set_ylabel("posterior weight", fontsize=11)
+plt.tight_layout()
+plt.show()
+```
+
+The observer correctly rules out $\lambda^{(3)}$ (the wrong reduced form) with
+probability
+one, but cannot distinguish $\lambda^{(1)}$ from $\lambda^{(2)}$ because they
+generate an
+identical price distribution.
+
+Nevertheless, the observer's **price expectations** converge
+to rational expectations because both structures imply the same reduced form
+$\bar\mu$.
+
+
+## Exercises
+
+```{exercise}
+:label: km_ex1
+
+**CARA portfolio utility and the stock-market interpretation.**
+
+Consider a two-state economy ($a_1 = 2$, $a_2 = 0.5$) where the informed agent has
+**CARA** (constant absolute risk aversion) preferences over portfolio wealth:
+
+$$
+u(W) = -e^{-\gamma W}, \quad W = x_2 + \bar{a}\, x_1.
+$$
+
+The agent chooses $x_1$ to maximize
+
+$$
+q\,u(W_1) + (1-q)\,u(W_2), \quad W_s = w - p\,x_1 + a_s\,x_1,
+$$
+
+subject to the budget constraint $p\,x_1 + x_2 = w$.
+
+Total supply of good 1 is $X_1 = 1$.
+
+1. Derive the first-order condition for the informed agent's optimal $x_1$.
+
+1. Use the market-clearing condition $x_1 = 1$ (the informed agent absorbs the
+ entire supply) to obtain an implicit equation for the equilibrium price
+ $p^*(q)$, and solve it numerically for $q \in (0,1)$ and several values of
+ $\gamma$.
+
+1. Show *analytically* that $p^*(q)$ admits the closed form
+
+ $$
+ p^*(q) = \frac{a_2 + R(q,\gamma)\, a_1}{1 + R(q,\gamma)},
+ \qquad R(q,\gamma) = \frac{q}{1-q}\, e^{-\gamma(a_1-a_2)},
+ $$
+
+ and verify that $p^*(q)$ is strictly increasing in $q$.
+```
+
+```{solution-start} km_ex1
+:class: dropdown
+```
+
+For the first-order condition, define $W_s = w + (a_s - p)\,x_1$ for
+$s = 1, 2$.
+
+Then the FOC is
+
+$$
+q\,(a_1 - p)\,\gamma\, e^{-\gamma W_1}
+= (1-q)\,(p - a_2)\,\gamma\, e^{-\gamma W_2},
+$$
+
+or equivalently (dividing by $\gamma$ and rearranging)
+
+$$
+q\,(a_1 - p)\, e^{-\gamma(a_1-p) x_1}
+ = (1-q)\,(p - a_2)\, e^{\gamma(p-a_2) x_1}.
+$$
+
+Setting $x_1 = 1$ (the informed agent absorbs all supply), this becomes a
+scalar root-finding problem in $p$:
+
+$$
+F(p;\,q,\gamma) \equiv
+ q\,(a_1-p)\,e^{-\gamma(a_1-p)} - (1-q)\,(p-a_2)\,e^{\gamma(p-a_2)} = 0.
+$$
+
+```{code-cell} ipython3
+from scipy.optimize import brentq
+
+def F_cara(p, q, a1, a2, γ, x1=1.0):
+ """Residual for the CARA equilibrium condition."""
+ return (q * (a1 - p) * np.exp(-γ * (a1 - p) * x1)
+ - (1 - q) * (p - a2) * np.exp(γ * (p - a2) * x1))
+
+a1, a2 = 2.0, 0.5
+q_grid = np.linspace(0.05, 0.95, 200)
+γ_values = [0.5, 1.0, 2.0, 5.0]
+colors_sol = plt.cm.plasma(np.linspace(0.15, 0.85, len(γ_values)))
+
+fig, ax = plt.subplots(figsize=(8, 5))
+for γ, color in zip(γ_values, colors_sol):
+ p_eq = [brentq(F_cara, a2, a1,
+ args=(q, a1, a2, γ))
+ for q in q_grid]
+ ax.plot(q_grid, p_eq, lw=2, color=color,
+ label=rf"$\gamma = {γ}$")
+
+ax.set_xlabel(r"posterior $q = \Pr(\bar a = a_1)$", fontsize=12)
+ax.set_ylabel("equilibrium price $p^*(q)$", fontsize=12)
+ax.set_title("CARA preferences: equilibrium prices", fontsize=12)
+ax.legend(fontsize=10)
+plt.tight_layout()
+plt.show()
+```
+
+The price is strictly increasing in $q$ for every $\gamma > 0$.
+
+For the closed form, start from the FOC at $x_1 = 1$, divide both sides by
+$(a_1 - p)(p - a_2)$, and combine the exponentials:
+
+$$
+\frac{q\,(a_1 - p)}{(1-q)\,(p - a_2)} = e^{\gamma(a_1 - a_2)}.
+$$
+
+Rearranging gives
+
+$$
+\frac{p - a_2}{a_1 - p} = \frac{q}{1-q}\, e^{-\gamma(a_1 - a_2)}
+\equiv R(q,\gamma),
+$$
+
+and solving the resulting linear equation in $p$ yields
+
+$$
+p^*(q) = \frac{a_2 + R(q,\gamma)\, a_1}{1 + R(q,\gamma)}.
+$$
+
+Since $R(q,\gamma)$ is strictly increasing in $q$ and
+$dp^*/dR = (a_1 - a_2)/(1 + R)^2 > 0$, the equilibrium price $p^*(q)$ is
+strictly increasing in $q$.
+
+This exercise uses the stock-market interpretation emphasized by
+{cite:t}`kihlstrom_mirman1975`.
+
+Portfolio wealth is $W = x_2 + \bar{a}\, x_1$, so $a x_1$ and $x_2$ are perfect
+substitutes in each state.
+
+Hence the elasticity of substitution between the two arguments of
+$u(a x_1, x_2)$ is infinite, corresponding to the $\sigma > 1$ side of
+{prf:ref}`ime_theorem_invertibility_conditions`.
+
+The difference is that this example is not the full equilibrium model the theorem analyzes, but rather a partial equilibrium model with a single informed agent and a fixed supply of the risky asset.
+
+```{solution-end}
+```
+
+```{exercise}
+:label: km_ex2
+
+In the Bayesian learning simulation, the speed of
+convergence to rational expectations is determined by the **Kullback-Leibler
+divergence**
+between the two reduced forms.
+
+The KL divergence $D_{KL}(\mu_1 \| \mu_2)$ from $g(\cdot \mid \mu_1)$ to
+$g(\cdot \mid \mu_2)$, for two normal distributions with means $\bar{p}_1$ and
+$\bar{p}_2$ and common variance $\sigma_p^2$, is
+
+$$
+D_{KL}(\mu_1 \| \mu_2) = \frac{(\bar{p}_1 - \bar{p}_2)^2}{2\sigma_p^2},
+$$
+
+which is symmetric in the two means under equal variances.
+
+1. For the "easy" case ($\bar{p}_1 = 2.0$, $\bar{p}_2 = 1.2$) and the "hard"
+ case
+($\bar{p}_1 = 2.0$, $\bar{p}_2 = 1.8$), compute $D_{KL}$ for $\sigma_p = 0.4$.
+
+1. Re-run the simulations from the lecture for both cases with $n=100$ paths.
+ For each
+path compute the first period $T_{0.99}$ at which $h_t \geq 0.99$. Plot
+histograms of
+$T_{0.99}$ for both cases.
+
+1. How does the median $T_{0.99}$ scale with $D_{KL}$? Verify numerically that
+roughly $T_{0.99} \approx C / D_{KL}$ for some constant $C$.
+```
+
+```{solution-start} km_ex2
+:class: dropdown
+```
+
+Here is one solution:
+
+```{code-cell} ipython3
+σ_p = 0.4
+
+def kl_normal(p1, p2, σ):
+ """Return the KL divergence for N(p1, σ^2) and N(p2, σ^2)."""
+ return (p1 - p2)**2 / (2 * σ**2)
+
+cases = [("Easy", 2.0, 1.2), ("Hard", 2.0, 1.8)]
+for name, p1, p2 in cases:
+ kl = kl_normal(p1, p2, σ_p)
+ print(f"{name} case: D_KL = {kl:.4f}")
+
+n_paths = 100
+
+fig, axes = plt.subplots(1, 2, figsize=(11, 4))
+for ax, (name, p1, p2) in zip(axes, cases):
+ kl = kl_normal(p1, p2, σ_p)
+ paths = simulate_bayesian_learning(p1, p2, σ_p, T=2000,
+ h0=0.5, n_paths=n_paths, seed=42)
+ # First period with posterior >= 0.99
+ T99 = []
+ for path in paths:
+ idx = np.where(path >= 0.99)[0]
+ T99.append(idx[0] if len(idx) > 0 else 2001)
+
+ median_T = np.median(T99)
+ ax.hist(T99, bins=20, color="steelblue", edgecolor="white", alpha=0.8)
+ ax.axvline(median_T, color="crimson", lw=2,
+ label=fr"Median $T_{{0.99}} = {median_T:.0f}$")
+ ax.set_title(
+ f"{name}: $D_{{KL}} = {kl:.4f}$, "
+ fr"$\widehat C = T_{{0.99}} D_{{KL}} \approx {median_T * kl:.1f}$",
+ fontsize=11
+ )
+ ax.set_xlabel(r"$T_{0.99}$", fontsize=12)
+ ax.set_ylabel("count", fontsize=11)
+ ax.legend(fontsize=10)
+
+plt.tight_layout()
+plt.show()
+```
+
+The median $T_{0.99}$ scales as approximately $C/D_{KL}$, confirming that
+learning is
+faster when the two reduced forms are more easily distinguished (large
+$D_{KL}$).
+
+```{solution-end}
+```
+
+```{exercise}
+:label: km_ex3
+
+{prf:ref}`ime_theorem_bayesian_convergence` requires the prior to assign
+positive probability to the true reduced-form class $\bar\mu$, equivalently to
+some structure that generates the true price distribution
+$g(\cdot \mid \bar\mu)$.
+
+In this exercise the true reduced form itself is excluded from the prior
+support, so we investigate what happens when no model in the prior generates the
+true price distribution.
+
+Simulate $T = 1,000$ periods of prices from $N(2.0, 0.4^2)$ but use a prior
+that places equal weight on two *wrong* models: $N(1.5, 0.4^2)$ and
+$N(2.3, 0.4^2)$.
+
+Plot the posterior weight on each model over time.
+
+Discuss your findings.
+```
+
+```{solution-start} km_ex3
+:class: dropdown
+```
+
+Here is one solution:
+
+```{code-cell} ipython3
+def simulate_misspecified(
+ T, p_bar_true, p_bar_wrong, σ_p, h0, n_paths, seed=0
+):
+ """Simulate learning under a misspecified two-model prior."""
+ rng = np.random.default_rng(seed)
+ h_paths = np.zeros((n_paths, T + 1, 2))
+ h_paths[:, 0, :] = h0
+
+ for path in range(n_paths):
+ h = np.array(h0, dtype=float)
+ prices = rng.normal(p_bar_true, σ_p, size=T)
+ for t, price in enumerate(prices):
+ likes = norm.pdf(price, loc=p_bar_wrong, scale=σ_p)
+ h = h * likes
+ h /= h.sum()
+ h_paths[path, t + 1, :] = h
+
+ return h_paths
+
+
+def predictive_density(weights, means, σ_p, p_grid):
+ """Return the predictive density under the current posterior weights."""
+ density = np.zeros_like(p_grid)
+ for weight, mean in zip(weights, means):
+ density += weight * norm.pdf(p_grid, loc=mean, scale=σ_p)
+ return density
+
+
+T = 1000
+p_true = 2.0
+p_wrong = np.array([1.5, 2.3])
+σ_p = 0.4
+h0 = np.array([0.5, 0.5])
+n_paths = 30
+
+h_misspec = simulate_misspecified(T, p_true, p_wrong, σ_p, h0, n_paths)
+
+kl_vals = (p_true - p_wrong)**2 / (2 * σ_p**2)
+for mean, kl in zip(p_wrong, kl_vals):
+ print(f"KL(true || N({mean:.1f}, σ^2)) = {kl:.4f}")
+
+t_grid = np.arange(T + 1)
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+
+labels = [r"$N(1.5, \sigma^2)$", r"$N(2.3, \sigma^2)$"]
+for ax, k, label in zip(axes, [0, 1], labels):
+ for path in h_misspec:
+ ax.plot(t_grid, path[:, k], alpha=0.2, lw=0.8, color="steelblue")
+ ax.plot(t_grid, np.median(h_misspec[:, :, k], axis=0),
+ color="navy", lw=2, label="median")
+ ax.set_title(f"Posterior weight on {label}", fontsize=11)
+ ax.set_xlabel("period $t$", fontsize=11)
+ ax.set_ylabel("posterior weight", fontsize=11)
+ ax.legend(fontsize=9)
+
+plt.tight_layout()
+plt.show()
+
+# Predictive density and mean along the median posterior path
+median_path = np.median(h_misspec, axis=0)
+p_grid = np.linspace(0.0, 3.5, 300)
+closer_idx = np.argmin(kl_vals)
+
+fig, ax = plt.subplots(figsize=(8, 4))
+colors = plt.cm.Blues(np.linspace(0.3, 1.0, 4))
+for t_snap, color in zip([0, 10, 100, T], colors):
+ dens = predictive_density(median_path[t_snap], p_wrong, σ_p, p_grid)
+ ax.plot(p_grid, dens, color=color, lw=2, label=f"t = {t_snap}")
+
+ax.plot(
+ p_grid,
+ norm.pdf(p_grid, loc=p_wrong[closer_idx], scale=σ_p),
+ "k--",
+ lw=2,
+ label="KL-best wrong model",
+)
+ax.set_xlabel("price $p$", fontsize=11)
+ax.set_ylabel("density", fontsize=11)
+ax.legend(fontsize=9)
+plt.tight_layout()
+plt.show()
+
+pred_mean = np.median(
+ h_misspec[:, :, 0] * p_wrong[0] + h_misspec[:, :, 1] * p_wrong[1], axis=0
+)
+print(f"True mean: {p_true}")
+print(f"Predictive mean at T={T}: {pred_mean[-1]:.4f}")
+print(f"Closer misspecified mean: {p_wrong[np.argmin(kl_vals)]:.1f}")
+```
+
+Here
+
+$$
+D_{KL}\bigl(N(2.0, 0.4^2)\,\|\,N(2.3, 0.4^2)\bigr)
+<
+D_{KL}\bigl(N(2.0, 0.4^2)\,\|\,N(1.5, 0.4^2)\bigr),
+$$
+
+so the model with mean $2.3$ is the KL-best approximation among the two wrong
+models, and in the simulation posterior weight concentrates on that model.
+
+Posterior odds are cumulative {doc}`likelihood ratios`.
+
+If we compare the two wrong Gaussian models $f$ and $g$, then under the true
+distribution $h$ the average log likelihood ratio satisfies
+
+$$
+\frac{1}{t} E_h[\log L_t] = K(h,g) - K(h,f).
+$$
+
+So if $f$ is KL-closer to $h$ than $g$ is, $\log L_t$ has positive drift and
+posterior odds tilt toward $f$.
+
+```{solution-end}
+```
diff --git a/lectures/networks.md b/lectures/networks.md
new file mode 100644
index 0000000..6ebac5f
--- /dev/null
+++ b/lectures/networks.md
@@ -0,0 +1,1336 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+ jupytext_version: 1.14.4
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+---
+
+# Networks
+
+```{code-cell} ipython3
+:tags: [hide-output]
+
+!pip install quantecon
+!pip install quantecon-book-networks==1.6
+```
+
+## Outline
+
+In recent years there has been rapid growth in a field called [network science](https://en.wikipedia.org/wiki/Network_science).
+
+Network science studies relationships between groups of objects.
+
+One important example is the [world wide web](https://en.wikipedia.org/wiki/World_Wide_Web#Linking)
+, where web pages are connected by hyperlinks.
+
+Another is the [human brain](https://en.wikipedia.org/wiki/Neural_circuit): studies of brain function emphasize the network of
+connections between nerve cells (neurons).
+
+[Artificial neural networks](https://en.wikipedia.org/wiki/Artificial_neural_network) are based on this idea, using data to build
+intricate connections between simple processing units.
+
+Epidemiologists studying [transmission of diseases](https://en.wikipedia.org/wiki/Network_medicine#Network_epidemics)
+like COVID-19 analyze interactions between groups of human hosts.
+
+In operations research, network analysis is used to study fundamental problems such as minimum cost flow, the traveling salesman, [shortest paths](https://en.wikipedia.org/wiki/Shortest_path_problem),
+and assignment.
+
+This lecture gives an introduction to economic and financial networks.
+
+Some parts of this lecture are drawn from the text
+https://networks.quantecon.org/ but the level of this lecture is more
+introductory.
+
+We will need the following imports.
+
+```{code-cell} ipython3
+import numpy as np
+import networkx as nx
+import matplotlib.pyplot as plt
+import pandas as pd
+import quantecon as qe
+
+import matplotlib.cm as cm
+import quantecon_book_networks.input_output as qbn_io
+import quantecon_book_networks.data as qbn_data
+
+import matplotlib.patches as mpatches
+```
+
+## Economic and financial networks
+
+Within economics, important examples of networks include
+
+* financial networks
+* production networks
+* trade networks
+* transport networks and
+* social networks
+
+Social networks affect trends in market sentiment and consumer decisions.
+
+The structure of financial networks helps to determine relative fragility of the financial system.
+
+The structure of production networks affects trade, innovation and the propagation of local shocks.
+
+To better understand such networks, let's look at some examples in more depth.
+
+
+### Example: Aircraft Exports
+
+The following figure shows international trade in large commercial aircraft in 2019 based on International Trade Data SITC Revision 2.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: "Commercial Aircraft Network \n"
+ name: aircraft_network
+tags: [hide-input]
+---
+ch1_data = qbn_data.introduction()
+export_figures = False
+
+DG = ch1_data['aircraft_network']
+pos = ch1_data['aircraft_network_pos']
+
+centrality = nx.eigenvector_centrality(DG)
+node_total_exports = qbn_io.node_total_exports(DG)
+edge_weights = qbn_io.edge_weights(DG)
+
+node_pos_dict = pos
+
+node_sizes = qbn_io.normalise_weights(node_total_exports,10000)
+edge_widths = qbn_io.normalise_weights(edge_weights,10)
+
+node_colors = qbn_io.colorise_weights(list(centrality.values()),color_palette=cm.viridis)
+node_to_color = dict(zip(DG.nodes,node_colors))
+edge_colors = []
+for src,_ in DG.edges:
+ edge_colors.append(node_to_color[src])
+
+fig, ax = plt.subplots(figsize=(10, 10))
+ax.axis('off')
+
+nx.draw_networkx_nodes(DG,
+ node_pos_dict,
+ node_color=node_colors,
+ node_size=node_sizes,
+ linewidths=2,
+ alpha=0.6,
+ ax=ax)
+
+nx.draw_networkx_labels(DG,
+ node_pos_dict,
+ ax=ax)
+
+nx.draw_networkx_edges(DG,
+ node_pos_dict,
+ edge_color=edge_colors,
+ width=edge_widths,
+ arrows=True,
+ arrowsize=20,
+ ax=ax,
+ arrowstyle='->',
+ node_size=node_sizes,
+ connectionstyle='arc3,rad=0.15')
+
+plt.show()
+```
+
+The circles in the figure are called **nodes** or **vertices** -- in this case they represent countries.
+
+The arrows in the figure are called **edges** or **links**.
+
+Node size is proportional to total exports and edge width is proportional to exports to the target country.
+
+(The data is for trade in commercial aircraft weighing at least 15,000kg and was sourced from CID Dataverse.)
+
+The figure shows that the US, France and Germany are major export hubs.
+
+In the discussion below, we learn to quantify such ideas.
+
+
+### Example: A Markov Chain
+
+Recall that, in our lecture on {ref}`Markov chains ` we studied a dynamic model of business cycles
+where the states are
+
+* "ng" = "normal growth"
+* "mr" = "mild recession"
+* "sr" = "severe recession"
+
+Let's examine the following figure
+
+```{image} _static/networks/mc.png
+:name: mc_networks
+:align: center
+```
+
+This is an example of a network, where the set of nodes $V$ equals the states:
+
+$$
+ V = \{ \text{"ng", "mr", "sr"} \}
+$$
+
+The edges between the nodes show the one month transition probabilities.
+
+
+## An introduction to graph theory
+
+Now we've looked at some examples, let's move on to theory.
+
+This theory will allow us to better organize our thoughts.
+
+The theoretical part of network science is constructed using a major branch of
+mathematics called [graph theory](https://en.wikipedia.org/wiki/Graph_theory).
+
+Graph theory can be complicated and we will cover only the basics.
+
+However, these concepts will already be enough for us to discuss interesting and
+important ideas on economic and financial networks.
+
+We focus on "directed" graphs, where connections are, in general, asymmetric
+(arrows typically point one way, not both ways).
+
+E.g.,
+
+* bank $A$ lends money to bank $B$
+* firm $A$ supplies goods to firm $B$
+* individual $A$ "follows" individual $B$ on a given social network
+
+("Undirected" graphs, where connections are symmetric, are a special
+case of directed graphs --- we just need to insist that each arrow pointing
+from $A$ to $B$ is paired with another arrow pointing from $B$ to $A$.)
+
+
+### Key definitions
+
+A **directed graph** consists of two things:
+
+1. a finite set $V$ and
+1. a collection of pairs $(u, v)$ where $u$ and $v$ are elements of $V$.
+
+The elements of $V$ are called the **vertices** or **nodes** of the graph.
+
+The pairs $(u,v)$ are called the **edges** of the graph and the set of all edges will usually be denoted by $E$
+
+Intuitively and visually, an edge $(u,v)$ is understood as an arrow from node $u$ to node $v$.
+
+(A neat way to represent an arrow is to record the location of the tail and
+head of the arrow, and that's exactly what an edge does.)
+
+In the aircraft export example shown in {numref}`aircraft_network`
+
+* $V$ is all countries included in the data set.
+* $E$ is all the arrows in the figure, each indicating some positive amount of aircraft exports from one country to another.
+
+Let's look at more examples.
+
+Two graphs are shown below, each with three nodes.
+
+```{figure} _static/networks/poverty_trap_1.png
+:name: poverty_trap_1
+
+Poverty Trap
+```
+
++++
+
+We now construct a graph with the same nodes but different edges.
+
+```{figure} _static/networks/poverty_trap_2.png
+:name: poverty_trap_2
+
+Poverty Trap
+```
+
++++
+
+For these graphs, the arrows (edges) can be thought of as representing
+positive transition probabilities over a given unit of time.
+
+In general, if an edge $(u, v)$ exists, then the node $u$ is called a
+**direct predecessor** of $v$ and $v$ is called a **direct successor** of $u$.
+
+Also, for $v \in V$,
+
+* the **in-degree** is $i_d(v) = $ the number of direct predecessors of $v$ and
+* the **out-degree** is $o_d(v) = $ the number of direct successors of $v$.
+
+
+### Digraphs in Networkx
+
+The Python package [Networkx](https://networkx.org/) provides a convenient
+data structure for representing directed graphs and implements many common
+routines for analyzing them.
+
+As an example, let us recreate {numref}`poverty_trap_2` using Networkx.
+
+To do so, we first create an empty `DiGraph` object:
+
+```{code-cell} ipython3
+G_p = nx.DiGraph()
+```
+
+Next we populate it with nodes and edges.
+
+To do this we write down a list of
+all edges, with *poor* represented by *p* and so on:
+
+```{code-cell} ipython3
+edge_list = [('p', 'p'),
+ ('m', 'p'), ('m', 'm'), ('m', 'r'),
+ ('r', 'p'), ('r', 'm'), ('r', 'r')]
+```
+
+Finally, we add the edges to our `DiGraph` object:
+
+```{code-cell} ipython3
+for e in edge_list:
+ u, v = e
+ G_p.add_edge(u, v)
+```
+
+Alternatively, we can use the method `add_edges_from`.
+
+```{code-cell} ipython3
+G_p.add_edges_from(edge_list)
+```
+
+Adding the edges automatically adds the nodes, so `G_p` is now a
+correct representation of our graph.
+
+We can verify this by plotting the graph via Networkx with the following code:
+
+```{code-cell} ipython3
+fig, ax = plt.subplots()
+nx.draw_spring(G_p, ax=ax, node_size=500, with_labels=True,
+ font_weight='bold', arrows=True, alpha=0.8,
+ connectionstyle='arc3,rad=0.25', arrowsize=20)
+plt.show()
+```
+
+The figure obtained above matches the original directed graph in {numref}`poverty_trap_2`.
+
+
+`DiGraph` objects have methods that calculate in-degree and out-degree
+of nodes.
+
+For example,
+
+```{code-cell} ipython3
+G_p.in_degree('p')
+```
+(strongly_connected)=
+### Communication
+
+Next, we study communication and connectedness, which have important
+implications for economic networks.
+
+Node $v$ is called **accessible** from node $u$ if either $u=v$ or there
+exists a sequence of edges that lead from $u$ to $v$.
+
+* in this case, we write $u \to v$
+
+(Visually, there is a sequence of arrows leading from $u$ to $v$.)
+
+For example, suppose we have a directed graph representing a production network, where
+
+* elements of $V$ are industrial sectors and
+* existence of an edge $(i, j)$ means that $i$ supplies products or services to $j$.
+
+Then $m \to \ell$ means that sector $m$ is an upstream supplier of sector $\ell$.
+
+Two nodes $u$ and $v$ are said to **communicate** if both $u \to v$ and $v \to u$.
+
+A graph is called **strongly connected** if all nodes communicate.
+
+For example, {numref}`poverty_trap_1` is strongly connected
+however in {numref}`poverty_trap_2` rich is not accessible from poor, thus it is not strongly connected.
+
+We can verify this by first constructing the graphs using Networkx and then using `nx.is_strongly_connected`.
+
+```{code-cell} ipython3
+fig, ax = plt.subplots()
+G1 = nx.DiGraph()
+
+G1.add_edges_from([('p', 'p'),('p','m'),('p','r'),
+ ('m', 'p'), ('m', 'm'), ('m', 'r'),
+ ('r', 'p'), ('r', 'm'), ('r', 'r')])
+
+nx.draw_networkx(G1, with_labels = True)
+```
+
+```{code-cell} ipython3
+nx.is_strongly_connected(G1) #checking if above graph is strongly connected
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots()
+G2 = nx.DiGraph()
+
+G2.add_edges_from([('p', 'p'),
+ ('m', 'p'), ('m', 'm'), ('m', 'r'),
+ ('r', 'p'), ('r', 'm'), ('r', 'r')])
+
+nx.draw_networkx(G2, with_labels = True)
+```
+
+```{code-cell} ipython3
+nx.is_strongly_connected(G2) #checking if above graph is strongly connected
+```
+
+## Weighted graphs
+
+We now introduce weighted graphs, where weights (numbers) are attached to each
+edge.
+
+
+### International private credit flows by country
+
+To motivate the idea, consider the following figure which shows flows of funds (i.e.,
+loans) between private banks, grouped by country of origin.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: "International Credit Network \n"
+ name: financial_network
+tags: [hide-input]
+---
+Z = ch1_data["adjacency_matrix"]["Z"]
+Z_visual= ch1_data["adjacency_matrix"]["Z_visual"]
+countries = ch1_data["adjacency_matrix"]["countries"]
+
+G = qbn_io.adjacency_matrix_to_graph(Z_visual, countries, tol=0.03)
+
+centrality = qbn_io.eigenvector_centrality(Z_visual, authority=False)
+node_total_exports = qbn_io.node_total_exports(G)
+edge_weights = qbn_io.edge_weights(G)
+
+node_pos_dict = nx.circular_layout(G)
+
+node_sizes = qbn_io.normalise_weights(node_total_exports,3000)
+edge_widths = qbn_io.normalise_weights(edge_weights,10)
+
+
+node_colors = qbn_io.colorise_weights(centrality)
+node_to_color = dict(zip(G.nodes,node_colors))
+edge_colors = []
+for src,_ in G.edges:
+ edge_colors.append(node_to_color[src])
+
+fig, ax = plt.subplots(figsize=(10, 10))
+ax.axis('off')
+
+nx.draw_networkx_nodes(G,
+ node_pos_dict,
+ node_color=node_colors,
+ node_size=node_sizes,
+ edgecolors='grey',
+ linewidths=2,
+ alpha=0.4,
+ ax=ax)
+
+nx.draw_networkx_labels(G,
+ node_pos_dict,
+ font_size=12,
+ ax=ax)
+
+nx.draw_networkx_edges(G,
+ node_pos_dict,
+ edge_color=edge_colors,
+ width=edge_widths,
+ arrows=True,
+ arrowsize=20,
+ alpha=0.8,
+ ax=ax,
+ arrowstyle='->',
+ node_size=node_sizes,
+ connectionstyle='arc3,rad=0.15')
+
+plt.show()
+```
+
+The country codes are given in the following table
+
+|Code| Country |Code| Country |Code| Country |Code| Country |
+|:--:|:--------------|:--:|:-------:|:--:|:-----------:|:--:|:--------------:|
+| AU | Australia | DE | Germany | CL | Chile | ES | Spain |
+| PT | Portugal | FR | France | TR | Turkey | GB | United Kingdom |
+| US | United States | IE | Ireland | AT | Austria | IT | Italy |
+| BE | Belgium | JP | Japan | SW | Switzerland | SE | Sweden |
+
+An arrow from Japan to the US indicates aggregate claims held by Japanese
+banks on all US-registered banks, as collected by the Bank of International
+Settlements (BIS).
+
+The size of each node in the figure is increasing in the
+total foreign claims of all other nodes on this node.
+
+The widths of the arrows are proportional to the foreign claims they represent.
+
+Notice that, in this network, an edge $(u, v)$ exists for almost every choice
+of $u$ and $v$ (i.e., almost every country in the network).
+
+(In fact, there are even more small arrows, which we have dropped for clarity.)
+
+Hence the existence of an edge from one node to another is not particularly informative.
+
+To understand the network, we need to record not just the existence or absence
+of a credit flow, but also the size of the flow.
+
+The correct data structure for recording this information is a "weighted
+directed graph".
+
++++
+
+### Definitions
+
+A **weighted directed graph** is a directed graph to which we have added a
+**weight function** $w$ that assigns a positive number to each edge.
+
+The figure above shows one weighted directed graph, where the weights are the size of fund flows.
+
+The following figure shows a weighted directed graph, with arrows
+representing edges of the induced directed graph.
+
+
+```{figure} _static/networks/weighted.png
+:name: poverty_trap_weighted
+
+Weighted Poverty Trap
+```
+
+
+The numbers next to the edges are the weights.
+
+In this case, you can think of the numbers on the arrows as transition
+probabilities for a household over, say, one year.
+
+We see that a rich household has a 10\% chance of becoming poor in one year.
+
+
+## Adjacency matrices
+
+Another way that we can represent weights, which turns out to be very
+convenient for numerical work, is via a matrix.
+
+The **adjacency matrix** of a weighted directed graph with nodes $\{v_1, \ldots, v_n\}$, edges $E$ and weight function $w$ is the matrix
+
+$$
+A = (a_{ij})_{1 \leq i,j \leq n}
+\quad \text{with} \quad
+a_{ij} =
+%
+\begin{cases}
+ w(v_i, v_j) & \text{ if } (v_i, v_j) \in E
+ \\
+ 0 & \text{ otherwise}.
+\end{cases}
+%
+$$
+
+Once the nodes in $V$ are enumerated, the weight function and
+adjacency matrix provide essentially the same information.
+
+For example, with $\{$poor, middle, rich$\}$ mapped to $\{1, 2, 3\}$ respectively,
+the adjacency matrix corresponding to the weighted directed graph in {numref}`poverty_trap_weighted` is
+
+$$
+\begin{pmatrix}
+ 0.9 & 0.1 & 0 \\
+ 0.4 & 0.4 & 0.2 \\
+ 0.1 & 0.1 & 0.8
+\end{pmatrix}.
+$$
+
+In QuantEcon's `DiGraph` implementation, weights are recorded via the
+keyword `weighted`:
+
+```{code-cell} ipython3
+A = ((0.9, 0.1, 0.0),
+ (0.4, 0.4, 0.2),
+ (0.1, 0.1, 0.8))
+A = np.array(A)
+G = qe.DiGraph(A, weighted=True) # store weights
+```
+
+One of the key points to remember about adjacency matrices is that taking the
+transpose _reverses all the arrows_ in the associated directed graph.
+
+
+For example, the following directed graph can be
+interpreted as a stylized version of a financial network, with nodes as banks
+and edges showing the flow of funds.
+
+```{code-cell} ipython3
+G4 = nx.DiGraph()
+
+G4.add_edges_from([('1','2'),
+ ('2','1'),('2','3'),
+ ('3','4'),
+ ('4','2'),('4','5'),
+ ('5','1'),('5','3'),('5','4')])
+pos = nx.circular_layout(G4)
+
+edge_labels={('1','2'): '100',
+ ('2','1'): '50', ('2','3'): '200',
+ ('3','4'): '100',
+ ('4','2'): '500', ('4','5'): '50',
+ ('5','1'): '150',('5','3'): '250', ('5','4'): '300'}
+
+nx.draw_networkx(G4, pos, node_color = 'none',node_size = 500)
+nx.draw_networkx_edge_labels(G4, pos, edge_labels=edge_labels)
+nx.draw_networkx_nodes(G4, pos, linewidths= 0.5, edgecolors = 'black',
+ node_color = 'none',node_size = 500)
+
+plt.show()
+```
+
+We see that bank 2 extends a loan of size 200 to bank 3.
+
+The corresponding adjacency matrix is
+
+$$
+A =
+\begin{pmatrix}
+ 0 & 100 & 0 & 0 & 0 \\
+ 50 & 0 & 200 & 0 & 0 \\
+ 0 & 0 & 0 & 100 & 0 \\
+ 0 & 500 & 0 & 0 & 50 \\
+ 150 & 0 & 250 & 300 & 0
+\end{pmatrix}.
+$$
+
+The transpose is
+
+$$
+A^\top =
+\begin{pmatrix}
+ 0 & 50 & 0 & 0 & 150 \\
+ 100 & 0 & 0 & 500 & 0 \\
+ 0 & 200 & 0 & 0 & 250 \\
+ 0 & 0 & 100 & 0 & 300 \\
+ 0 & 0 & 0 & 50 & 0
+\end{pmatrix}.
+$$
+
+The corresponding network is visualized in the following figure which shows the network of liabilities after the loans have been granted.
+
+Both of these networks (original and transpose) are useful for analyzing financial markets.
+
+```{code-cell} ipython3
+G5 = nx.DiGraph()
+
+G5.add_edges_from([('1','2'),('1','5'),
+ ('2','1'),('2','4'),
+ ('3','2'),('3','5'),
+ ('4','3'),('4','5'),
+ ('5','4')])
+
+edge_labels={('1','2'): '50', ('1','5'): '150',
+ ('2','1'): '100', ('2','4'): '500',
+ ('3','2'): '200', ('3','5'): '250',
+ ('4','3'): '100', ('4','5'): '300',
+ ('5','4'): '50'}
+
+nx.draw_networkx(G5, pos, node_color = 'none',node_size = 500)
+nx.draw_networkx_edge_labels(G5, pos, edge_labels=edge_labels)
+nx.draw_networkx_nodes(G5, pos, linewidths= 0.5, edgecolors = 'black',
+ node_color = 'none',node_size = 500)
+
+plt.show()
+```
+
+In general, every nonnegative $n \times n$ matrix $A = (a_{ij})$ can be
+viewed as the adjacency matrix of a weighted directed graph.
+
+To build the graph we set $V = 1, \ldots, n$ and take the edge set $E$ to be
+all $(i,j)$ such that $a_{ij} > 0$.
+
+For the weight function we set $w(i, j) = a_{ij}$ for all edges $(i,j)$.
+
+We call this graph the weighted directed graph induced by $A$.
+
+
+## Properties
+
+Consider a weighted directed graph with adjacency matrix $A$.
+
+Let $a^k_{ij}$ be element $i,j$ of $A^k$, the $k$-th power of $A$.
+
+The following result is useful in many applications:
+
+````{prf:theorem}
+:label: graph_theory_property1
+
+For distinct nodes $i, j$ in $V$ and any integer $k$, we have
+
+$$
+a^k_{i j} > 0
+\quad \text{if and only if} \quad
+\text{ $j$ is accessible from $i$}.
+$$
+
+````
+
++++
+
+The above result is obvious when $k=1$ and a proof of the general case can be
+found in {cite}`sargent2022economic`.
+
+Now recall from the eigenvalues lecture that a
+nonnegative matrix $A$ is called {ref}`irreducible` if for each $(i,j)$ there is an integer $k \geq 0$ such that $a^{k}_{ij} > 0$.
+
+From the preceding theorem, it is not too difficult (see
+{cite}`sargent2022economic` for details) to get the next result.
+
+````{prf:theorem}
+:label: graph_theory_property2
+
+For a weighted directed graph the following statements are equivalent:
+
+1. The directed graph is strongly connected.
+2. The adjacency matrix of the graph is irreducible.
+
+````
+
++++
+
+We illustrate the above theorem with a simple example.
+
+Consider the following weighted directed graph.
+
+
+```{image} _static/networks/properties.png
+:name: properties_graph
+
+```
+
++++
+
+We first create the above network as a Networkx `DiGraph` object.
+
+```{code-cell} ipython3
+G6 = nx.DiGraph()
+
+G6.add_edges_from([('1','2'),('1','3'),
+ ('2','1'),
+ ('3','1'),('3','2')])
+```
+
+Then we construct the associated adjacency matrix A.
+
+```{code-cell} ipython3
+A = np.array([[0,0.7,0.3], # adjacency matrix A
+ [1,0,0],
+ [0.4,0.6,0]])
+```
+
+```{code-cell} ipython3
+:tags: [hide-input]
+
+def is_irreducible(P):
+ n = len(P)
+ result = np.zeros((n, n))
+ for i in range(n):
+ result += np.linalg.matrix_power(P, i)
+ return np.all(result > 0)
+```
+
+```{code-cell} ipython3
+is_irreducible(A) # check irreducibility of A
+```
+
+```{code-cell} ipython3
+nx.is_strongly_connected(G6) # check connectedness of graph
+```
+
+## Network centrality
+
+When studying networks of all varieties, a recurring topic is the relative
+"centrality" or "importance" of different nodes.
+
+Examples include
+
+* ranking of web pages by search engines
+* determining the most important bank in a financial network (which one a
+ central bank should rescue if there is a financial crisis)
+* determining the most important industrial sector in an economy.
+
+In what follows, a **centrality measure** associates to each weighted directed
+graph a vector $m$ where the $m_i$ is interpreted as the centrality (or rank)
+of node $v_i$.
+
+### Degree centrality
+
+Two elementary measures of "importance" of a node in a given directed
+graph are its in-degree and out-degree.
+
+Both of these provide a centrality measure.
+
+In-degree centrality is a vector containing the in-degree of each node in
+the graph.
+
+Consider the following simple example.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Sample Graph
+ name: sample_gph_1
+---
+G7 = nx.DiGraph()
+
+G7.add_nodes_from(['1','2','3','4','5','6','7'])
+
+G7.add_edges_from([('1','2'),('1','6'),
+ ('2','1'),('2','4'),
+ ('3','2'),
+ ('4','2'),
+ ('5','3'),('5','4'),
+ ('6','1'),
+ ('7','4'),('7','6')])
+pos = nx.planar_layout(G7)
+
+nx.draw_networkx(G7, pos, node_color='none', node_size=500)
+nx.draw_networkx_nodes(G7, pos, linewidths=0.5, edgecolors='black',
+ node_color='none',node_size=500)
+
+plt.show()
+```
+
+The following code displays the in-degree centrality of all nodes.
+
+```{code-cell} ipython3
+iG7 = [G7.in_degree(v) for v in G7.nodes()] # computing in-degree centrality
+
+for i, d in enumerate(iG7):
+ print(i+1, d)
+```
+
+Consider the international credit network displayed in {numref}`financial_network`.
+
+The following plot displays the in-degree centrality of each country.
+
+```{code-cell} ipython3
+D = qbn_io.build_unweighted_matrix(Z)
+indegree = D.sum(axis=0)
+```
+
+```{code-cell} ipython3
+def centrality_plot_data(countries, centrality_measures):
+ df = pd.DataFrame({'code': countries,
+ 'centrality':centrality_measures,
+ 'color': qbn_io.colorise_weights(centrality_measures).tolist()
+ })
+ return df.sort_values('centrality')
+```
+
+```{code-cell} ipython3
+fig, ax = plt.subplots()
+
+df = centrality_plot_data(countries, indegree)
+
+ax.bar('code', 'centrality', data=df, color=df["color"], alpha=0.6)
+
+patch = mpatches.Patch(color=None, label='in degree', visible=False)
+ax.legend(handles=[patch], fontsize=12, loc="upper left", handlelength=0, frameon=False)
+
+ax.set_ylim((0,20))
+
+plt.show()
+```
+
+Unfortunately, while in-degree and out-degree centrality are simple to
+calculate, they are not always informative.
+
+In {numref}`financial_network`, an edge exists between almost every node,
+so the in- or out-degree based centrality ranking fails to effectively separate the countries.
+
+This can be seen in the above graph as well.
+
+Another example is the task of a web search engine, which ranks pages
+by relevance whenever a user enters a search.
+
+Suppose web page A has twice as many inbound links as page B.
+
+In-degree centrality tells us that page A deserves a higher rank.
+
+But in fact, page A might be less important than page B.
+
+To see why, suppose that the links to A are from pages that receive almost no traffic,
+while the links to B are from pages that receive very heavy traffic.
+
+In this case, page B probably receives more visitors, which in turn suggests
+that page B contains more valuable (or entertaining) content.
+
+Thinking about this point suggests that importance might be *recursive*.
+
+This means that the importance of a given node depends on the importance of
+other nodes that link to it.
+
+As another example, we can imagine a production network where the importance of a
+given sector depends on the importance of the sectors that it supplies.
+
+This reverses the order of the previous example: now the importance of a given
+node depends on the importance of other nodes that *it links to*.
+
+The next centrality measures will have these recursive features.
+
+
+### Eigenvector centrality
+
+Suppose we have a weighted directed graph with adjacency matrix $A$.
+
+For simplicity, we will suppose that the nodes $V$ of the graph are just the
+integers $1, \ldots, n$.
+
+Let $r(A)$ denote the {ref}`spectral radius` of $A$.
+
+The **eigenvector centrality** of the graph is defined as the $n$-vector $e$ that solves
+
+$$
+\begin{aligned}
+ e = \frac{1}{r(A)} A e.
+\end{aligned}
+$$ (ev_central)
+
+In other words, $e$ is the dominant eigenvector of $A$ (the eigenvector of the
+largest eigenvalue --- see the discussion of the {ref}`Perron-Frobenius theorem` in the eigenvalue lecture.
+
+To better understand {eq}`ev_central`, we write out the full expression
+for some element $e_i$
+
+$$
+\begin{aligned}
+ e_i = \frac{1}{r(A)} \sum_{1 \leq j \leq n} a_{ij} e_j
+\end{aligned}
+$$ (eq_eicen)
+
+
+Note the recursive nature of the definition: the centrality obtained by node
+$i$ is proportional to a sum of the centrality of all nodes, weighted by
+the *rates of flow* from $i$ into these nodes.
+
+A node $i$ is highly ranked if
+1. there are many edges leaving $i$,
+2. these edges have large weights, and
+3. the edges point to other highly ranked nodes.
+
+Later, when we study demand shocks in production networks, there will be a more
+concrete interpretation of eigenvector centrality.
+
+We will see that, in production networks, sectors with high eigenvector
+centrality are important *suppliers*.
+
+In particular, they are activated by a wide array of demand shocks once orders
+flow backwards through the network.
+
+To compute eigenvector centrality we can use the following function.
+
+```{code-cell} ipython3
+def eigenvector_centrality(A, k=40, authority=False):
+ """
+ Computes the dominant eigenvector of A. Assumes A is
+ primitive and uses the power method.
+
+ """
+ A_temp = A.T if authority else A
+ n = len(A_temp)
+ r = np.max(np.abs(np.linalg.eigvals(A_temp)))
+ e = r**(-k) * (np.linalg.matrix_power(A_temp, k) @ np.ones(n))
+ return e / np.sum(e)
+```
+
+Let's compute eigenvector centrality for the graph generated in {numref}`sample_gph_1`.
+
+```{code-cell} ipython3
+A = nx.to_numpy_array(G7) # compute adjacency matrix of graph
+```
+
+```{code-cell} ipython3
+e = eigenvector_centrality(A)
+n = len(e)
+
+for i in range(n):
+ print(i+1,e[i])
+```
+
+While nodes $2$ and $4$ had the highest in-degree centrality, we can see that nodes $1$ and $2$ have the
+highest eigenvector centrality.
+
+Let's revisit the international credit network in {numref}`financial_network`.
+
+```{code-cell} ipython3
+eig_central = eigenvector_centrality(Z)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Eigenvector centrality
+ name: eigenvctr_centrality
+---
+fig, ax = plt.subplots()
+
+df = centrality_plot_data(countries, eig_central)
+
+ax.bar('code', 'centrality', data=df, color=df["color"], alpha=0.6)
+
+patch = mpatches.Patch(color=None, visible=False)
+ax.legend(handles=[patch], fontsize=12, loc="upper left", handlelength=0, frameon=False)
+
+plt.show()
+```
+
+Countries that are rated highly according to this rank tend to be important
+players in terms of supply of credit.
+
+Japan takes the highest rank according to this measure, although
+countries with large financial sectors such as Great Britain and France are
+not far behind.
+
+The advantage of eigenvector centrality is that it measures a node's importance while considering the importance of its neighbours.
+
+A variant of eigenvector centrality is at the core of Google's PageRank algorithm, which is used to rank web pages.
+
+The main principle is that links from important nodes (as measured by degree centrality) are worth more than links from unimportant nodes.
+
+
+### Katz centrality
+
+One problem with eigenvector centrality is that $r(A)$ might be zero, in which
+case $1/r(A)$ is not defined.
+
+For this and other reasons, some researchers prefer another measure of
+centrality for networks called Katz centrality.
+
+Fixing $\beta$ in $(0, 1/r(A))$, the **Katz centrality** of a weighted
+directed graph with adjacency matrix $A$ is defined as the vector $\kappa$
+that solves
+
+$$
+\kappa_i = \beta \sum_{j=0}^{n-1} a_{ij} \kappa_j + 1
+\qquad \text{for all } i \in \{0, \ldots, n-1\}.
+$$ (katz_central)
+
+Here $\beta$ is a parameter that we can choose.
+
+In vector form we can write
+
+$$
+\kappa = \mathbf 1 + \beta A \kappa
+$$ (katz_central_vec)
+
+where $\mathbf 1$ is a column vector of ones.
+
+The intuition behind this centrality measure is similar to that provided for
+eigenvector centrality: high centrality is conferred on $i$ when it is linked
+to by nodes that themselves have high centrality.
+
+Provided that $0 < \beta < 1/r(A)$, Katz centrality is always finite and well-defined
+because then $r(\beta A) < 1$.
+
+This means that {eq}`katz_central_vec` has the unique solution
+
+$$
+\kappa = (I - \beta A)^{-1} \mathbf{1}
+$$
+
+
+This follows from the {ref}`Neumann series theorem`.
+
+The parameter $\beta$ is used to ensure that $\kappa$ is finite
+
+When $r(A)<1$, we use $\beta=1$ as the default for Katz centrality computations.
+
+
+### Authorities vs hubs
+
+Search engine designers recognize that web pages can be important in two
+different ways.
+
+Some pages have high **hub centrality**, meaning that they link to valuable
+sources of information (e.g., news aggregation sites).
+
+Other pages have high **authority centrality**, meaning that they contain
+valuable information, as indicated by the number and significance of incoming
+links (e.g., websites of respected news organizations).
+
+Similar ideas can and have been applied to economic networks (often using
+different terminology).
+
+The eigenvector centrality and Katz centrality measures we discussed above
+measure hub centrality.
+
+(Nodes have high centrality if they point to other nodes with high centrality.)
+
+If we care more about authority centrality, we can use the same definitions
+except that we take the transpose of the adjacency matrix.
+
+This works because taking the transpose reverses the direction of the arrows.
+
+(Now nodes will have high centrality if they receive links from other nodes
+with high centrality.)
+
+For example, the **authority-based eigenvector centrality** of a weighted
+directed graph with adjacency matrix $A$ is the vector $e$ solving
+
+$$
+e = \frac{1}{r(A)} A^\top e.
+$$ (eicena0)
+
+The only difference from the original definition is that $A$ is replaced by
+its transpose.
+
+(Transposes do not affect the spectral radius of a matrix so we wrote $r(A)$ instead of $r(A^\top)$.)
+
+Element-by-element, this is given by
+
+$$
+e_j = \frac{1}{r(A)} \sum_{1 \leq i \leq n} a_{ij} e_i
+$$ (eicena)
+
+We see $e_j$ will be high if many nodes with high authority rankings link to $j$.
+
+The following figure shows the authority-based eigenvector centrality ranking for the international
+credit network shown in {numref}`financial_network`.
+
+```{code-cell} ipython3
+ecentral_authority = eigenvector_centrality(Z, authority=True)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Eigenvector authority
+ name: eigenvector_centrality
+---
+fig, ax = plt.subplots()
+
+df = centrality_plot_data(countries, ecentral_authority)
+
+ax.bar('code', 'centrality', data=df, color=df["color"], alpha=0.6)
+
+patch = mpatches.Patch(color=None, visible=False)
+ax.legend(handles=[patch], fontsize=12, loc="upper left", handlelength=0, frameon=False)
+
+plt.show()
+```
+
+Highly ranked countries are those that attract large inflows of credit, or
+credit inflows from other major players.
+
+In this case the US clearly dominates the rankings as a target of interbank credit.
+
+
+## Further reading
+
+We apply the ideas discussed in this lecture to:
+
+Textbooks on economic and social networks include {cite}`jackson2010social`,
+{cite}`easley2010networks`, {cite}`borgatti2018analyzing`,
+{cite}`sargent2022economic` and {cite}`goyal2023networks`.
+
+
+Within the realm of network science, the texts
+by {cite}`newman2018networks`, {cite}`menczer2020first` and
+{cite}`coscia2021atlas` are excellent.
+
+
+## Exercises
+
+```{exercise-start}
+:label: networks_ex1
+```
+
+Here is a mathematical exercise for those who like proofs.
+
+Let $(V, E)$ be a directed graph and write $u \sim v$ if $u$ and $v$ communicate.
+
+Show that $\sim$ is an [equivalence relation](https://en.wikipedia.org/wiki/Equivalence_relation) on $V$.
+
+```{exercise-end}
+```
+
+```{solution-start} networks_ex1
+:class: dropdown
+```
+
+**Reflexivity:**
+
+Trivially, $u = v \Rightarrow u \rightarrow v$.
+
+Thus, $u \sim u$.
+
+**Symmetry:**
+Suppose, $u \sim v$
+
+$\Rightarrow u \rightarrow v$ and $v \rightarrow u$.
+
+By definition, this implies $v \sim u$.
+
+**Transitivity:**
+
+Suppose, $u \sim v$ and $v \sim w$
+
+This implies, $u \rightarrow v$ and $v \rightarrow u$ and also $v \rightarrow w$ and $w \rightarrow v$.
+
+Thus, we can conclude $u \rightarrow v \rightarrow w$ and $w \rightarrow v \rightarrow u$.
+
+Which means $u \sim w$.
+
+```{solution-end}
+```
+
+```{exercise-start}
+:label: networks_ex2
+```
+
+Consider a directed graph $G$ with the set of nodes
+
+$$
+V = \{0,1,2,3,4,5,6,7\}
+$$
+
+and the set of edges
+
+$$
+E = \{(0, 1), (0, 3), (1, 0), (2, 4), (3, 2), (3, 4), (3, 7), (4, 3), (5, 4), (5, 6), (6, 3), (6, 5), (7, 0)\}
+$$
+
+1. Use `Networkx` to draw graph $G$.
+
+2. Find the associated adjacency matrix $A$ for $G$.
+
+3. Use the functions defined above to compute in-degree centrality, out-degree centrality and eigenvector centrality
+ of G.
+
+```{exercise-end}
+```
+
+```{solution-start} networks_ex2
+:class: dropdown
+```
+
+```{code-cell} ipython3
+# First, let's plot the given graph
+
+G = nx.DiGraph()
+
+G.add_nodes_from(np.arange(8)) # adding nodes
+
+G.add_edges_from([(0,1),(0,3), # adding edges
+ (1,0),
+ (2,4),
+ (3,2),(3,4),(3,7),
+ (4,3),
+ (5,4),(5,6),
+ (6,3),(6,5),
+ (7,0)])
+
+nx.draw_networkx(G, pos=nx.circular_layout(G), node_color='gray', node_size=500, with_labels=True)
+
+plt.show()
+```
+
+```{code-cell} ipython3
+A = nx.to_numpy_array(G) #find adjacency matrix associated with G
+
+A
+```
+
+```{code-cell} ipython3
+oG = [G.out_degree(v) for v in G.nodes()] # computing in-degree centrality
+
+for i, d in enumerate(oG):
+ print(i, d)
+```
+
+```{code-cell} ipython3
+e = eigenvector_centrality(A) # computing eigenvector centrality
+n = len(e)
+
+for i in range(n):
+ print(i+1, e[i])
+```
+
+```{solution-end}
+```
+
+```{exercise-start}
+:label: networks_ex3
+```
+
+Consider a graph $G$ with $n$ nodes and $n \times n$ adjacency matrix $A$.
+
+Let $S = \sum_{k=0}^{n-1} A^k$
+
+We can say for any two nodes $i$ and $j$, $j$ is accessible from $i$ if and only if
+$S_{ij} > 0$.
+
+Devise a function `is_accessible` that checks if any two nodes of a given graph are accessible.
+
+Consider the graph in {ref}`networks_ex2` and use this function to check if
+
+1. $1$ is accessible from $2$
+2. $6$ is accessible from $3$
+
+```{exercise-end}
+```
+
+```{solution-start} networks_ex3
+:class: dropdown
+```
+
+```{code-cell} ipython3
+def is_accessible(G,i,j):
+ A = nx.to_numpy_array(G)
+ n = len(A)
+ result = np.zeros((n, n))
+ for i in range(n):
+ result += np.linalg.matrix_power(A, i)
+ if result[i,j]>0:
+ return True
+ else:
+ return False
+```
+
+```{code-cell} ipython3
+G = nx.DiGraph()
+
+G.add_nodes_from(np.arange(8)) # adding nodes
+
+G.add_edges_from([(0,1),(0,3), # adding edges
+ (1,0),
+ (2,4),
+ (3,2),(3,4),(3,7),
+ (4,3),
+ (5,4),(5,6),
+ (6,3),(6,5),
+ (7,0)])
+```
+
+```{code-cell} ipython3
+is_accessible(G, 2, 1)
+```
+
+```{code-cell} ipython3
+is_accessible(G, 3, 6)
+```
+
+```{solution-end}
+```
diff --git a/lectures/phillips_drifts_volatilities.md b/lectures/phillips_drifts_volatilities.md
new file mode 100644
index 0000000..117478e
--- /dev/null
+++ b/lectures/phillips_drifts_volatilities.md
@@ -0,0 +1,3413 @@
+---
+jupytext:
+ text_representation:
+ extension: .md
+ format_name: myst
+ format_version: 0.13
+ jupytext_version: 1.16.7
+kernelspec:
+ display_name: Python 3 (ipykernel)
+ language: python
+ name: python3
+---
+
+(phillips_drifts_volatilities)=
+```{raw} jupyter
+
+```
+
+# Drifts and Volatilities
+
+```{index} single: Phillips Curve; Drifts and Volatilities
+```
+
+```{contents} Contents
+:depth: 2
+```
+
+## Overview
+
+The lectures in this section have told a story about how a government's *model*
+of the Phillips curve, and the *policy* it induces, can drift over time.
+
+In {doc}`phillips_learning` and {doc}`phillips_escaping_nash` a government that
+fits and refits an approximating Phillips curve is repeatedly pushed away from a
+bad {doc}`self-confirming equilibrium ` along an
+*escape route*, while {doc}`phillips_priors` and {doc}`phillips_lost_conquest`
+use drifting beliefs to interpret the rise and fall of American inflation.
+
+Those lectures were mostly about *theory*.
+
+This lecture turns to the *data*.
+
+It studies {cite:t}`CogleySargent2005`, which asks a deceptively simple
+question:
+
+> When we look at postwar U.S. time series on inflation, unemployment, and
+> interest rates, do we see evidence that the dynamics have *drifted*?
+
+Tim Cogley and Thomas Sargent began this work as an empirical companion to the
+*Conquest* book {cite}`Sargent1999` and the escape-route papers
+{cite}`ChoWilliamsSargent2002`.
+
+It is also a response to searching comments by {cite:t}`Sims2001comment` and
+{cite:t}`Stock2001comment` on an earlier paper {cite}`CogleySargent2001`, and it
+grew into a friendly debate with {cite:t}`SimsZha2006` and
+{cite:t}`BernankeMihov1998` about a question that organizes this whole section:
+
+*Was the Great Inflation of the 1970s and its conquest in the 1980s a story of
+bad policy, or of bad luck?*
+
+To let the data speak to that question we need a statistical model flexible
+enough to accommodate *both* answers.
+
+That model is a *Bayesian vector autoregression whose coefficients drift as
+random walks and whose shock variances evolve as stochastic volatilities*.
+
+Fitting it requires a Markov chain Monte Carlo algorithm that combines the
+{doc}`Kalman filter `, the forward-filter/backward-sample smoother of
+{cite:t}`CarterKohn1994`, and the stochastic-volatility sampler of
+{cite:t}`Jacquier1994`.
+
+Readers who want background will find companion-form vector autoregressions in
+{doc}`var_dmd`, the Kalman smoother in {doc}`kalman_2`, and Bayesian inference
+for state-space models by MCMC in {doc}`ar1_bayes` and {doc}`ar1_turningpts`.
+
+We work through the data transformation, prior, sampler, and main empirical
+results.
+
+Let's start with some imports and the URL for the data.
+
+```{code-cell} ipython3
+import time
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from IPython.display import display, Math
+from scipy import linalg
+from scipy.special import expit
+from scipy.stats import invwishart
+
+
+data_url = 'https://github.com/QuantEcon/data-lectures/raw/main/lectures/NEWQDATA.csv'
+```
+
+## Bad policy or bad luck?
+
+Two respectable views compete to explain the American Great Inflation — the same
+two stories, triumph versus vindication, that open {doc}`phillips_two_stories`.
+
+The **bad policy** view is the one dramatized throughout this section and in the
+*Conquest* book {cite}`Sargent1999`.
+
+Something about Arthur Burns's *model* of the economy, his *patience*, or his
+inability to *commit* to a better rule led the Federal Reserve to administer
+monetary policy in a way that produced the greatest peacetime inflation in U.S.
+history, while an improved model, more patience, or greater discipline led Paul
+Volcker to conquer it {cite}`DeLong1997,Taylor1997comment`.
+
+On this view, what changed between the 1970s and the 1980s was the *systematic
+part* of policy, namely the way the Fed's interest-rate setting responded to
+inflation and unemployment.
+
+The **bad luck** view says something quite different.
+
+What distinguished the Burns and Volcker eras was not their models or policies,
+but the *shocks* that hit the economy.
+
+On this view the *coefficients* of a reduced-form description of the economy
+were essentially constant, and what changed was the *size* of the disturbances,
+namely the *volatility*.
+
+{cite:t}`BernankeMihov1998` and {cite:t}`SimsZha2006` marshaled evidence for
+this second view, in part by applying classical tests that *failed to reject*
+the hypothesis that VAR coefficients were time invariant.
+
+How can we discriminate?
+
+A constant-coefficient, constant-volatility VAR can generate unusually large
+realized shocks, but it cannot represent systematic changes in their variance
+over time.
+
+A model with drifting coefficients but constant volatility can also mistake
+changing volatility for coefficient drift.
+
+So Cogley and Sargent build a model that has room for *both* channels at once,
+and they let a Bayesian posterior sort out how much of each the data call for.
+
+(csdv-model)=
+## A VAR with drifting coefficients and stochastic volatility
+
+Let the variables be ordered as nominal interest, transformed unemployment, and
+inflation,
+
+$$
+y_t = \begin{bmatrix} i_t & u_t & \pi_t \end{bmatrix}^\top.
+$$
+
+(Here $u_t$ is not the raw unemployment rate but its logit, we define this transformation in the data section below.)
+
+The measurement equation is a VAR with two lags and date-specific coefficients,
+
+```{math}
+:label: csdv_measurement
+y_t = X_t^\top \theta_t + \varepsilon_t,
+\qquad
+X_t^\top = I_3 \otimes \begin{bmatrix} 1 & y_{t-1}^\top & y_{t-2}^\top \end{bmatrix}.
+```
+
+Each equation has an intercept and six lag coefficients, so $\theta_t$ contains
+$3(1+2\times 3)=21$ elements.
+
+A two-lag VAR can be rewritten as a one-lag system by stacking $y_t$ and
+$y_{t-1}$ into a single vector; the matrix that multiplies this stacked vector
+is the **companion matrix**, and the rewritten system is the VAR in
+**companion form**.
+
+The following function builds the companion matrix from the stacked
+coefficients.
+
+```{code-cell} ipython3
+n_variables = 3
+n_lags = 2
+n_regressors = 1 + n_variables * n_lags
+n_coefficients = n_variables * n_regressors
+
+
+def companion_matrix(θ):
+ """Return the intercept and companion matrix for one coefficient vector."""
+ equation_rows = np.asarray(θ, dtype=float).reshape(
+ n_variables, n_regressors
+ )
+ intercept = np.r_[equation_rows[:, 0], np.zeros(n_variables)]
+ companion = np.zeros((n_variables * n_lags, n_variables * n_lags))
+ companion[:n_variables] = equation_rows[:, 1:]
+ companion[n_variables:, :n_variables] = np.eye(n_variables)
+ return intercept, companion
+
+
+def design_matrix(regressors):
+ """Return the observation matrix X_t prime for one date."""
+ return np.kron(np.eye(n_variables), np.asarray(regressors, dtype=float))
+```
+
+The coefficient vector follows a driftless random walk,
+
+```{math}
+:label: csdv_transition
+\theta_t = \theta_{t-1} + v_t,
+\qquad
+v_t \sim N(0,Q).
+```
+
+A prior over how fast coefficients drift plays a mirror-image role in
+{doc}`phillips_priors`: there it is the *government's* prior about a drifting
+Phillips curve that shapes the policy it chooses, whereas here it is the
+*econometrician's* prior in a posterior about drifting reduced-form dynamics.
+
+The companion system is stable when every companion-matrix eigenvalue lies
+strictly inside the unit circle.
+
+For an AR(1), this is $|\rho|<1$; equivalently, the zero of $1-\rho z$ lies
+outside the unit circle.
+
+Cogley and Sargent rule out explosive paths by retaining a path only when the
+companion matrix is stable at every date and using the truncated prior
+
+```{math}
+:label: csdv_stability
+p(\theta^T,Q) \propto I(\theta^T) f(\theta^T \mid Q) f(Q),
+```
+
+where $I(\theta^T)=1$ denotes a stable path.
+
+This restriction encodes the belief that the economy did not in fact follow an
+explosive path.
+
+This restriction also tilts the marginal prior for $Q$ toward values that are
+less likely to generate explosive coefficient paths.
+
+The code below applies the stability restriction to an entire trajectory
+
+```{code-cell} ipython3
+def companion_roots(θ_path):
+ """Return all companion roots along a path with shape (21, T)."""
+ θ_path = np.asarray(θ_path, dtype=float)
+ if θ_path.ndim == 1:
+ θ_path = θ_path[:, None]
+ companions = np.stack(
+ [
+ companion_matrix(θ_path[:, t])[1]
+ for t in range(θ_path.shape[1])
+ ]
+ )
+ return np.linalg.eigvals(companions)
+
+
+def is_stable(θ_path):
+ """Test whether every companion root is strictly inside the unit circle."""
+ return bool(np.max(np.abs(companion_roots(θ_path))) < 1)
+```
+
+The reduced-form innovation covariance changes over time according to
+
+```{math}
+:label: csdv_covariance
+\varepsilon_t = R_t^{1/2}\xi_t,
+\qquad
+\xi_t \sim N(0,I_3),
+\qquad
+R_t = B^{-1} H_t B^{-1\prime},
+```
+
+where
+
+$$
+B =
+\begin{bmatrix}
+1 & 0 & 0 \\
+\beta_{21} & 1 & 0 \\
+\beta_{31} & \beta_{32} & 1
+\end{bmatrix},
+\qquad
+H_t = \operatorname{diag}(h_{1t},h_{2t},h_{3t}).
+$$
+
+The diagonal elements $h_{it}$ let the size of each orthogonalized shock wax and
+wane over time.
+
+The next two functions construct the triangular factor and the reduced-form
+innovation covariance.
+
+```{code-cell} ipython3
+def b_matrix(β):
+ """Construct B from β_21, β_31, and β_32."""
+ matrix = np.eye(n_variables)
+ matrix[1, 0], matrix[2, 0], matrix[2, 1] = np.asarray(β, dtype=float)
+ return matrix
+
+
+def innovation_covariance(h, β):
+ """Construct R_t from one vector of orthogonalized variances."""
+ inverse = np.linalg.inv(b_matrix(β))
+ return inverse @ np.diag(h) @ inverse.T
+```
+
+Each diagonal volatility is a geometric random walk,
+
+```{math}
+:label: csdv_volatility
+\log h_{it} = \log h_{i,t-1} + \sigma_i \eta_{it},
+\qquad
+\eta_{it} \sim N(0,1).
+```
+
+The standardized measurement innovations, coefficient innovations, and
+volatility innovations are mutually independent.
+
+Setting $Q=0$ produces constant coefficients with drifting volatility, while
+holding $H_t$ fixed produces drifting coefficients with constant volatility.
+
+The posterior contains the full paths $\theta^T$ and $H^T$ together with $Q$,
+$\beta$, and $(\sigma_1,\sigma_2,\sigma_3)$.
+
+This posterior has thousands of dimensions, so later sections build a sampler
+that updates one group of parameters at a time, holding the rest fixed.
+
+## The data
+
+We begin with {cite:t}`CogleySargent2005`'s quarterly U.S. dataset ending in 2000Q4.
+
+Inflation is the log difference of the seasonally adjusted CPI for all urban
+consumers, point sampled in the third month of each quarter.
+
+Unemployment is the quarterly average of the seasonally adjusted civilian
+unemployment rate and enters the VAR as $0.01\log[u/(1-u)]$, a logit
+transformation that maps a bounded rate into an unconstrained variable.
+
+The nominal interest rate is the log of one plus the three-month Treasury-bill
+rate, averaged over daily observations in the first month of each quarter and
+expressed as a quarterly fraction.
+
+The data starts in 1948Q2 because its first inflation observation is
+already differenced, so two VAR lags make 1948Q4 the first usable regression
+date.
+
+The following cell performs every transformation and constructs the VAR(2) data
+directly from the series.
+
+```{code-cell} ipython3
+def prepare_data(source, ordering=('i', 'u', 'pi')):
+ """Transform a quarterly table and construct the VAR data."""
+ if isinstance(source, str):
+ table = pd.read_csv(source)
+ else:
+ table = source.copy()
+ variables = {
+ 'i': table['y3'].to_numpy(dtype=float),
+ 'u': 0.01 * np.log(
+ table['ur'].to_numpy(dtype=float)
+ / (1 - table['ur'].to_numpy(dtype=float))
+ ),
+ 'pi': table['dp'].to_numpy(dtype=float),
+ }
+ if sorted(ordering) != ['i', 'pi', 'u']:
+ raise ValueError("ordering must be a permutation of ('i', 'u', 'pi')")
+ raw_y = np.column_stack([variables[name] for name in ordering])
+ raw_dates = table['date'].to_numpy(dtype=float)
+ regressors = np.ones((len(table) - n_lags, n_regressors))
+ for lag in range(1, n_lags + 1):
+ left = 1 + n_variables * (lag - 1)
+ regressors[:, left:left + n_variables] = raw_y[n_lags-lag:-lag]
+ targets = raw_y[n_lags:]
+ dates = raw_dates[n_lags:]
+ n_training = 4 * 11 - n_lags - 1
+ return {
+ 'raw_dates': raw_dates,
+ 'raw_y': raw_y,
+ 'prior_dates': dates[:n_training],
+ 'prior_y': targets[:n_training],
+ 'prior_x': regressors[:n_training],
+ 'dates': dates[n_training:],
+ 'y': targets[n_training:],
+ 'x': regressors[n_training:],
+ }
+
+
+data = prepare_data(data_url)
+
+data_summary = pd.Series(
+ {
+ 'ordering': 'interest, unemployment, inflation',
+ 'prior sample': '1948Q4--1958Q4',
+ 'prior observations': len(data['prior_dates']),
+ 'posterior sample': '1959Q1--2000Q4',
+ 'posterior observations': len(data['dates']),
+ 'VAR lags': n_lags,
+ 'coefficient dimension': n_coefficients,
+ },
+ name='value',
+)
+
+data_summary.to_frame()
+```
+
+The early observations calibrate the prior, while the remaining observations
+form the posterior sample.
+
+Let's view the data in familiar economic units.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Observed $\pi_t$, $u_t$, and $i_t$, 1948Q2--2000Q4
+ name: fig-csdv-historical-data
+---
+dates_raw = data['raw_dates']
+interest = 400 * np.expm1(data['raw_y'][:, 0])
+unemployment = 100 * expit(100 * data['raw_y'][:, 1])
+inflation = 400 * data['raw_y'][:, 2]
+
+fig, axes = plt.subplots(3, 1, figsize=(9, 7), sharex=True)
+axes[0].plot(dates_raw, inflation, lw=2)
+axes[0].set_ylabel('inflation (annual %)')
+axes[1].plot(dates_raw, unemployment, lw=2)
+axes[1].set_ylabel('unemployment (%)')
+axes[2].plot(dates_raw, interest, lw=2)
+axes[2].set_ylabel('interest (annual %)')
+axes[2].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+Inflation and nominal interest rates rise together through the 1970s and peak
+around 1980, unemployment peaks only after disinflation begins, and all three
+series are calmer in the 1990s.
+
+These observations locate the episode but cannot tell whether changed
+systematic dynamics or unusually large shocks produced it, which is the
+distinction the model is built to examine.
+
+## Priors
+
+The priors are independent across parameters and deliberately weak so that, in
+Cogley and Sargent's phrase, "the data are free to speak."
+
+They are calibrated from a time-invariant VAR fitted to the short 1948--1958
+training sample.
+
+The initial coefficient prior is a stable, truncated Gaussian,
+
+$$
+p(\theta_0) \propto I(\theta_0)N(\bar\theta,\bar P),
+$$
+
+where $\bar\theta$ and $\bar P$ come from a constant-coefficient seemingly
+unrelated regression fitted to the 1948--1958 training sample.
+
+Because all three equations have identical regressors, the SUR coefficient
+estimates equal equation-by-equation OLS estimates, which lets us implement the
+calibration compactly.
+
+```{code-cell} ipython3
+def sur_prior(y, x):
+ """Calibrate the Gaussian coefficient prior from a constant VAR."""
+ xx_inverse = np.linalg.inv(x.T @ x)
+ coefficients = xx_inverse @ x.T @ y
+ residuals = y - x @ coefficients
+ residual_covariance = np.cov(residuals, rowvar=False, ddof=1)
+ θ = coefficients.T.reshape(-1)
+ covariance = np.kron(residual_covariance, xx_inverse)
+ return θ, covariance, residual_covariance
+
+
+θ_bar, p_bar, r_bar = sur_prior(data['prior_y'], data['prior_x'])
+
+assert is_stable(θ_bar)
+```
+
+The coefficient-drift covariance has the inverse-Wishart prior
+
+```{math}
+:label: csdv_q_prior
+Q \sim IW_{21}\left(T_0,T_0\bar Q\right),
+\qquad
+T_0 = 22,
+\qquad
+\bar Q = \gamma^2 \bar P,
+\qquad
+\gamma^2 = 3.5\times 10^{-4}.
+```
+
+The convention in {eq}`csdv_q_prior` lists degrees of freedom first and the
+inverse-Wishart scale matrix second.
+
+Because $T_0$ is only one greater than the dimension of $\theta_t$, the prior is
+proper but has no finite mean.
+
+The matrix $\bar Q$ is therefore a conservative scale calibration rather than
+the expectation of $Q$.
+
+This calibration favors slow coefficient drift before the posterior sees the
+main sample.
+
+The remaining priors are
+
+$$
+\begin{aligned}
+\log h_{i0} &\sim N(\log \bar R_{ii},10), \\
+\beta &\sim N(0,10000I_3), \\
+\sigma_i^2 &\sim IG\left(\frac{1}{2},\frac{0.01^2}{2}\right),
+\end{aligned}
+$$
+
+where $\bar R$ is the residual covariance from the training-sample regression.
+
+The next function gathers these hyperparameters.
+
+```{code-cell} ipython3
+def calibrate_prior(model_data, γ_squared=3.5e-4):
+ """Return every calibrated prior for one variable ordering."""
+ θ_mean, θ_covariance, residual_covariance = sur_prior(
+ model_data['prior_y'], model_data['prior_x']
+ )
+ degrees_freedom = n_coefficients + 1
+ q_center = γ_squared * θ_covariance
+ return {
+ 'θ_mean': θ_mean,
+ 'θ_covariance': θ_covariance,
+ 'q_center': q_center,
+ 'q_scale': degrees_freedom * q_center,
+ 'q_degrees_freedom': degrees_freedom,
+ 'log_h_mean': np.log(np.diag(residual_covariance)),
+ 'log_h_variance': 10.0,
+ 'β_mean': np.zeros(3),
+ 'β_variance': 10000.0,
+ 'σ_degrees_freedom': 1.0,
+ 'σ_scale': 0.01**2,
+ 'γ_squared': γ_squared,
+ }
+
+
+prior = calibrate_prior(data)
+q_bar = prior['q_center']
+
+prior_summary = pd.Series(
+ {
+ 'dim(θ)': n_coefficients,
+ 'T0': prior['q_degrees_freedom'],
+ 'γ squared': prior['γ_squared'],
+ 'trace(Q bar)': np.trace(q_bar),
+ 'log-h prior variance': prior['log_h_variance'],
+ 'β prior variance': prior['β_variance'],
+ 'σ-squared IG shape': prior['σ_degrees_freedom'] / 2,
+ 'σ-squared IG scale': prior['σ_scale'] / 2,
+ },
+ name='value',
+)
+
+prior_summary.to_frame()
+```
+
+(csdv-sampler)=
+## A Metropolis-within-Gibbs sampler
+
+We simulate the posterior by cycling through five parameter blocks used by
+{cite:t}`CogleySargent2005`.
+
+One pass through all five blocks is called a sweep, and the sampler runs many
+sweeps to build up the posterior sample.
+
+Cogley and Sargent simulate the unrestricted posterior and then discard a
+complete MCMC realization whenever its coefficient path is explosive.
+
+But in this case only stable sweeps contribute realizations
+to the retained restricted-posterior sample.
+
+We instead impose stability inside the coefficient-path block with the
+elliptical slice sampler of {cite:t}`MurrayAdamsMacKay2010`.
+
+This changes the transition kernel, but not its target posterior.
+
+1. An auxiliary Gaussian coefficient path is drawn by a Carter--Kohn
+ forward-filter/backward-sample step and used to update the stable path by
+ elliptical slice sampling.
+
+2. The drift covariance $Q$ is drawn from an inverse-Wishart distribution
+ conditional on the coefficient innovations.
+
+3. The volatility innovation variances $\sigma_i^2$ are drawn from inverse-gamma
+ distributions conditional on the volatility increments.
+
+4. The covariance parameters $\beta$ are drawn from two transformed Gaussian
+ regressions among the VAR residuals.
+
+5. The volatility paths $H^T$ are drawn one date at a time by a
+ Jacquier--Polson--Rossi Metropolis step.
+
+This ordering matters: every block in a sweep is paired with the values on which
+it was actually conditioned.
+
+### Coefficient path
+
+Conditional on $R^T$ and $Q$, a forward Kalman filter followed by the
+Carter--Kohn backward simulator draws the entire coefficient path
+{cite}`CarterKohn1994`.
+
+The forward pass is a Kalman filter, while the backward pass samples
+$\theta_T,\theta_{T-1},\ldots,\theta_0$ in reverse with each state conditioned
+on the draw that follows it.
+
+Sampling $\theta_0$ is essential.
+
+It supplies all $T$ random-walk increments to the conjugate update for $Q$.
+
+```{code-cell} ipython3
+def covariance_root(matrix):
+ """Return a numerically stable lower covariance factor."""
+ matrix = 0.5 * (matrix + matrix.T)
+ scale = max(1.0, np.max(np.abs(np.diag(matrix))))
+ return np.linalg.cholesky(matrix + 1e-12 * scale * np.eye(len(matrix)))
+
+
+def draw_coefficient_path(
+ y, x, q, h, β, prior, rng, return_mean=False
+):
+ """Draw θ_0,...,θ_T and optionally return its smoothing mean."""
+ periods = len(y)
+ filtered_mean = np.empty((periods + 1, n_coefficients))
+ filtered_covariance = np.empty(
+ (periods + 1, n_coefficients, n_coefficients)
+ )
+ predicted_covariance = np.empty_like(filtered_covariance)
+ filtered_mean[0] = prior['θ_mean']
+ filtered_covariance[0] = prior['θ_covariance']
+ predicted_covariance[0] = prior['θ_covariance']
+ for t in range(1, periods + 1):
+ observation = design_matrix(x[t - 1])
+ prediction_covariance = filtered_covariance[t - 1] + q
+ r_t = innovation_covariance(h[t], β)
+ forecast_covariance = (
+ observation @ prediction_covariance @ observation.T + r_t
+ )
+ gain = linalg.solve(
+ forecast_covariance,
+ (prediction_covariance @ observation.T).T,
+ assume_a='pos',
+ ).T
+ mean = filtered_mean[t - 1]
+ mean = mean + gain @ (y[t - 1] - observation @ mean)
+ covariance = (
+ prediction_covariance
+ - gain @ observation @ prediction_covariance
+ )
+ covariance = 0.5 * (covariance + covariance.T)
+ filtered_mean[t] = mean
+ filtered_covariance[t] = covariance
+ predicted_covariance[t] = prediction_covariance
+ path = np.empty((n_coefficients, periods + 1))
+ if not return_mean:
+ path[:, -1] = (
+ filtered_mean[-1]
+ + covariance_root(filtered_covariance[-1])
+ @ rng.standard_normal(n_coefficients)
+ )
+ for t in range(periods - 1, -1, -1):
+ smoother = linalg.solve(
+ predicted_covariance[t + 1],
+ filtered_covariance[t].T,
+ assume_a='pos',
+ ).T
+ mean = filtered_mean[t] + smoother @ (
+ path[:, t + 1] - filtered_mean[t]
+ )
+ covariance = (
+ filtered_covariance[t]
+ - smoother @ predicted_covariance[t + 1] @ smoother.T
+ )
+ path[:, t] = mean + covariance_root(covariance) @ (
+ rng.standard_normal(n_coefficients)
+ )
+ return path
+
+ smoothed_mean = np.empty_like(path)
+ centered_draw = np.empty_like(path)
+ smoothed_mean[:, -1] = filtered_mean[-1]
+ centered_draw[:, -1] = covariance_root(filtered_covariance[-1]) @ (
+ rng.standard_normal(n_coefficients)
+ )
+ for t in range(periods - 1, -1, -1):
+ smoother = linalg.solve(
+ predicted_covariance[t + 1],
+ filtered_covariance[t].T,
+ assume_a='pos',
+ ).T
+ smoothed_mean[:, t] = filtered_mean[t] + smoother @ (
+ smoothed_mean[:, t + 1] - filtered_mean[t]
+ )
+ covariance = (
+ filtered_covariance[t]
+ - smoother @ predicted_covariance[t + 1] @ smoother.T
+ )
+ centered_draw[:, t] = (
+ smoother @ centered_draw[:, t + 1]
+ + covariance_root(covariance) @ rng.standard_normal(n_coefficients)
+ )
+ return smoothed_mean + centered_draw, smoothed_mean
+
+
+def draw_stable_coefficient_path(
+ current, y, x, q, h, β, prior, rng, max_contractions=100
+):
+ """Elliptical-slice update of the stability-truncated Gaussian path."""
+ if not is_stable(current):
+ raise ValueError('the elliptical-slice update needs a stable path')
+ gaussian_draw, mean = draw_coefficient_path(
+ y, x, q, h, β, prior, rng, return_mean=True
+ )
+ current_centered = current - mean
+ innovation = gaussian_draw - mean
+ angle = rng.uniform(0, 2 * np.pi)
+ lower = angle - 2 * np.pi
+ upper = angle
+ for contractions in range(max_contractions + 1):
+ proposal = (
+ mean
+ + current_centered * np.cos(angle)
+ + innovation * np.sin(angle)
+ )
+ if is_stable(proposal):
+ return proposal, contractions
+ if angle < 0:
+ lower = angle
+ else:
+ upper = angle
+ angle = rng.uniform(lower, upper)
+ raise RuntimeError('elliptical-slice stability bracket did not contract')
+```
+
+### Drift covariance
+
+Conditional on the coefficient increments, $Q$ has an inverse-Wishart full
+conditional.
+
+The retained simulation draws $Q$ from a scale matrix formed by its prior scale
+and all $T$ squared increments from $\theta_0$ through $\theta_T$.
+
+```{code-cell} ipython3
+def draw_q(θ_path, prior, rng):
+ """Draw Q conditional on the sampled coefficient path."""
+ increments = np.diff(θ_path, axis=1)
+ scale = prior['q_scale'] + increments @ increments.T
+ degrees_freedom = prior['q_degrees_freedom'] + increments.shape[1]
+ return invwishart.rvs(df=degrees_freedom, scale=scale, random_state=rng)
+```
+
+### Volatility parameters and paths
+
+Conditional on the volatility increments, each $\sigma_i^2$ has an inverse-gamma
+full conditional.
+
+Conditional on the VAR residuals and $H^T$, the free elements of $B$ are drawn
+from two Gaussian regressions.
+
+Conditional on the orthogonalized residuals, each volatility state is updated
+with the single-site Metropolis step of {cite:t}`Jacquier1994`.
+
+The random-walk neighbors determine the Gaussian proposal for a log volatility,
+while the corresponding orthogonalized residual determines whether that proposal
+is accepted.
+
+The following code implements the conditional updates, including the different
+endpoint proposals.
+
+```{code-cell} ipython3
+def var_residuals(y, x, θ_path):
+ """Return residuals with shape (T, 3)."""
+ if θ_path.shape[1] == len(y) + 1:
+ θ_path = θ_path[:, 1:]
+ if θ_path.shape[1] != len(y):
+ raise ValueError('θ_path must contain T or T + 1 states')
+ coefficients = θ_path.T.reshape(len(y), n_variables, n_regressors)
+ fitted = np.einsum('tk,tnk->tn', x, coefficients)
+ return y - fitted
+
+
+def draw_σ(h, prior, rng):
+ """Draw the three log-volatility innovation standard deviations."""
+ increments = np.diff(np.log(h), axis=0)
+ shape = (prior['σ_degrees_freedom'] + increments.shape[0]) / 2
+ scales = (prior['σ_scale'] + np.sum(increments**2, axis=0)) / 2
+ σ_squared = scales / rng.gamma(shape, 1.0, size=n_variables)
+ return np.sqrt(σ_squared)
+
+
+def draw_β(residuals, h, prior, rng):
+ """Draw the free elements of B from transformed Gaussian regressions."""
+ β = np.empty(3)
+ offset = 0
+ for equation in range(1, n_variables):
+ standardized = residuals / np.sqrt(h[1:, equation])[:, None]
+ dependent = standardized[:, equation]
+ regressors = -standardized[:, :equation]
+ prior_precision = np.eye(equation) / prior['β_variance']
+ covariance = np.linalg.inv(prior_precision + regressors.T @ regressors)
+ prior_slice = prior['β_mean'][offset:offset + equation]
+ mean = covariance @ (
+ prior_precision @ prior_slice + regressors.T @ dependent
+ )
+ β[offset:offset + equation] = (
+ mean + covariance_root(covariance) @ rng.standard_normal(equation)
+ )
+ offset += equation
+ return β
+
+
+def accept_volatility(proposal, current, residual, rng):
+ """Apply the Jacquier--Polson--Rossi likelihood acceptance step."""
+ log_ratio = (
+ -0.5 * np.log(proposal)
+ - residual**2 / (2 * proposal)
+ + 0.5 * np.log(current)
+ + residual**2 / (2 * current)
+ )
+ return proposal if np.log(rng.random()) <= min(0.0, log_ratio) else current
+
+
+def draw_volatility_path(h, residuals, β, σ, prior, rng):
+ """Update all stochastic-volatility states one date at a time."""
+ periods = len(residuals)
+ orthogonalized = (b_matrix(β) @ residuals.T).T
+ updated = np.empty_like(h)
+ for equation in range(n_variables):
+ variance = σ[equation]**2
+ initial_variance = (
+ prior['log_h_variance'] * variance
+ / (variance + prior['log_h_variance'])
+ )
+ initial_mean = initial_variance * (
+ prior['log_h_mean'][equation] / prior['log_h_variance']
+ + np.log(h[1, equation]) / variance
+ )
+ updated[0, equation] = np.exp(
+ initial_mean + np.sqrt(initial_variance) * rng.standard_normal()
+ )
+ for t in range(1, periods):
+ mean = 0.5 * (
+ np.log(updated[t - 1, equation]) + np.log(h[t + 1, equation])
+ )
+ proposal = np.exp(
+ mean + np.sqrt(variance / 2) * rng.standard_normal()
+ )
+ updated[t, equation] = accept_volatility(
+ proposal,
+ h[t, equation],
+ orthogonalized[t - 1, equation],
+ rng,
+ )
+ proposal = np.exp(
+ np.log(updated[-2, equation])
+ + σ[equation] * rng.standard_normal()
+ )
+ updated[-1, equation] = accept_volatility(
+ proposal,
+ h[-1, equation],
+ orthogonalized[-1, equation],
+ rng,
+ )
+ return updated
+```
+
+### Complete sampler
+
+The restricted posterior assigns zero density to coefficient paths that are
+explosive at any date, including $\theta_0$.
+
+Stack the whole coefficient path $\theta_0,\ldots,\theta_T$ into $z$, and
+collect the other parameter blocks in $\lambda=(Q,H^T,\beta,\sigma)$.
+
+Conditional on $\lambda$ and the data $Y^T$, the unrestricted Carter--Kohn
+distribution is Gaussian with smoothing mean $m$ and covariance $C$,
+
+$$
+z\mid \lambda,Y^T \sim N(m,C).
+$$
+
+Let $\mathcal A$ be the set of paths whose companion roots are strictly inside
+the unit circle at every date, so that the restricted full conditional is this
+Gaussian truncated to $\mathcal A$,
+
+$$
+\pi_{\mathcal A}(z\mid\lambda,Y^T)
+= \frac{N(z;m,C)\,\mathbb{1}_{\mathcal A}(z)}
+ {\mathbb{P}\{z\in\mathcal A\mid\lambda,Y^T\}}.
+$$
+
+That normalizing probability is difficult to compute, but the elliptical
+transition never evaluates it.
+
+Cogley and Sargent instead simulate the unrestricted joint posterior and keep
+a realization only when its whole path is stable, which is valid because
+conditioning the unrestricted posterior on $z\in\mathcal A$ reproduces the
+restricted posterior exactly.
+
+If $a$ denotes the unrestricted probability that a path is stable, this
+rejection scheme needs roughly $1/a$ complete sweeps, including the four other
+parameter updates, for every retained draw.
+
+The elliptical slice sampler avoids that waste by moving within the stable
+region instead of restarting from an arbitrary draw.
+
+It starts from the current stable path $z^{(c)}$ and a fresh Carter--Kohn
+draw $\widetilde z\sim N(m,C)$, whose centered version
+$\nu=\widetilde z-m$ traces an ellipse together with $z^{(c)}$,
+
+$$
+z(\phi)
+=m+(z^{(c)}-m)\cos\phi+\nu\sin\phi,
+\qquad 0\leq\phi<2\pi,
+$$
+
+so that $z(0)=z^{(c)}$ and $z(\pi/2)=\widetilde z$.
+
+The algorithm draws an angle uniformly from the full circle and, whenever
+$z(\phi)$ is explosive, shrinks the bracket to the side containing the
+known-stable angle $\phi=0$ before drawing again.
+
+Because companion roots vary continuously with the coefficients, a nonzero
+interval around $\phi=0$ is always stable, so this bracket search always
+terminates.
+
+Each rejected angle costs only a linear combination and a companion-root
+check, far cheaper than another Kalman filter and backward simulation.
+
+The transition is valid because rotating the pair $(z^{(c)}-m,\nu)$ by any
+angle leaves their joint Gaussian density unchanged, so the search moves
+along a fixed orbit on which every point is equally likely under the
+unrestricted density.
+
+The stability indicator plays the role of the likelihood in a standard
+elliptical slice update, and since it equals one at the current point, every
+accepted angle is automatically a stable one and no separate slice-height
+draw is needed.
+
+Marginalizing out the auxiliary path shows that this transition leaves the
+truncated Gaussian $N(z;m,C)\mathbb{1}_{\mathcal A}(z)$ invariant, exactly the
+property a valid transition kernel needs.
+
+The other four blocks require no such adjustment.
+
+Conditional on a stable $z$, the stability indicator is constant in
+$Q,H^T,\beta$, and $\sigma$, so it cancels from each of their full
+conditionals and leaves the same updates as the unrestricted sampler.
+
+$Q$'s conditional is unchanged in form, but its marginal posterior still tilts
+toward less explosive drift because every draw of $Q$ is conditioned on a
+stable coefficient path.
+
+Together, the elliptical transition for $z$ and the unchanged updates for
+$Q,H^T,\beta$, and $\sigma$ target the same stability-restricted posterior as
+Cogley and Sargent's original rejection sampler, at a fraction of the cost.
+
+The next function composes these five blocks and includes a
+stochastic-volatility warm-up.
+
+```{code-cell} ipython3
+def initial_volatilities(y, prior):
+ """Construct the sampler's initial volatility path."""
+ changes = np.diff(y, axis=0)
+ centered = changes - changes.mean(axis=0)
+ log_h = np.empty((len(y) + 1, n_variables))
+ log_h[:2] = prior['log_h_mean']
+ log_h[2:] = np.log(np.maximum(centered**2, np.finfo(float).tiny))
+ return np.exp(log_h)
+
+
+def run_sampler(
+ y,
+ x,
+ prior,
+ n_sweeps=1_000,
+ burn=500,
+ thin=1,
+ seed=42,
+ warmup=200,
+ max_contractions=100,
+ stable=True,
+ fixed_q=None,
+ retain=('S0D', 'SD', 'QD', 'HD', 'CD', 'VD', 'stable_draw'),
+ progress_every=0,
+):
+ """Run a Gibbs sampler for the unrestricted or stable posterior.
+
+ For the stable posterior, an elliptical-slice transition updates the FFBS
+ path inside its stability-truncated Gaussian full conditional. Passing
+ ``fixed_q`` holds Q at that value (for example a matrix of zeros) instead of
+ drawing it, which nests the constant-coefficient model.
+ """
+ if not (0 <= burn < n_sweeps and thin >= 1):
+ raise ValueError('require 0 <= burn < n_sweeps and thin >= 1')
+ if (n_sweeps - burn) % thin:
+ raise ValueError('(n_sweeps - burn) must be divisible by thin')
+ valid_retain = {'S0D', 'SD', 'QD', 'HD', 'CD', 'VD', 'stable_draw'}
+ unknown = set(retain) - valid_retain
+ if unknown:
+ raise ValueError(f'unknown retained arrays: {sorted(unknown)}')
+
+ started = time.perf_counter()
+ rng = np.random.default_rng(seed)
+ h = initial_volatilities(y, prior)
+ β = prior['β_mean'].copy()
+ warm_θ = np.repeat(prior['θ_mean'][:, None], len(y), axis=1)
+ warm_residuals = var_residuals(y, x, warm_θ)
+ for _ in range(warmup):
+ σ = draw_σ(h, prior, rng)
+ β = draw_β(warm_residuals, h, prior, rng)
+ h = draw_volatility_path(
+ h, warm_residuals, β, σ, prior, rng
+ )
+
+ q = (
+ prior['q_center'].copy()
+ if fixed_q is None
+ else np.array(fixed_q, dtype=float)
+ )
+ θ = np.repeat(
+ prior['θ_mean'][:, None], len(y) + 1, axis=1
+ )
+ if stable and not is_stable(θ):
+ raise ValueError('the prior mean does not provide a stable start')
+ slice_contractions = 0
+ maximum_slice_contractions = 0
+
+ retained = {name: [] for name in retain}
+ saved_stability = []
+ for sweep in range(1, n_sweeps + 1):
+ if stable:
+ θ, contractions = draw_stable_coefficient_path(
+ θ,
+ y,
+ x,
+ q,
+ h,
+ β,
+ prior,
+ rng,
+ max_contractions=max_contractions,
+ )
+ else:
+ θ = draw_coefficient_path(y, x, q, h, β, prior, rng)
+ contractions = 0
+ slice_contractions += contractions
+ maximum_slice_contractions = max(
+ maximum_slice_contractions, contractions
+ )
+ if fixed_q is None:
+ q = draw_q(θ, prior, rng)
+ residuals = var_residuals(y, x, θ)
+ σ = draw_σ(h, prior, rng)
+ β = draw_β(residuals, h, prior, rng)
+ h = draw_volatility_path(
+ h, residuals, β, σ, prior, rng
+ )
+
+ if sweep > burn and (sweep - burn) % thin == 0:
+ path_is_stable = is_stable(θ)
+ saved_stability.append(path_is_stable)
+ values = {
+ 'S0D': θ[:, 0],
+ 'SD': θ[:, 1:],
+ 'QD': q,
+ 'HD': h,
+ 'CD': β,
+ 'VD': σ,
+ 'stable_draw': path_is_stable,
+ }
+ for name in retained:
+ retained[name].append(np.asarray(values[name]).copy())
+
+ if progress_every and sweep % progress_every == 0:
+ elapsed = time.perf_counter() - started
+ print(
+ f'{sweep:,}/{n_sweeps:,} sweeps; '
+ f'{slice_contractions:,} slice contractions; '
+ f'{elapsed / 60:.1f} minutes',
+ flush=True,
+ )
+
+ stack_axis = {
+ 'S0D': 1,
+ 'SD': 2,
+ 'QD': 2,
+ 'HD': 2,
+ 'CD': 1,
+ 'VD': 1,
+ 'stable_draw': 0,
+ }
+ result = {
+ name: np.stack(values, axis=stack_axis[name])
+ for name, values in retained.items()
+ }
+ result['diagnostics'] = {
+ 'sampler_version': 3,
+ 'seed': int(seed),
+ 'stable_restriction': bool(stable),
+ 'n_sweeps': int(n_sweeps),
+ 'burn': int(burn),
+ 'thin': int(thin),
+ 'warmup': int(warmup),
+ 'retained_draws': int((n_sweeps - burn) // thin),
+ 'slice_contractions': int(slice_contractions),
+ 'mean_slice_contractions': float(slice_contractions / n_sweeps),
+ 'maximum_slice_contractions': int(maximum_slice_contractions),
+ 'retained_stability_rate': float(np.mean(saved_stability)),
+ 'elapsed_seconds': float(time.perf_counter() - started),
+ }
+ return result
+```
+
+The executable version below uses 1,000 sweeps, discards the first 500, and
+retains the remaining 500.
+
+It uses the complete historical sample and the ordering $(i,u,\pi)$.
+
+Instead of running a large MCMC experiment, we intentionally keep the sampler
+run small so that it finishes in a reasonable time for a lecture.
+
+This short run illustrates the method but not a numerical replication.
+
+However the main qualitative features of the posterior are close to those reported by Cogley and Sargent.
+
+```{code-cell} ipython3
+posterior = run_sampler(
+ data['y'],
+ data['x'],
+ prior,
+ n_sweeps=1_000,
+ burn=500,
+ thin=1,
+ seed=42,
+ warmup=200,
+ stable=True,
+ progress_every=0,
+)
+
+def validate_posterior_arrays(result, periods):
+ """Check posterior shapes, finiteness, positivity, and stability."""
+ draws = result['diagnostics']['retained_draws']
+ expected = {
+ 'S0D': (n_coefficients, draws),
+ 'SD': (n_coefficients, periods, draws),
+ 'QD': (n_coefficients, n_coefficients, draws),
+ 'HD': (periods + 1, n_variables, draws),
+ 'CD': (3, draws),
+ 'VD': (3, draws),
+ 'stable_draw': (draws,),
+ }
+ assert {name: result[name].shape for name in expected} == expected
+ assert all(np.all(np.isfinite(result[name])) for name in expected)
+ assert np.all(result['HD'] > 0)
+ assert np.all(result['VD'] > 0)
+ assert np.all(result['stable_draw'])
+ return expected
+
+
+expected_shapes = validate_posterior_arrays(posterior, len(data['dates']))
+```
+
+(csdv-results)=
+## What the data say
+
+We summarize the posterior by its mean coefficient path $\mathbb{E}(\theta_t\mid T)$ and
+mean covariance path $\mathbb{E}(R_t\mid T)$, and then interpret them
+in the context of the question we asked.
+
+### The rate and structure of drift
+
+The trace of $Q$ measures the total rate of coefficient drift, with
+$\operatorname{tr}(Q)=0$ corresponding to constant coefficients.
+
+The histogram shows the retained $Q$ draws and the prior scale.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior $\operatorname{tr}(Q)$ and prior $\operatorname{tr}(\bar Q)$
+ name: fig-csdv-drift-rate
+---
+trace_q = np.trace(posterior['QD'], axis1=0, axis2=1)
+fig, ax = plt.subplots()
+ax.hist(trace_q, bins=30, histtype='step', lw=2)
+ax.axvline(np.trace(q_bar), color='C1', lw=2,
+ label=r'prior $\mathrm{tr}(\bar Q)$')
+ax.set_xlabel(r'$\mathrm{tr}(Q)$')
+ax.set_ylabel('frequency')
+ax.legend()
+plt.show()
+```
+
+The posterior drift rate lies well above the conservative prior calibration,
+indicating more coefficient variation than that calibration anticipated.
+
+This is not a formal comparison with a fixed-coefficient model because the
+continuous prior assigns no point mass to $Q=0$.
+
+Within the fitted TVP-VAR, this variation is attributed to changing systematic
+relationships; it does not identify policy as the cause.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior mean VAR coefficients $E(\theta_t\mid T)$
+ name: fig-csdv-coefficient-paths
+---
+θ_mean = posterior['SD'].mean(axis=2)
+mean_path_root_modulus = np.max(
+ np.abs(companion_roots(θ_mean)), axis=1
+)
+assert np.all(mean_path_root_modulus < 1)
+coefficient_labels = (
+ 'constant',
+ r'$i_{t-1}$',
+ r'$u_{t-1}$',
+ r'$\pi_{t-1}$',
+ r'$i_{t-2}$',
+ r'$u_{t-2}$',
+ r'$\pi_{t-2}$',
+)
+equation_labels = (
+ 'interest equation',
+ 'unemployment equation',
+ 'inflation equation',
+)
+
+
+def plot_equation_coefficients(axes, dates, θ_path):
+ """Plot seven labeled coefficients for each VAR equation."""
+ first_lines = None
+ for equation, (ax, label) in enumerate(zip(axes, equation_labels)):
+ start = equation * len(coefficient_labels)
+ stop = start + len(coefficient_labels)
+ lines = ax.plot(dates, θ_path[start:stop].T, lw=2)
+ for line, coefficient_label in zip(lines, coefficient_labels):
+ line.set_label(coefficient_label)
+ if first_lines is None:
+ first_lines = lines
+ ax.axhline(0, color='0.65', lw=1)
+ ax.set_xlabel('year')
+ ax.set_ylabel(f'{label} coefficient')
+ return first_lines
+
+
+fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharex=True)
+coefficient_lines = plot_equation_coefficients(
+ axes,
+ data['dates'],
+ θ_mean,
+)
+fig.legend(
+ coefficient_lines,
+ coefficient_labels,
+ loc='lower center',
+ ncol=4,
+)
+plt.tight_layout(rect=(0, 0.17, 1, 1))
+plt.show()
+```
+
+The unemployment equation is comparatively stable, whereas several inflation
+equation coefficients move strongly through the 1970s and turn around near
+1980.
+
+The drift is therefore concentrated in how inflation propagates rather than
+spread evenly across the VAR, and individual lag coefficients should not be
+given structural interpretations because paired lags can offset one another.
+
+We summarize drift for the ordering $(i,u,\pi)$, which has the smallest stable
+posterior mean $\operatorname{tr}(Q)$.
+
+```{code-cell} ipython3
+trace_q = np.trace(posterior['QD'], axis1=0, axis2=1)
+q_mean = posterior['QD'].mean(axis=2)
+drift_summary = {
+ r'\text{Posterior mean } \operatorname{tr}(Q)': np.trace(q_mean),
+ r'\text{Posterior mean largest eigenvalue}': (
+ np.linalg.eigvalsh(q_mean)[-1]
+ ),
+ r'\text{Prior } \operatorname{tr}(\bar Q)': np.trace(q_bar),
+}
+drift_rows = ' \\\\\n'.join(
+ f'{label} & {value:.4f}' for label, value in drift_summary.items()
+)
+display(Math(rf'''
+\begin{{array}}{{lr}}
+\text{{Quantity}} & \text{{Estimate}} \\
+\hline
+{drift_rows}
+\end{{array}}
+'''))
+```
+
+Cogley and Sargent estimated every ordering.
+
+Their posterior means show that
+the ordering changes magnitudes but does not remove drift:
+
+| Ordering | Stable $\operatorname{tr}(Q)$ | Stable $\max(\lambda)$ | Unrestricted $\operatorname{tr}(Q)$ | Unrestricted $\max(\lambda)$ |
+|---|---:|---:|---:|---:|
+| $(i,\pi,u)$ | 0.055 | 0.025 | 0.056 | 0.027 |
+| $(i,u,\pi)$ | 0.047 | 0.023 | 0.059 | 0.031 |
+| $(\pi,i,u)$ | 0.064 | 0.031 | 0.082 | 0.044 |
+| $(\pi,u,i)$ | 0.062 | 0.031 | 0.088 | 0.051 |
+| $(u,i,\pi)$ | 0.057 | 0.026 | 0.051 | 0.028 |
+| $(u,\pi,i)$ | 0.055 | 0.024 | 0.072 | 0.035 |
+
+For the minimum-$Q$ ordering, removing the stability restriction raises the
+posterior mean drift rate, as the table above shows.
+
+The analysis that follows adopts the $(i,u,\pi)$ ordering, which places the
+nominal interest rate first and inflation last.
+
+Diagonalizing the posterior mean of $Q$ reveals that the drift is low
+dimensional.
+
+The following eigendecomposition summarizes the posterior mean of $Q$.
+
+```{code-cell} ipython3
+q_mean = posterior['QD'].mean(axis=2)
+q_eigenvalues = np.linalg.eigvalsh(q_mean)[::-1]
+q_cumulative = np.cumsum(q_eigenvalues) / q_eigenvalues.sum()
+
+drift_structure = pd.DataFrame(
+ {
+ 'eigenvalue': q_eigenvalues[:3],
+ 'cumulative share': q_cumulative[:3],
+ },
+ index=pd.Index(range(1, 4), name='principal component'),
+)
+drift_structure.round(4)
+```
+
+The first three principal components account for the large majority of total
+coefficient drift even though the VAR contains 21 coefficients.
+
+### The evolution of volatility
+
+We first ask how the *size* of the shocks changed.
+
+Equation {eq}`csdv_covariance` can be averaged over draws without constructing a
+four-dimensional covariance array.
+
+```{code-cell} ipython3
+def mean_innovation_covariance(h_draws, β_draws):
+ """Compute E(R_t | T) with working memory proportional to T times D."""
+ n_draws = h_draws.shape[2]
+ matrices = np.broadcast_to(np.eye(3), (n_draws, 3, 3)).copy()
+ matrices[:, 1, 0] = β_draws[0]
+ matrices[:, 2, 0] = β_draws[1]
+ matrices[:, 2, 1] = β_draws[2]
+ inverses = np.linalg.solve(
+ matrices,
+ np.broadcast_to(np.eye(3), matrices.shape),
+ )
+ h = h_draws[1:]
+ mean = np.empty((h.shape[0], 3, 3))
+ for row in range(3):
+ for column in range(row + 1):
+ value = np.zeros(h.shape[0])
+ for shock in range(3):
+ weights = inverses[:, row, shock] * inverses[:, column, shock]
+ value += h[:, shock, :] @ weights
+ mean[:, row, column] = value / n_draws
+ mean[:, column, row] = mean[:, row, column]
+ return mean
+
+
+r_mean = mean_innovation_covariance(posterior['HD'], posterior['CD'])
+```
+
+The next plot shows the innovation standard deviations and correlations.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Standard deviations and correlations implied by $E(R_t\mid T)$
+ name: fig-csdv-volatility-correlation
+---
+variances = ((0, 'Nominal interest'), (2, 'Inflation'), (1, 'Unemployment'))
+correlations = (
+ (0, 1, 'Interest--unemployment'),
+ (0, 2, 'Interest--inflation'),
+ (2, 1, 'Inflation--unemployment'),
+)
+fig, axes = plt.subplots(3, 2, figsize=(9, 8), sharex=True)
+for row, (index, label) in enumerate(variances):
+ axes[row, 0].plot(
+ data['dates'], 10000 * np.sqrt(r_mean[:, index, index]), lw=2
+ )
+ axes[row, 0].set_title(label)
+for row, (left, right, label) in enumerate(correlations):
+ scale = np.sqrt(r_mean[:, left, left] * r_mean[:, right, right])
+ axes[row, 1].plot(data['dates'], r_mean[:, left, right] / scale, lw=2)
+ axes[row, 1].set_title(label)
+axes[1, 0].set_ylabel(
+ r'innovation standard deviation $\times 10^4$'
+)
+axes[1, 1].set_ylabel('correlation')
+axes[-1, 0].set_xlabel('year')
+axes[-1, 1].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+The interest-rate and inflation innovation standard deviations peak sharply
+around 1980, whereas unemployment innovation volatility declines more gradually
+toward the end of the sample.
+
+All three innovation correlations also move most abruptly around 1980, so both
+the size and the joint composition of reduced-form shocks changed.
+
+These movements give the changing-shocks, or bad-luck, explanation an important
+role, although the interest-rate innovation is not itself a structural
+monetary-policy shock.
+
+The log determinant of the posterior mean covariance matrix summarizes the
+generalized one-step innovation variance {cite}`Whittle1953`.
+
+The following transformation summarizes generalized innovation variance.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Generalized innovation variance $\log |E(R_t\mid T)|$
+ name: fig-csdv-total-variance
+---
+sign, logdet_r = np.linalg.slogdet(r_mean)
+assert np.all(sign > 0)
+fig, ax = plt.subplots()
+ax.plot(data['dates'], logdet_r, lw=2)
+ax.set_xlabel('year')
+ax.set_ylabel(r'$\log |E(R_t\mid T)|$')
+plt.show()
+```
+
+Because a less negative log determinant means greater joint innovation
+variance, the two-step rise to an exceptional 1981 peak and the long subsequent
+decline mark a large shock episode followed by the Great Moderation documented
+by {cite:t}`KimNelson1999` and {cite:t}`McConnellPerezQuiros2000`.
+
+This is the clearest aggregate evidence for changing luck, but it cannot explain
+the coefficient-based changes in inflation dynamics examined next.
+
+### Core inflation and the natural rate
+
+To study the systematic dynamics, write the VAR at date $t$ in companion form as
+
+$$
+z_t = \mu_{t\mid T} + A_{t\mid T}z_{t-1} + e_t.
+$$
+
+At date $t$, the local mean $m_t$ is the fixed point to which the companion-form
+VAR would converge if its coefficients remained fixed at their posterior means
+and future innovations were zero.
+
+```{math}
+:label: csdv_local_means
+m_t = (I-A_{t\mid T})^{-1}\mu_{t\mid T},
+\qquad
+\bar\pi_t = 4s_\pi m_t,
+\qquad
+\bar u_t = \frac{\exp(100s_u m_t)}{1+\exp(100s_u m_t)}.
+```
+
+Here $s_\pi$ and $s_u$ select inflation and unemployment from $m_t$, the factor
+four annualizes inflation, and the inverse-logit transformation returns
+unemployment to its observed rate.
+
+These are date-specific steady states of frozen local systems rather than
+unconditional means of the globally drifting process.
+
+Core inflation is the long-horizon inflation forecast implied by freezing the
+date-$t$ coefficients, while the natural rate is the corresponding long-run
+unemployment anchor rather than a natural interest rate or a forecast of next
+quarter's unemployment.
+
+Freezing the current coefficients and projecting forward is exactly the
+*anticipated-utility* device used by the learning governments of
+{doc}`phillips_learning` and {doc}`phillips_escaping_nash`, who act as if their
+current beliefs will never be revised.
+
+The following implementation annualizes core inflation and reverses the
+archive's unemployment transformation.
+
+```{code-cell} ipython3
+def local_means(θ_path):
+ """Compute local core inflation and the natural unemployment rate."""
+ θ_path = np.asarray(θ_path, dtype=float)
+ if θ_path.ndim == 1:
+ θ_path = θ_path[:, None]
+ core = np.empty(θ_path.shape[1])
+ natural = np.empty(θ_path.shape[1])
+ for t in range(θ_path.shape[1]):
+ intercept, companion = companion_matrix(θ_path[:, t])
+ mean = np.linalg.solve(np.eye(6) - companion, intercept)
+ core[t] = 4 * mean[2]
+ natural[t] = expit(100 * mean[1])
+ return core, natural
+
+
+core_inflation, natural_rate = local_means(θ_mean)
+```
+
+We plot the fourth-quarter observation from each year.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Local means $\bar\pi_t$ and $\bar u_t$
+ name: fig-csdv-local-means
+---
+def annual_indices(dates, start=1960):
+ """Return the final observation in each year from a start date."""
+ year_values = np.floor(dates + 1e-8).astype(int)
+ return np.array(
+ [
+ np.flatnonzero(year_values == year)[-1]
+ for year in np.unique(year_values)
+ if year >= start
+ ]
+ )
+
+
+years = np.floor(data['dates'] + 1e-8).astype(int)
+annual = annual_indices(data['dates'])
+
+fig, ax = plt.subplots()
+ax.plot(data['dates'][annual], 100 * core_inflation[annual], 'o-', lw=2,
+ markersize=3, label='core inflation')
+ax.plot(data['dates'][annual], 100 * natural_rate[annual], '+-', lw=2,
+ markersize=5, label='natural rate')
+ax.set_xlabel('year')
+ax.set_ylabel('percent')
+ax.legend()
+plt.show()
+```
+
+```{code-cell} ipython3
+core_summary = pd.Series(
+ {
+ 'early-1960s mean core inflation (%)': (
+ 100 * core_inflation[(years >= 1960) & (years <= 1964)].mean()
+ ),
+ 'peak core inflation (%)': 100 * core_inflation[annual].max(),
+ '1985--2000 mean core inflation (%)': (
+ 100 * core_inflation[(years >= 1985) & (years <= 2000)].mean()
+ ),
+ },
+ name='estimate',
+)
+core_summary.to_frame().round(2)
+```
+
+Core inflation climbs from a low level in the early 1960s to a high peak near
+1980 and then falls back, while the natural unemployment rate rises more
+smoothly from about 5 to 6.5 percent before returning toward 4 percent.
+
+Because both lines are determined by the fitted coefficients rather than by
+realized shocks, their persistent shifts point to a changing systematic
+component.
+
+The following calculation summarizes their comovement over the full posterior
+sample.
+
+```{code-cell} ipython3
+core_natural_correlation = np.corrcoef(core_inflation, natural_rate)[0, 1]
+
+core_natural_correlation
+```
+
+The strong positive quarterly correlation between $\bar\pi_t$ and $\bar u_t$
+says that the model-implied long-run inflation and unemployment anchors share
+a broad cycle, not that current inflation and current unemployment must move
+together.
+
+Because $(I-A_t)^{-1}$ amplifies small coefficient changes when the largest root
+is close to one, long-run means are intrinsically more sensitive than
+short-horizon forecasts.
+
+### Inflation persistence
+
+We now ask whether the *systematic* dynamics drifted on top of the moving
+volatilities.
+
+The main summary is inflation persistence, measured by the normalized spectrum
+of inflation at frequency zero.
+
+The spectral density of inflation at date $t$ is
+
+```{math}
+:label: csdv_spectrum
+f_{\pi\pi}(\omega,t)
+=
+\frac{1}{2\pi}
+s_\pi
+(I-A_{t\mid T}e^{-i\omega})^{-1}
+\mathcal R_t
+(I-A_{t\mid T}^\top e^{i\omega})^{-1}
+s_\pi^\top,
+```
+
+where $\mathcal R_t$ embeds $\mathbb{E}(R_t\mid T)$ in the companion system.
+
+Low-frequency power depends on both the autoregressive coefficients and the
+innovation covariance.
+
+The next function evaluates inflation power at any frequency measured in cycles
+per quarter and also returns inflation variance at date $t$.
+
+```{code-cell} ipython3
+def inflation_spectrum(θ, covariance, frequencies):
+ """Compute inflation power and its variance-normalized counterpart."""
+ _, companion = companion_matrix(θ)
+ innovation = np.zeros((6, 6))
+ innovation[:3, :3] = covariance
+ selector = np.zeros(6)
+ selector[2] = 1
+ stationary = linalg.solve_discrete_lyapunov(companion, innovation)
+ variance = float(selector @ stationary @ selector)
+ power = np.empty(len(frequencies))
+ for index, frequency in enumerate(frequencies):
+ phase = np.exp(-2j * np.pi * frequency)
+ transfer = np.linalg.solve(np.eye(6) - companion * phase, np.eye(6))
+ power[index] = np.real(
+ selector @ transfer @ innovation @ transfer.conj().T @ selector
+ ) / (2 * np.pi)
+ return power, power / variance
+
+```
+
+The normalized spectrum divides by inflation variance at date $t$,
+
+```{math}
+:label: csdv_normalized_spectrum
+g_{\pi\pi}(\omega,t)
+=
+\frac{f_{\pi\pi}(\omega,t)}
+{\int_{-\pi}^{\pi}f_{\pi\pi}(\omega,t)d\omega},
+```
+
+so $g_{\pi\pi}(0,t)$ is an autocorrelation-based persistence measure.
+
+The normalization removes a common scale factor from $R_t$, but it can still
+depend on the relative variances and covariances in $R_t$.
+
+The first figure isolates frequency zero as a one-dimensional persistence
+summary.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Normalized zero-frequency spectrum $g_{\pi\pi}(0,t)$
+ name: fig-csdv-inflation-persistence
+---
+zero_frequency = np.array([0.0])
+inflation_persistence = np.array([
+ inflation_spectrum(θ_mean[:, t], r_mean[t], zero_frequency)[1][0]
+ for t in range(len(data['dates']))
+])
+
+fig, ax = plt.subplots()
+ax.plot(
+ data['dates'][annual],
+ inflation_persistence[annual],
+ 'o-',
+ lw=2,
+ markersize=3,
+)
+ax.set_xlabel('year')
+ax.set_ylabel(r'$g_{\pi\pi}(0,t)$')
+plt.show()
+```
+
+```{code-cell} ipython3
+persistence_summary = pd.Series(
+ {
+ '1960--64 mean': inflation_persistence[
+ (years >= 1960) & (years <= 1964)
+ ].mean(),
+ '1970--79 mean': inflation_persistence[
+ (years >= 1970) & (years <= 1979)
+ ].mean(),
+ '1985--2000 mean': inflation_persistence[
+ (years >= 1985) & (years <= 2000)
+ ].mean(),
+ 'peak': inflation_persistence[annual].max(),
+ 'peak year': years[annual][np.argmax(inflation_persistence[annual])],
+ },
+ name='estimate',
+)
+persistence_summary.to_frame().round(3)
+```
+
+Normalized zero-frequency power rises sharply from a low level in the early
+1960s to a high peak around 1980 and then falls back below one for most of the
+remaining sample.
+
+The post-1980 collapse cannot be explained by a proportional rescaling of all
+innovations, although the normalized statistic can still depend on the
+composition of $R_t$.
+
+For comparison, an $AR(1)$ with coefficient $\rho$ has normalized zero-frequency
+power $(1+\rho)/[2\pi(1-\rho)]$.
+
+Values between 2 and 10 correspond to $\rho$ between approximately $0.85$ and
+$0.97$.
+
+The zero-frequency path omits the rest of the frequency distribution.
+
+The following heatmaps show how raw and variance-normalized inflation power move
+over both time and frequency.
+
+```{code-cell} ipython3
+def inflation_spectrum_surface(θ_path, covariance_path, frequencies):
+ """Evaluate the inflation spectrum at each date on a frequency grid."""
+ raw = np.empty((len(frequencies), θ_path.shape[1]))
+ normalized = np.empty_like(raw)
+ for date in range(θ_path.shape[1]):
+ raw[:, date], normalized[:, date] = inflation_spectrum(
+ θ_path[:, date],
+ covariance_path[date],
+ frequencies,
+ )
+ return raw, normalized
+
+
+spectrum_frequencies = np.linspace(0, 0.5, 41)
+raw_spectrum, normalized_spectrum = inflation_spectrum_surface(
+ θ_mean,
+ r_mean,
+ spectrum_frequencies,
+)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Inflation spectra $f_{\pi\pi}(\omega,t)$ and $g_{\pi\pi}(\omega,t)$
+ name: fig-csdv-inflation-spectra
+---
+spectrum_start = data['dates'] >= 1960
+fig, axes = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
+surfaces = (
+ (raw_spectrum, 'raw spectrum', 'log10 power'),
+ (normalized_spectrum, 'normalized spectrum', 'log10 normalized power'),
+)
+for ax, (surface, title, color_label) in zip(axes, surfaces):
+ image = ax.pcolormesh(
+ data['dates'][spectrum_start],
+ spectrum_frequencies,
+ np.log10(surface[:, spectrum_start]),
+ shading='auto',
+ )
+ ax.set_xlabel('year')
+ ax.set_ylabel(f'{title}\ncycles per quarter')
+ fig.colorbar(image, ax=ax, label=color_label)
+plt.tight_layout()
+plt.show()
+```
+
+The raw spectrum is brightest near frequency zero around 1980 because shocks
+and persistence are both elevated, while the normalized spectrum retains a
+broad low-frequency ridge through the 1970s that recedes after 1980.
+
+The ridge that survives normalization shows that the Great Inflation was not
+only a high-volatility episode because inflation shocks were also propagated
+more persistently.
+
+Point estimates do not reveal how strongly the data locate these paths.
+
+We therefore compute the same annual features from every retained draw.
+
+```{code-cell} ipython3
+def posterior_feature_draws(result, indices):
+ """Compute selected local means and persistence for retained draws."""
+ n_draws = result['SD'].shape[2]
+ shape = (len(indices), n_draws)
+ core = np.empty(shape)
+ natural = np.empty(shape)
+ persistence = np.empty(shape)
+ for draw in range(n_draws):
+ core[:, draw], natural[:, draw] = local_means(
+ result['SD'][:, indices, draw]
+ )
+ for row, date in enumerate(indices):
+ covariance = innovation_covariance(
+ result['HD'][date + 1, :, draw],
+ result['CD'][:, draw],
+ )
+ persistence[row, draw] = inflation_spectrum(
+ result['SD'][:, date, draw],
+ covariance,
+ zero_frequency,
+ )[1][0]
+ return {
+ 'core': 100 * core,
+ 'natural': 100 * natural,
+ 'persistence': persistence,
+ }
+
+
+historical_feature_draws = posterior_feature_draws(posterior, annual)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior medians and pointwise 90 percent intervals for $\bar\pi_t$, $\bar u_t$, and $g_{\pi\pi}(0,t)$
+ name: fig-csdv-feature-uncertainty
+---
+fig, axes = plt.subplots(3, 1, figsize=(9, 8), sharex=True)
+feature_specs = (
+ ('core', 'core inflation (%)'),
+ ('natural', 'natural rate (%)'),
+ ('persistence', r'$g_{\pi\pi}(0,t)$'),
+)
+for ax, (key, ylabel) in zip(axes, feature_specs):
+ lower, median, upper = np.quantile(
+ historical_feature_draws[key],
+ (0.05, 0.5, 0.95),
+ axis=1,
+ )
+ line, = ax.plot(data['dates'][annual], median, lw=2)
+ ax.fill_between(
+ data['dates'][annual],
+ lower,
+ upper,
+ color=line.get_color(),
+ alpha=0.2,
+ )
+ ax.set_ylabel(ylabel)
+axes[-1].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+The solid curves are posterior medians, and the shaded regions are pointwise 90
+percent intervals rather than simultaneous bands for entire paths.
+
+The median core-inflation and persistence paths preserve the rise and
+post-1980 fall, but their right-skewed bands widen markedly through the 1970s
+and around 1980.
+
+The natural-rate median moves more smoothly, with broad uncertainty around 1980
+and again at the sample endpoint, so the direction of the historical movement
+is clearer than its exact magnitude.
+
+The following plot compares the timing of core inflation and persistence.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $\bar\pi_t$ and $g_{\pi\pi}(0,t)$
+ name: fig-csdv-core-persistence
+---
+core_persistence_correlation = np.corrcoef(
+ core_inflation, inflation_persistence
+)[0, 1]
+
+fig, ax = plt.subplots()
+ax.plot(data['dates'][annual], 100 * core_inflation[annual], 'o-',
+ lw=2, markersize=3, label='core inflation (%)')
+ax.plot(data['dates'][annual], inflation_persistence[annual], 'x-',
+ lw=2, markersize=4, label='normalized spectrum at zero')
+ax.set_xlabel('year')
+ax.legend()
+plt.show()
+```
+
+Core inflation and normalized persistence rise together through the late 1960s
+and 1970s and collapse almost simultaneously after 1980, although core remains
+positive after persistence falls below one.
+
+Their quarterly correlation of 0.909 summarizes common timing rather than a
+causal relationship because the two measures have different units and are
+nonlinear summaries of the same fitted coefficients.
+
+Within the fitted TVP-VAR, the joint movement supports an important role for
+changing propagation as well as changing shock volatility.
+
+A direct comparison with fixed coefficients requires fitting the restricted
+$Q=0$ model.
+
+The fall in persistence during the Volcker disinflation conflicts with
+escape-route models in which persistence grows along a transition from high to
+low inflation {cite}`Sargent1999,ChoWilliamsSargent2002`.
+
+That tension helped motivate later learning models in which policymakers became
+reluctant to disinflate during the 1970s and then changed course.
+
+### Monetary policy activism
+
+Cogley and Sargent summarize systematic policy with a forward-looking Taylor
+rule,
+
+```{math}
+:label: csdv_policy_rule
+i_t = \beta_0
++ \beta_1 \mathbb{E}_t\bar\pi_{t,t+h_\pi}
++ \beta_2 \mathbb{E}_t\bar u_{t,t+h_u}
++ \beta_3 i_{t-1}
++ \nu_t.
+```
+
+They define the activism coefficient as $\mathcal A_t=\beta_1/(1-\beta_3)$ and
+call policy active when $\mathcal A_t\geq 1$.
+
+At each date, population two-stage least squares projections implied by the
+local VAR produce the policy-rule coefficients.
+
+The benchmark horizons are $h_\pi=4$ quarters and $h_u=2$ quarters, reflecting
+conventional views about monetary-policy lags.
+
+The following function projects the short rate on the model-implied inflation
+and unemployment forecasts using the stationary second moments of each local
+VAR.
+
+```{code-cell} ipython3
+def policy_rule_coefficients(θ, covariance, h_pi=4, h_u=2):
+ """Return the local policy-rule coefficients."""
+ _, companion = companion_matrix(θ)
+ innovation = np.zeros((6, 6))
+ innovation[:3, :3] = covariance
+ stationary_covariance = linalg.solve_discrete_lyapunov(
+ companion, innovation
+ )
+ selectors = np.eye(6)
+ companion_power = np.eye(6)
+ inflation_loading = np.zeros(6)
+ unemployment_loading = np.zeros(6)
+ for horizon in range(1, max(h_pi, h_u) + 1):
+ companion_power = companion_power @ companion
+ if horizon <= h_pi:
+ inflation_loading += selectors[2] @ companion_power
+ if horizon <= h_u:
+ unemployment_loading += selectors[1] @ companion_power
+ inflation_loading /= h_pi
+ unemployment_loading /= h_u
+ loadings = np.vstack(
+ (inflation_loading, unemployment_loading, selectors[0])
+ )
+ regressor_covariance = loadings @ stationary_covariance @ loadings.T
+ cross_covariance = (
+ loadings @ stationary_covariance @ companion.T @ selectors[0]
+ )
+ return np.linalg.solve(regressor_covariance, cross_covariance)
+
+
+policy_coefficients = np.array(
+ [
+ policy_rule_coefficients(θ_mean[:, t], r_mean[t])
+ for t in range(len(data['dates']))
+ ]
+)
+inflation_response = policy_coefficients[:, 0]
+interest_persistence = policy_coefficients[:, 2]
+policy_margin = np.where(
+ np.abs(interest_persistence) < 1,
+ inflation_response + interest_persistence - 1,
+ np.nan,
+)
+```
+
+For $|\beta_3|<1$, the policy margin
+$\mathcal{M}_t=\beta_{1t}+\beta_{3t}-1$ is nonnegative exactly when
+$\mathcal A_t\geq1$ and avoids division by $1-\beta_3$.
+
+Because $\beta_3$ multiplies the lagged interest rate, the long-run response
+sums $\beta_1(1+\beta_3+\beta_3^2+\cdots)$ and exists only when
+$|\beta_3|<1$.
+
+The figures leave other dates blank because $\mathcal A_t$ has no finite
+long-run interpretation there.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Policy margin $\mathcal{M}_t=\beta_{1t}+\beta_{3t}-1$
+ name: fig-csdv-policy-activism
+---
+fig, ax = plt.subplots()
+ax.plot(data['dates'], policy_margin, lw=2)
+ax.axhline(0, color='0.45', lw=1)
+ax.set_xlabel('year')
+ax.set_ylabel(r'policy margin $\mathcal{M}_t$')
+plt.show()
+```
+
+The policy margin falls below zero through much of the 1970s and then moves
+decisively above zero after the early 1980s.
+
+This timing is consistent with a policy-regime contribution to the Great
+Inflation.
+
+Blank intervals mark dates at which the policy margin is omitted by the rule
+above.
+
+The following scatter plots use fourth-quarter observations.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $\mathcal{M}_t$ versus $\bar\pi_t$ and $g_{\pi\pi}(0,t)$
+ name: fig-csdv-activism-correlations
+---
+displayed_annual = annual[np.isfinite(policy_margin[annual])]
+
+fig, axes = plt.subplots(1, 2, figsize=(9, 4))
+pairs = (
+ (100 * core_inflation, 'core inflation (%)'),
+ (inflation_persistence, 'normalized spectrum at zero'),
+)
+for ax, (feature, label) in zip(axes, pairs):
+ ax.scatter(
+ policy_margin[displayed_annual],
+ feature[displayed_annual],
+ s=18,
+ )
+ ax.axvline(0, color='0.45', lw=1)
+ ax.set_xlabel(r'policy margin $\mathcal{M}_t$')
+ ax.set_ylabel(label)
+plt.tight_layout()
+plt.show()
+```
+
+High core-inflation and persistence observations cluster near or below the zero
+margin, whereas positive policy margins cluster at low values of both measures.
+
+The displayed fourth-quarter observations establish a historical association
+rather than a causal policy effect.
+
+The policy-rule coefficients are weakly identified at some dates, so a path
+based on posterior mean inputs understates uncertainty.
+
+We therefore calculate activism from every retained draw in 1975, 1985, and
+1995.
+
+```{code-cell} ipython3
+selected_years = (1975, 1985, 1995)
+selected_dates = [np.flatnonzero(years == year)[-1] for year in selected_years]
+
+activism_draws = {
+ year: np.empty(posterior['SD'].shape[2]) for year in selected_years
+}
+stable_response_by_year = {
+ year: np.empty(posterior['SD'].shape[2], dtype=bool)
+ for year in selected_years
+}
+for draw in range(posterior['SD'].shape[2]):
+ for year, date in zip(selected_years, selected_dates):
+ covariance = innovation_covariance(
+ posterior['HD'][date + 1, :, draw],
+ posterior['CD'][:, draw],
+ )
+ rule = policy_rule_coefficients(
+ posterior['SD'][:, date, draw], covariance
+ )
+ activism_draws[year][draw] = rule[0] / (1 - rule[2])
+ stable_response_by_year[year][draw] = np.abs(rule[2]) < 1
+```
+
+These draws give posterior probabilities of active policy at each date and of a
+rise in activism after 1975, conditional on $|\beta_3|<1$.
+
+```{code-cell} ipython3
+activism_events = (
+ activism_draws[1975] > 1,
+ activism_draws[1985] > 1,
+ activism_draws[1995] > 1,
+ activism_draws[1985] > activism_draws[1975],
+ activism_draws[1995] > activism_draws[1975],
+)
+activism_conditions = (
+ stable_response_by_year[1975],
+ stable_response_by_year[1985],
+ stable_response_by_year[1995],
+ stable_response_by_year[1985] & stable_response_by_year[1975],
+ stable_response_by_year[1995] & stable_response_by_year[1975],
+)
+assert all(condition.any() for condition in activism_conditions)
+activism_probability_values = np.array(
+ [
+ event[condition].mean()
+ for event, condition in zip(activism_events, activism_conditions)
+ ]
+)
+
+activism_probability_index = (
+ 'P(A_1975 > 1)',
+ 'P(A_1985 > 1)',
+ 'P(A_1995 > 1)',
+ 'P(A_1985 > A_1975)',
+ 'P(A_1995 > A_1975)',
+)
+activism_probabilities = pd.DataFrame(
+ {'conditional estimate': activism_probability_values},
+ index=activism_probability_index,
+)
+stable_response_shares = pd.Series(
+ {year: draws.mean() for year, draws in stable_response_by_year.items()},
+ name='share with |beta_3| < 1',
+)
+
+display(activism_probabilities.round(3))
+stable_response_shares.to_frame().round(3)
+```
+
+The first three probabilities condition on $|\beta_3|<1$ at that date, while
+the comparisons require this condition at both dates.
+
+The central draw distributions expose the overlap and skewness behind the
+probability estimates.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Central posterior draws of $\mathcal A_t$ conditional on $|\beta_3|<1$ in 1975, 1985, and 1995
+ name: fig-csdv-activism-distributions
+---
+stable_activism_draws = {
+ year: activism_draws[year][stable_response_by_year[year]]
+ for year in selected_years
+}
+pooled_activism = np.concatenate(tuple(stable_activism_draws.values()))
+activism_limits = np.quantile(pooled_activism, (0.05, 0.95))
+activism_bins = np.linspace(*activism_limits, 31)
+
+fig, ax = plt.subplots()
+for year in selected_years:
+ central_draws = stable_activism_draws[year]
+ central_draws = central_draws[
+ (central_draws >= activism_limits[0])
+ & (central_draws <= activism_limits[1])
+ ]
+ ax.hist(
+ central_draws,
+ bins=activism_bins,
+ histtype='step',
+ lw=2,
+ label=str(year),
+ )
+ax.axvline(1, color='0.45', lw=1)
+ax.set_xlabel('activism coefficient')
+ax.set_ylabel('retained draws')
+ax.legend()
+plt.show()
+```
+
+The 1975 distribution is concentrated around or below one, while the 1985 and
+1995 distributions shift substantially to the right but remain broad, skewed,
+and overlapping.
+
+The posterior therefore favors a passive-to-active shift after 1975 but does
+not sharply distinguish 1985 from 1995.
+
+The figure plots draws with $|\beta_3|<1$ inside the pooled 5th and 95th
+percentiles, while the probability calculations use every such draw.
+
+(csdv-updated-evidence)=
+## Another quarter-century of evidence
+
+The sample ends in 2000Q4, so it misses the financial crisis, the
+zero-interest-rate period, the pandemic, and the 2021--2022 inflation surge.
+
+The question is whether these episodes alter the earlier evidence about drift,
+volatility, inflation persistence, and systematic policy.
+
+### The new observations
+
+To examine those observations, we append current data after 2000Q4 while leaving pre-2000Q4 sample unchanged.
+
+This splice prevents revisions to pre-2001 CPI and unemployment data from being
+mistaken for information in the additional quarter-century.
+
+We download seasonally adjusted [CPI][fred-cpi], seasonally adjusted
+[unemployment][fred-unemployment], and the [three-month Treasury-bill
+rate][fred-interest] from FRED.
+
+[fred-cpi]: https://fred.stlouisfed.org/series/CPIAUCSL
+[fred-unemployment]: https://fred.stlouisfed.org/series/UNRATE
+[fred-interest]: https://fred.stlouisfed.org/series/TB3MS
+
+The transformations and within-quarter timing remain unchanged: CPI comes from
+the third month, unemployment is a three-month average, and the interest rate
+comes from the first month.
+
+```{code-cell} ipython3
+fred_url = (
+ 'https://fred.stlouisfed.org/graph/fredgraph.csv?'
+ 'id=CPIAUCSL%2CUNRATE%2CTB3MS'
+)
+fred_monthly = pd.read_csv(
+ fred_url,
+ parse_dates=['observation_date'],
+).set_index('observation_date')
+
+```
+
+The [BLS notes](https://www.bls.gov/web/empsit/cpsee_e12.pdf) that the October
+2025 unemployment observation was not collected during the federal government
+shutdown, leaving 2025Q4 without a complete three-month average.
+
+The code below therefore ends the updated sample at the last quarter whose three
+monthly unemployment readings are all present, detected automatically rather
+than hard-coded; at the time of writing that quarter is 2025Q3.
+
+```{code-cell} ipython3
+def fred_quarterly_table(unemployment_monthly):
+ """Construct transformed quarterly observations from current FRED data."""
+ interest = fred_monthly.loc[
+ fred_monthly.index.month.isin((1, 4, 7, 10)), 'TB3MS'
+ ].copy()
+ interest.index = interest.index.to_period('Q').start_time
+
+ cpi = fred_monthly.loc[
+ fred_monthly.index.month.isin((3, 6, 9, 12)), 'CPIAUCSL'
+ ].copy()
+ cpi.index = cpi.index.to_period('Q').start_time
+
+ unemployment = unemployment_monthly.resample('QS').mean()
+ quarterly = pd.concat(
+ {
+ 'interest': interest,
+ 'unemployment': unemployment,
+ 'cpi': cpi,
+ },
+ axis=1,
+ )
+ quarterly['y3'] = np.log1p(quarterly['interest'] / 400)
+ quarterly['ur'] = quarterly['unemployment'] / 100
+ quarterly['dp'] = np.log(quarterly['cpi']).diff()
+ quarterly['date'] = (
+ quarterly.index.year + (quarterly.index.quarter - 1) / 4
+ )
+ columns = ['date', 'y3', 'ur', 'dp']
+ return quarterly.loc[:, columns].dropna()
+
+
+unemployment_monthly = fred_monthly['UNRATE']
+unemployment_counts = unemployment_monthly.resample('QS').count()
+latest_quarterly = fred_quarterly_table(unemployment_monthly)
+
+# Extend the archived sample dynamically: append post-2000Q4 quarters only up
+# to the first one missing any of its three monthly unemployment readings, so
+# no endpoint is hard-coded. An isolated missing month -- for example October
+# 2025, which the federal shutdown left uncollected -- caps the sample at the
+# preceding complete quarter, and a normally incomplete current quarter caps it
+# at the last finished one.
+quarterly_unfilled = latest_quarterly.loc['2001-01-01':]
+month_counts = unemployment_counts.reindex(quarterly_unfilled.index)
+incomplete_quarters = month_counts.index[month_counts < 3]
+if len(incomplete_quarters):
+ first_incomplete_quarter = incomplete_quarters[0]
+ complete_extension = quarterly_unfilled.loc[
+ quarterly_unfilled.index < first_incomplete_quarter
+ ]
+else:
+ complete_extension = quarterly_unfilled
+
+cs_sample = pd.read_csv(data_url)
+overlap_date = pd.Timestamp('2000-10-01')
+cs_sample_overlap = cs_sample.iloc[-1]
+latest_overlap = latest_quarterly.loc[overlap_date]
+
+
+def scaled_observation(row):
+ """Return observable units for a transformed quarterly row."""
+ return pd.Series(
+ {
+ 'interest rate (annual %)': 400 * np.expm1(row['y3']),
+ 'unemployment (%)': 100 * row['ur'],
+ 'inflation (annual %)': 400 * row['dp'],
+ }
+ )
+
+
+cs_sample_scaled = scaled_observation(cs_sample_overlap)
+latest_scaled = scaled_observation(latest_overlap)
+splice_audit = pd.DataFrame(
+ {
+ 'Cogley-Sargent (2005) 2000Q4': cs_sample_scaled,
+ 'latest revised 2000Q4 data': latest_scaled,
+ 'current minus Cogley-Sargent (2005)': (
+ latest_scaled - cs_sample_scaled
+ ),
+ 'first appended 2001Q1': scaled_observation(
+ complete_extension.iloc[0]
+ ),
+ }
+)
+extended_observations = pd.concat(
+ (cs_sample, complete_extension.reset_index(drop=True)),
+ ignore_index=True,
+)
+
+
+def quarter_label(timestamp):
+ """Format a timestamp as year and quarter."""
+ return str(timestamp.to_period('Q'))
+
+
+display(splice_audit.round(3))
+```
+
+No missing value is filled or otherwise treated as observed.
+
+The overlap table shows any break created by joining the archived data to the
+newly downloaded data.
+
+The appended 2001Q1 inflation rate uses the latest revised CPI level for both
+2000Q4 and 2001Q1, so one log difference never combines observations from two
+data releases.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Observed $\pi_t$, $u_t$, and $i_t$, 1959Q1--2025Q3
+ name: fig-csdv-updated-data
+---
+fig, axes = plt.subplots(3, 1, figsize=(9, 7), sharex=True)
+extended_plot = extended_observations[
+ extended_observations['date'] >= 1959
+]
+updated_sample_end = complete_extension['date'].iloc[-1]
+axes[0].plot(
+ extended_plot['date'],
+ 400 * extended_plot['dp'],
+ lw=2,
+)
+axes[0].set_ylabel('inflation (annual %)')
+axes[1].plot(
+ extended_plot['date'],
+ 100 * extended_plot['ur'],
+ lw=2,
+)
+axes[1].set_ylabel('unemployment (%)')
+axes[2].plot(
+ extended_plot['date'],
+ 400 * np.expm1(extended_plot['y3']),
+ lw=2,
+)
+axes[2].set_ylabel('interest (annual %)')
+axes[2].set_xlabel('year')
+for index, ax in enumerate(axes):
+ cs_sample_label = 'Cogley-Sargent (2005)' if index == 0 else None
+ sample_end_label = 'updated sample end' if index == 0 else None
+ ax.axvline(
+ 2000.75,
+ color='0.65',
+ ls='--',
+ lw=1,
+ label=cs_sample_label,
+ )
+ ax.axvline(
+ updated_sample_end,
+ color='0.45',
+ ls=':',
+ lw=1,
+ label=sample_end_label,
+ )
+axes[0].legend()
+plt.tight_layout()
+plt.show()
+```
+
+```{code-cell} ipython3
+endpoint_observation = complete_extension.iloc[-1]
+pd.Series(
+ {
+ 'annualized quarterly inflation (%)': (
+ 400 * endpoint_observation['dp']
+ ),
+ 'unemployment (%)': 100 * endpoint_observation['ur'],
+ 'three-month interest rate (%)': (
+ 400 * np.expm1(endpoint_observation['y3'])
+ ),
+ },
+ name=quarter_label(complete_extension.index[-1]),
+).to_frame().round(3)
+```
+
+The extension adds the financial-crisis contraction, a pandemic quarterly
+unemployment spike to 13.0 percent, and a short 2021--2022 inflation surge
+alongside two long stretches of near-zero interest rates.
+
+Because inflation is annualized from quarterly changes, isolated movements look
+especially large in this panel, but the recent surge is still visibly much
+shorter than the sustained 1970s rise.
+
+These observations provide a demanding test of whether the model assigns recent
+extremes to shock volatility or persistent dynamics, while the near-zero rate
+also weakens short-rate measures of policy after 2008.
+
+We fit the stable TVP-VAR through 2025Q3, the last complete quarter, using the
+Treasury-bill measure above.
+
+```{code-cell} ipython3
+def append_extension(extension):
+ """Append a transformed FRED extension to the Cogley-Sargent sample."""
+ extension = extension.reset_index(drop=True)
+ table = pd.concat((cs_sample, extension), ignore_index=True)
+ assert table['date'].is_unique
+ assert np.allclose(np.diff(table['date']), 0.25)
+ return table
+
+
+def fit_updated_model(table):
+ """Calibrate and fit the stable drifting VAR to one table."""
+ model_data = prepare_data(table)
+ model_prior = calibrate_prior(model_data)
+ result = run_sampler(
+ model_data['y'],
+ model_data['x'],
+ model_prior,
+ n_sweeps=5_000,
+ burn=2_500,
+ thin=5,
+ seed=42,
+ warmup=500,
+ stable=True,
+ progress_every=500,
+ )
+ validate_posterior_arrays(result, len(model_data['dates']))
+ trace = np.trace(result['QD'], axis1=0, axis2=1)
+ return {
+ 'data': model_data,
+ 'prior': model_prior,
+ 'posterior': result,
+ 'trace': trace,
+ }
+
+
+updated_fit = fit_updated_model(append_extension(complete_extension))
+
+latest_data_table = latest_quarterly.loc[
+ (latest_quarterly.index >= pd.Timestamp('1948-04-01'))
+ & (latest_quarterly.index <= complete_extension.index[-1])
+].reset_index(drop=True)
+assert np.allclose(np.diff(latest_data_table['date']), 0.25)
+latest_data_fit = fit_updated_model(latest_data_table)
+```
+
+### Did coefficient drift continue?
+
+```{code-cell} ipython3
+def updated_drift_summary(fit):
+ """Summarize the updated coefficient-drift distribution."""
+ trace = fit['trace']
+ q_mean = fit['posterior']['QD'].mean(axis=2)
+ eigenvalues = np.linalg.eigvalsh(q_mean)[::-1]
+ return {
+ 'posterior mean tr(Q)': trace.mean(),
+ 'share in first three eigen-directions': (
+ eigenvalues[:3].sum() / eigenvalues.sum()
+ ),
+ }
+
+
+updated_label = quarter_label(complete_extension.index[-1])
+updated_summary = pd.DataFrame(
+ {
+ 'Cogley-Sargent (2005) + extension': (
+ updated_drift_summary(updated_fit)
+ ),
+ 'latest revised data for full sample': updated_drift_summary(
+ latest_data_fit
+ ),
+ }
+)
+updated_summary.round(3)
+```
+
+The drift-rate distributions and coefficient paths provide the same views used
+for the {cite:t}`CogleySargent2005` sample.
+
+The dashed vertical line in updated time-series figures marks the
+{cite:t}`CogleySargent2005` sample's 2000Q4 endpoint, but each updated path
+comes from a full re-estimation rather than from attaching new points to an
+unchanged historical estimate.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior $\operatorname{tr}(Q)$ and prior $\operatorname{tr}(\bar Q)$ through 2025Q3
+ name: fig-csdv-updated-drift-rate
+---
+fig, ax = plt.subplots()
+ax.hist(updated_fit['trace'], bins=30, histtype='step', lw=2)
+ax.axvline(
+ np.trace(updated_fit['prior']['q_center']),
+ color='C1',
+ lw=2,
+ label=r'prior $\mathrm{tr}(\bar Q)$',
+)
+ax.set_xlabel(r'$\mathrm{tr}(Q)$')
+ax.set_ylabel('frequency')
+ax.legend()
+plt.show()
+```
+
+In both data constructions, the full-sample posterior drift rate remains above
+its conservative prior calibration.
+
+Within the TVP-VAR this indicates non-negligible full-sample drift, but it is
+not a formal comparison of $Q=0$ with $Q>0$.
+
+Its magnitude depends on whether the historical observations come from the
+archived dataset or the latest revisions.
+
+Because $Q$ is a single variance parameter for the full 1959--2025 path, this
+histogram does not by itself show that drift accelerated after 2000.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior mean VAR coefficients $E(\theta_t\mid T)$ through 2025Q3
+ name: fig-csdv-updated-coefficient-paths
+---
+fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharex=True)
+updated_dates = updated_fit['data']['dates']
+updated_coefficients = updated_fit['posterior']['SD'].mean(axis=2)
+coefficient_lines = plot_equation_coefficients(
+ axes,
+ updated_dates,
+ updated_coefficients,
+)
+for ax in axes:
+ ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+fig.legend(
+ coefficient_lines,
+ coefficient_labels,
+ loc='lower center',
+ ncol=4,
+)
+plt.tight_layout(rect=(0, 0.17, 1, 1))
+plt.show()
+```
+
+Several coefficient paths continue moving gradually after 2000, with the
+largest changes again concentrated in the inflation equation rather than
+appearing as abrupt financial-crisis or pandemic breaks.
+
+The opposing movements of some paired lag coefficients also show why their
+combined dynamic implications are more informative than any single line.
+
+The first three eigen-directions account for the large majority of the
+posterior mean drift variance, so the estimated movement remains low
+dimensional.
+
+The stability restriction now applies to a longer path, which makes the
+posterior drift rate a property of the full 1959--2025 sample.
+
+### Volatility after the Great Moderation
+
+We next separate changes in the sizes and correlations of innovations from
+changes in the VAR dynamics.
+
+```{code-cell} ipython3
+def model_features(fit):
+ """Compute features at posterior mean parameters."""
+ result = fit['posterior']
+ θ = result['SD'].mean(axis=2)
+ mean_root_modulus = np.max(
+ np.abs(companion_roots(θ)), axis=1
+ )
+ if not np.all(mean_root_modulus < 1):
+ raise ValueError('mean coefficient path is unstable')
+ covariance = mean_innovation_covariance(result['HD'], result['CD'])
+ core, natural = local_means(θ)
+ persistence = np.array([
+ inflation_spectrum(
+ θ[:, t], covariance[t], zero_frequency
+ )[1][0]
+ for t in range(len(fit['data']['dates']))
+ ])
+ sign, logdet = np.linalg.slogdet(covariance)
+ assert np.all(sign > 0)
+ policy_coefficients = np.array(
+ [
+ policy_rule_coefficients(θ[:, t], covariance[t])
+ for t in range(len(fit['data']['dates']))
+ ]
+ )
+ inflation_response = policy_coefficients[:, 0]
+ interest_persistence = policy_coefficients[:, 2]
+ policy_margin = np.where(
+ np.abs(interest_persistence) < 1,
+ inflation_response + interest_persistence - 1,
+ np.nan,
+ )
+ return {
+ 'θ': θ,
+ 'maximum_companion_root': mean_root_modulus.max(),
+ 'covariance': covariance,
+ 'core': core,
+ 'natural': natural,
+ 'persistence': persistence,
+ 'logdet': logdet,
+ 'policy_margin': policy_margin,
+ }
+
+
+updated_features = model_features(updated_fit)
+latest_data_features = model_features(latest_data_fit)
+
+pd.Series(
+ {
+ 'Cogley-Sargent (2005) + extension': (
+ updated_features['maximum_companion_root']
+ ),
+ 'latest revised data for full sample': (
+ latest_data_features['maximum_companion_root']
+ ),
+ },
+ name='maximum companion-root modulus of mean path',
+).to_frame().round(4)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Standard deviations and correlations implied by $E(R_t\mid T)$ through 2025Q3
+ name: fig-csdv-updated-volatility-correlation
+---
+fig, axes = plt.subplots(3, 2, figsize=(9, 8), sharex=True)
+dates = updated_fit['data']['dates']
+covariance = updated_features['covariance']
+for row, (index, _) in enumerate(variances):
+ axes[row, 0].plot(
+ dates,
+ 10000 * np.sqrt(covariance[:, index, index]),
+ lw=2,
+ )
+for row, (left, right, _) in enumerate(correlations):
+ scale = np.sqrt(
+ covariance[:, left, left] * covariance[:, right, right]
+ )
+ axes[row, 1].plot(
+ dates,
+ covariance[:, left, right] / scale,
+ lw=2,
+ )
+for row, (_, label) in enumerate(variances):
+ axes[row, 0].set_title(label)
+for row, (_, _, label) in enumerate(correlations):
+ axes[row, 1].set_title(label)
+for ax in axes.flat:
+ ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+axes[1, 0].set_ylabel(
+ r'innovation standard deviation $\times 10^4$'
+)
+axes[1, 1].set_ylabel('correlation')
+axes[-1, 0].set_xlabel('year')
+axes[-1, 1].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+The Volcker transition still dominates interest-rate innovation volatility, the
+financial crisis produces the largest inflation-volatility spike, and the
+pandemic uniquely dominates unemployment volatility.
+
+The pandemic also drives a sharp fall in the inflation--unemployment
+correlation, so recent episodes changed the mix of reduced-form shocks as well
+as their size.
+
+This is direct evidence for a bad-luck component, although reduced-form
+innovations do not establish that the underlying disturbances were structurally
+exogenous.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $\log |E(R_t\mid T)|$ through 2025Q3
+ name: fig-csdv-updated-total-variance
+---
+fig, ax = plt.subplots()
+ax.plot(
+ updated_fit['data']['dates'],
+ updated_features['logdet'],
+ lw=2,
+)
+ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+ax.set_xlabel('year')
+ax.set_ylabel(r'$\log |E(R_t\mid T)|$')
+plt.show()
+```
+
+Joint innovation variance rises during the financial crisis, reaches a deep
+Great Moderation trough in the 2010s, and then jumps in 2020 above even its 1981
+peak before falling rapidly.
+
+```{code-cell} ipython3
+def decimal_quarter_label(value):
+ """Format a decimal quarterly date."""
+ year = int(np.floor(value + 1e-8))
+ quarter = int(round(4 * (value - year))) + 1
+ return f'{year}Q{quarter}'
+
+
+def volatility_summary(fit, features):
+ """Summarize innovation-volatility peaks and endpoints."""
+ dates = fit['data']['dates']
+ standard_deviation = 10000 * np.sqrt(
+ np.diagonal(features['covariance'], axis1=1, axis2=2)
+ )
+ names = ('interest', 'unemployment', 'inflation')
+ return pd.DataFrame(
+ {
+ 'peak quarter': [
+ decimal_quarter_label(
+ dates[np.argmax(standard_deviation[:, index])]
+ )
+ for index in range(n_variables)
+ ],
+ },
+ index=names,
+ )
+
+
+updated_volatility = volatility_summary(updated_fit, updated_features)
+updated_volatility
+```
+
+All three innovation standard deviations are well below their peaks at the
+2025Q3 endpoint, which supports a large but transient pandemic-shock
+interpretation inside this model.
+
+### Core inflation and the natural rate after 2000
+
+The same local-mean calculation distinguishes temporary inflation from a shift
+in the model's long-horizon forecast.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Local means $\bar\pi_t$ and $\bar u_t$ through 2025Q3
+ name: fig-csdv-updated-local-means
+---
+updated_dates = updated_fit['data']['dates']
+updated_annual = annual_indices(updated_dates)
+fig, ax = plt.subplots()
+ax.plot(
+ updated_dates[updated_annual],
+ 100 * updated_features['core'][updated_annual],
+ 'o-',
+ lw=2,
+ markersize=3,
+ label='core inflation',
+)
+ax.plot(
+ updated_dates[updated_annual],
+ 100 * updated_features['natural'][updated_annual],
+ '+-',
+ lw=2,
+ markersize=5,
+ label='natural rate',
+)
+ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+ax.set_xlabel('year')
+ax.set_ylabel('percent')
+ax.legend()
+plt.show()
+```
+
+The 1970s core-inflation peak is lower here than in the
+{cite:t}`CogleySargent2005` figure because later observations revise the
+smoothed history, so values on both sides of the dashed line belong to one
+updated fit.
+
+After 2000 the core rate stays mostly between about 2 and 3 percent and
+rises only modestly after 2020, even though observed inflation moves much more
+sharply.
+
+The natural-rate line is the model's locally implied long-run unemployment
+anchor, which is why quarterly unemployment can jump to 13.0 percent in 2020
+while this line remains near 5 percent.
+
+The monthly peak was 14.8 percent.
+
+The final annual point represents 2025Q3 because the sample ends before the
+fourth quarter.
+
+```{code-cell} ipython3
+def endpoint_features(features):
+ """Return economically scaled endpoint features."""
+ return pd.Series(
+ {
+ 'core inflation (%)': 100 * features['core'][-1],
+ 'natural rate (%)': 100 * features['natural'][-1],
+ 'normalized persistence': features['persistence'][-1],
+ 'log generalized innovation variance': features['logdet'][-1],
+ }
+ )
+
+
+updated_endpoints = endpoint_features(updated_features).to_frame(
+ name=updated_label
+).T
+latest_data_endpoints = endpoint_features(
+ latest_data_features
+).to_frame(name=updated_label).T
+data_revision_sensitivity = pd.concat(
+ {
+ 'Cogley-Sargent (2005) + extension': updated_endpoints,
+ 'latest revised data for full sample': latest_data_endpoints,
+ }
+)
+data_revision_sensitivity.round(3)
+```
+
+The endpoint core-inflation, natural-rate, and persistence summaries are similar
+across the two data constructions.
+
+The draw-wise annual paths show how uncertainty evolves inside the updated fit.
+
+```{code-cell} ipython3
+updated_feature_draws = posterior_feature_draws(
+ updated_fit['posterior'],
+ updated_annual,
+)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Posterior medians and pointwise 90 percent intervals for $\bar\pi_t$, $\bar u_t$, and $g_{\pi\pi}(0,t)$ through 2025Q3
+ name: fig-csdv-updated-feature-uncertainty
+---
+fig, axes = plt.subplots(3, 1, figsize=(9, 8), sharex=True)
+for ax, (key, ylabel) in zip(axes, feature_specs):
+ lower, median, upper = np.quantile(
+ updated_feature_draws[key],
+ (0.05, 0.5, 0.95),
+ axis=1,
+ )
+ line, = ax.plot(updated_dates[updated_annual], median, lw=2)
+ ax.fill_between(
+ updated_dates[updated_annual],
+ lower,
+ upper,
+ color=line.get_color(),
+ alpha=0.2,
+ )
+ ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+ ax.set_ylabel(ylabel)
+axes[-1].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+The solid lines are posterior medians, and the shaded regions are pointwise 90
+percent intervals rather than simultaneous whole-path bands.
+
+After 2000 the core-inflation median is comparatively flat and the persistence
+median remains low, whereas the natural-rate band is broad and widens again at
+the endpoint.
+
+The endpoint rows summarize uncertainty in the three nonlinear features.
+
+```{code-cell} ipython3
+def endpoint_feature_intervals(draws):
+ """Summarize draw-wise uncertainty in endpoint model features."""
+ labels = {
+ 'core': 'core inflation (%)',
+ 'natural': 'natural rate (%)',
+ 'persistence': 'normalized persistence',
+ }
+ rows = {}
+ for key, label in labels.items():
+ values = draws[key][-1]
+ rows[label] = {
+ 'median': np.median(values),
+ '5th percentile': np.quantile(values, 0.05),
+ '95th percentile': np.quantile(values, 0.95),
+ }
+ return pd.DataFrame.from_dict(rows, orient='index')
+
+
+updated_endpoint_intervals = endpoint_feature_intervals(
+ updated_feature_draws
+)
+updated_endpoint_intervals.round(3)
+```
+
+At 2025Q3 both the core-inflation and natural-rate medians remain
+historically moderate, as the table above shows.
+
+Their intervals remain broad, and normalized persistence retains a substantial
+upper tail, so the classification of recent inflation as temporary is not
+certain.
+
+### Did inflation become persistent again?
+
+Here persistence means propagation of an inflation innovation into future
+inflation, rather than the number of quarters in which observed inflation
+remains high.
+
+The normalized spectrum reduces sensitivity to a common change in shock scale,
+although it still depends on the relative variances and covariances in $R_t$.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $g_{\pi\pi}(0,t)$ through 2025Q3
+ name: fig-csdv-updated-persistence
+---
+fig, ax = plt.subplots()
+ax.plot(
+ updated_dates[updated_annual],
+ updated_features['persistence'][updated_annual],
+ 'o-',
+ lw=2,
+ markersize=3,
+)
+ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+ax.set_xlabel('year')
+ax.set_ylabel(r'$g_{\pi\pi}(0,t)$')
+plt.show()
+```
+
+The estimated $g_{\pi\pi}(0,t)$ path recreates the rise to a 1980 peak and the
+subsequent collapse, but it stays near 0.2--0.5 after 2000 and rises only
+slightly after 2020.
+
+A sequence of large reduced-form innovations can keep observed inflation high
+for several quarters without generating the strong propagation estimated for
+the 1970s.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $\bar\pi_t$ and $g_{\pi\pi}(0,t)$ through 2025Q3
+ name: fig-csdv-updated-core-persistence
+---
+fig, axes = plt.subplots(
+ 2,
+ 1,
+ figsize=(9, 6),
+ sharex=True,
+)
+axes[0].plot(
+ updated_dates[updated_annual],
+ 100 * updated_features['core'][updated_annual],
+ 'o-',
+ lw=2,
+ markersize=3,
+)
+axes[1].plot(
+ updated_dates[updated_annual],
+ updated_features['persistence'][updated_annual],
+ 'x-',
+ lw=2,
+ markersize=4,
+)
+for ax in axes:
+ ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+axes[0].set_ylabel('core inflation (%)')
+axes[1].set_ylabel('normalized spectrum at zero')
+axes[1].set_xlabel('year')
+plt.tight_layout()
+plt.show()
+```
+
+Core inflation recovers from its mid-2010s low toward its earlier level, while
+persistence remains in its low post-1980 range instead of rising with it.
+
+This divergence separates a modest rise in the model's long-run inflation rate
+from a return to 1970s-style propagation.
+
+The full spectrum shows where the difference comes from.
+
+```{code-cell} ipython3
+updated_raw_spectrum, updated_normalized_spectrum = (
+ inflation_spectrum_surface(
+ updated_features['θ'],
+ updated_features['covariance'],
+ spectrum_frequencies,
+ )
+)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: $f_{\pi\pi}(\omega,t)$ and $g_{\pi\pi}(\omega,t)$ through 2025Q3
+ name: fig-csdv-updated-spectra
+---
+fig, axes = plt.subplots(
+ 1,
+ 2,
+ figsize=(9, 4),
+ sharey=True,
+ constrained_layout=True,
+)
+updated_surfaces = (
+ (updated_raw_spectrum, 'raw spectrum', 'log10 power'),
+ (
+ updated_normalized_spectrum,
+ 'normalized spectrum',
+ 'log10 normalized power',
+ ),
+)
+for ax, (surface, title, color_label) in zip(axes, updated_surfaces):
+ image = ax.pcolormesh(
+ updated_dates,
+ spectrum_frequencies,
+ np.log10(surface),
+ shading='auto',
+ )
+ ax.axvline(2000.75, color='0.85', ls='--', lw=1)
+ ax.set_xlabel('year')
+ ax.set_ylabel(f'{title}\ncycles per quarter')
+ fig.colorbar(image, ax=ax, label=color_label)
+plt.show()
+```
+
+The financial crisis and pandemic appear as bright, broad bands in
+$f_{\pi\pi}(\omega,t)$ because reduced-form innovation variance increased.
+
+The normalized spectrum $g_{\pi\pi}(\omega,t)$ lacks a post-2000 low-frequency
+ridge comparable to the 1970s.
+
+Together, the panels weigh against a return to 1970s-style persistence without
+making the normalized statistic independent of $R_t$.
+
+```{code-cell} ipython3
+def episode_summary(fit, features):
+ """Average selected features over economically distinct episodes."""
+ years = np.floor(fit['data']['dates'] + 1e-8).astype(int)
+ periods = {
+ '1970--1979': (years >= 1970) & (years <= 1979),
+ '1985--2000': (years >= 1985) & (years <= 2000),
+ '2001--2019': (years >= 2001) & (years <= 2019),
+ '2020--2022': (years >= 2020) & (years <= 2022),
+ '2023--2025Q3': years >= 2023,
+ }
+ rows = {}
+ for label, mask in periods.items():
+ rows[label] = {
+ 'core inflation (%)': 100 * features['core'][mask].mean(),
+ 'normalized persistence': features['persistence'][mask].mean(),
+ 'log generalized innovation variance': (
+ features['logdet'][mask].mean()
+ ),
+ }
+ return pd.DataFrame.from_dict(rows, orient='index')
+
+
+updated_episodes = episode_summary(updated_fit, updated_features)
+updated_episodes.round(3)
+```
+
+The episode averages contrast the high volatility of 2020--2022 with the lower
+normalized persistence of the post-2000 decades.
+
+### Can recent policy activism be measured?
+
+When $|\beta_3|<1$, the policy margin gives the same active-policy
+classification after 2000 without division by $1-\beta_3$.
+
+In the following figures, the gray line marks the active-policy threshold
+$\mathcal{M}_t=0$.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Policy margin $\mathcal{M}_t$ through 2025Q3
+ name: fig-csdv-updated-activism
+---
+fig, ax = plt.subplots()
+updated_policy_margin = updated_features['policy_margin']
+ax.plot(updated_dates, updated_policy_margin, lw=2)
+ax.axhline(0, color='0.45', lw=1)
+ax.axvline(2000.75, color='0.65', ls='--', lw=1)
+ax.set_xlabel('year')
+ax.set_ylabel(r'policy margin $\mathcal{M}_t$')
+plt.show()
+```
+
+Among the displayed post-2000 dates, the margin is mostly positive and moves
+toward zero in the mid-2010s and around 2020.
+
+Blank intervals omit dates for which the fitted interest-rate response does not
+settle down.
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Post-2000 $\mathcal{M}_t$ versus $\bar\pi_t$ and $g_{\pi\pi}(0,t)$
+ name: fig-csdv-updated-activism-correlations
+---
+fig, axes = plt.subplots(1, 2, figsize=(9, 4))
+recent_annual = annual_indices(updated_dates, start=2001)
+displayed_recent = recent_annual[
+ np.isfinite(updated_policy_margin[recent_annual])
+]
+axes[0].scatter(
+ updated_policy_margin[displayed_recent],
+ 100 * updated_features['core'][displayed_recent],
+ s=18,
+)
+axes[1].scatter(
+ updated_policy_margin[displayed_recent],
+ updated_features['persistence'][displayed_recent],
+ s=18,
+)
+for ax in axes:
+ ax.axvline(0, color='0.45', lw=1)
+ ax.set_xlabel(r'policy margin $\mathcal{M}_t$')
+axes[0].set_ylabel('core inflation (%)')
+axes[1].set_ylabel('normalized spectrum at zero')
+plt.tight_layout()
+plt.show()
+```
+
+The remaining post-2000 observations do not establish a stable causal relation
+between the policy margin, core inflation, and persistence.
+
+We also examine the central draw distribution at the 2025Q3 endpoint.
+
+```{code-cell} ipython3
+def endpoint_policy_margin_draws(fit):
+ """Return endpoint margins and stability indicators."""
+ result = fit['posterior']
+ margins = np.empty(result['SD'].shape[2])
+ stable_response = np.empty(result['SD'].shape[2], dtype=bool)
+ for draw in range(len(margins)):
+ covariance = innovation_covariance(
+ result['HD'][-1, :, draw],
+ result['CD'][:, draw],
+ )
+ rule = policy_rule_coefficients(
+ result['SD'][:, -1, draw], covariance
+ )
+ margins[draw] = rule[0] + rule[2] - 1
+ stable_response[draw] = np.abs(rule[2]) < 1
+ return margins, stable_response
+
+
+updated_margin_draws, stable_response_draws = (
+ endpoint_policy_margin_draws(updated_fit)
+)
+assert np.any(stable_response_draws)
+stable_margin_draws = updated_margin_draws[stable_response_draws]
+updated_margin_summary = pd.Series(
+ {
+ 'median M': np.median(stable_margin_draws),
+ '5th percentile of M': np.quantile(stable_margin_draws, 0.05),
+ '95th percentile of M': np.quantile(stable_margin_draws, 0.95),
+ 'P(M >= 0 | |beta_3| < 1)': np.mean(stable_margin_draws >= 0),
+ 'share with |beta_3| < 1': stable_response_draws.mean(),
+ },
+ name=updated_label,
+).to_frame().T
+updated_margin_summary.round(3)
+```
+
+```{code-cell} ipython3
+---
+mystnb:
+ figure:
+ caption: Central 90 percent of posterior draws for $\mathcal{M}_t$ conditional on $|\beta_3|<1$ at 2025Q3
+ name: fig-csdv-updated-activism-distributions
+---
+updated_margin_limits = np.quantile(
+ stable_margin_draws,
+ (0.05, 0.95),
+)
+updated_margin_bins = np.linspace(*updated_margin_limits, 31)
+
+fig, ax = plt.subplots()
+central_draws = stable_margin_draws[
+ (stable_margin_draws >= updated_margin_limits[0])
+ & (stable_margin_draws <= updated_margin_limits[1])
+]
+ax.hist(
+ central_draws,
+ bins=updated_margin_bins,
+ histtype='step',
+ lw=2,
+)
+ax.axvline(0, color='0.45', lw=1)
+ax.set_xlabel(r'policy margin $\mathcal{M}_t$')
+ax.set_ylabel('draw count')
+plt.show()
+```
+
+The table reports the share of draws with a stable interest-rate response and
+the active-policy probability among those draws.
+
+An interval spanning zero indicates that the endpoint classification remains
+uncertain.
+
+The plot displays only draws between the 5th and 95th percentiles, and the zero
+lower bound and unconventional policy further weaken the interpretation of this
+short-rate projection after 2008.
+
+### What the additional observations change
+
+The extra quarter-century adds a dramatic volatility episode without a
+post-2000 low-frequency ridge comparable to the 1970s.
+
+The pandemic is the dominant aggregate uncertainty episode, but the recent
+inflation surge does not reproduce the 1970s low-frequency persistence ridge.
+
+Within the updated TVP-VAR, the full-sample drift-rate posterior remains above
+its conservative prior calibration; a fixed-coefficient comparison would
+require a separate $Q=0$ model.
+
+The 2025Q3 natural-rate and policy-margin estimates remain imprecise, especially
+because not every posterior draw satisfies $|\beta_3|<1$.
+
+(csdv-verdict)=
+## Bad policy or bad luck? A verdict
+
+The Bayesian VAR delivers a nuanced answer to the question that opened this
+lecture.
+
+- *Volatilities drifted:* the size of the shocks changed enormously, with a
+ Volcker-era spike and a subsequent Great Moderation, so the bad-luck story
+ captures something real.
+
+- *The fitted TVP-VAR attributes variation to coefficients too:* inflation
+ persistence and core inflation rose through the 1970s and fell in the 1980s,
+ although a formal fixed-versus-drifting comparison requires a separate $Q=0$
+ model.
+
+- *The new observations do not overturn that distinction:* the pandemic
+ produces an extreme volatility episode, while the recent inflation surge does
+ not recreate the persistence of the 1970s.
+
+The reduced-form VAR cannot by itself prove that changes in Federal Reserve
+beliefs caused the coefficient drift, because private behavior and other omitted
+mechanisms can also change reduced-form dynamics.
+
+There is one more twist, and it loops us back to the theory of this section.
+
+The escape-route models of {doc}`phillips_learning` and
+{doc}`phillips_escaping_nash` predict that inflation persistence should *grow*
+along a disinflation as a learning government becomes reluctant to abandon a
+high-inflation self-confirming equilibrium.
+
+The data show the opposite because persistence fell as inflation came down after
+1980.
+
+It helped motivate later learning models, including
+{cite:t}`CogleySargentConquest2005` and {cite:t}`Primiceri2006`, in which
+policymakers' reluctance to disinflate in the 1970s and their eventual
+conversion generate persistence that first rises and then falls.
+
+The friendly debate with Sims, Zha, Bernanke, and Mihov thus did more than
+adjudicate a historical question.
+
+It sharpened the theoretical models of learning and drift that run through this
+section, from the {doc}`self-confirming equilibria `
+of the *Conquest* book to the
+{doc}`drifting Fed beliefs ` used to interpret later
+inflation.
+
+## Exercises
+
+```{exercise}
+:label: csdv_ex1
+
+For an $AR(1)$ process, normalized zero-frequency power is
+
+$$
+g(0)=\frac{1+\rho}{2\pi(1-\rho)}.
+$$
+
+Compute $g(0)$ for $\rho=0$, $0.85$, and $0.97$, and use the persistence
+path above to interpret how inflation dynamics changed around 1980.
+```
+
+```{solution-start} csdv_ex1
+:class: dropdown
+```
+
+```{code-cell} ipython3
+ρ = np.array([0.0, 0.85, 0.97])
+g0 = (1 + ρ) / (2 * np.pi * (1 - ρ))
+pd.Series(g0, index=ρ, name='normalized power at zero').to_frame()
+```
+
+White noise has $g(0)=1/(2\pi)$, while values between roughly 2 and 10
+correspond to highly persistent autoregressions with coefficients between about
+$0.85$ and $0.97$.
+
+Thus, the rise in zero-frequency power during the Great Inflation and its
+decline after 1980 represent a large change in persistence.
+
+```{solution-end}
+```
+
+```{exercise}
+:label: csdv_ex2
+
+Throughout this lecture we noted that a formal contrast between drifting and
+constant coefficients requires refitting the model with $Q=0$, the pure
+"bad luck" special case in which the VAR coefficients are frozen and only the
+stochastic volatilities $H_t$ move.
+
+The sampler already supports this: pass
+`fixed_q=np.zeros((n_coefficients, n_coefficients))` to `run_sampler` to hold
+$Q$ at zero instead of drawing it.
+
+Fit this constant-coefficient model to the Cogley--Sargent sample and plot its
+normalized zero-frequency spectrum $g_{\pi\pi}(0,t)$ against the
+drifting-coefficient path from {numref}`fig-csdv-inflation-persistence`.
+
+What happens to the 1970s rise and post-1980 fall in measured persistence, and
+what does that tell you about whether drifting *volatility* alone can account
+for the persistence dynamics?
+```
+
+```{solution-start} csdv_ex2
+:class: dropdown
+```
+
+With $Q=0$ the elliptical-slice update returns a coefficient path that is
+constant across time, so $A_{t\mid T}$ no longer moves and the only remaining
+source of time variation in $g_{\pi\pi}(0,t)$ is the drifting covariance $R_t$.
+
+```{code-cell} ipython3
+constant_posterior = run_sampler(
+ data['y'],
+ data['x'],
+ prior,
+ n_sweeps=1_000,
+ burn=500,
+ thin=1,
+ seed=42,
+ warmup=200,
+ stable=True,
+ fixed_q=np.zeros((n_coefficients, n_coefficients)),
+)
+
+constant_θ = constant_posterior['SD'].mean(axis=2)
+constant_R = mean_innovation_covariance(
+ constant_posterior['HD'], constant_posterior['CD']
+)
+constant_persistence = np.array([
+ inflation_spectrum(constant_θ[:, t], constant_R[t], zero_frequency)[1][0]
+ for t in range(len(data['dates']))
+])
+
+fig, ax = plt.subplots()
+ax.plot(data['dates'][annual], inflation_persistence[annual], 'o-',
+ lw=2, markersize=3, label='drifting coefficients')
+ax.plot(data['dates'][annual], constant_persistence[annual], 's-',
+ lw=2, markersize=3, label=r'constant coefficients ($Q=0$)')
+ax.set_xlabel('year')
+ax.set_ylabel(r'$g_{\pi\pi}(0,t)$')
+ax.legend()
+plt.show()
+```
+
+With a fixed $A$ the persistence measure still moves, because the *composition*
+of $R_t$ — the relative sizes of the three orthogonal shocks — changes even
+after its overall scale is normalized out.
+
+In fact the constant-coefficient path also climbs to a peak around 1980, as
+inflation innovations grow large relative to the others, so the bad-luck channel
+alone can manufacture much of the *rise*.
+
+What it cannot reproduce is the *fall*: after 1980 the constant-coefficient
+persistence stays elevated, around 2 to 3, for the rest of the sample, whereas
+the drifting-coefficient path collapses back below one.
+
+Freezing $A$ at its full-sample average leaves inflation propagating almost as
+strongly in the 1990s as in the 1970s.
+
+So drifting volatility alone accounts for part of the run-up but none of the
+Volcker-era disinflation of persistence — the post-1980 collapse is evidence
+about the *systematic* dynamics, which is exactly why both channels are needed
+to answer the bad-policy-or-bad-luck question.
+
+```{solution-end}
+```
diff --git a/sync/ledger.yml b/sync/ledger.yml
index 008e779..ae46cf6 100644
--- a/sync/ledger.yml
+++ b/sync/ledger.yml
@@ -355,6 +355,16 @@ lectures/ifp_opi.md:
rewrites:
- from: _admonition/gpu.md
to: _static/_shared/_admonition/gpu.md
+lectures/information_market_equilibrium.md:
+ canonical: intermediate
+ sources:
+ - series: intermediate
+ path: information_market_equilibrium.md
+ digest: 6032176f53b45f763a5f01ca7e06235205c70924d61c1d8fdd729adedda1118e
+ promoted_at:
+ intermediate: 0abf4031c9db8f7a50abb2b21af78929ca39e179
+ assets: []
+ rewrites: []
lectures/inventory_q.md:
canonical: intermediate
interim: true
@@ -591,6 +601,31 @@ lectures/mccall_q.md:
dp-test: 8d26b9349f55f6f185395c82279518bd3ba6a51c
assets: []
rewrites: []
+lectures/networks.md:
+ canonical: intro
+ sources:
+ - series: intro
+ path: networks.md
+ digest: ab39506f1a33f0755a8ddc5b5dc06cbc73f4bd5bee5051d09c1feaba21482df8
+ promoted_at:
+ intro: 8a60bd2728658689b9389e0c0b5ef4e31568d81e
+ assets:
+ - _static/networks/mc.png
+ - _static/networks/poverty_trap_1.png
+ - _static/networks/poverty_trap_2.png
+ - _static/networks/properties.png
+ - _static/networks/weighted.png
+ rewrites:
+ - from: /_static/lecture_specific/networks/poverty_trap_1.png
+ to: _static/networks/poverty_trap_1.png
+ - from: /_static/lecture_specific/networks/poverty_trap_2.png
+ to: _static/networks/poverty_trap_2.png
+ - from: /_static/lecture_specific/networks/properties.png
+ to: _static/networks/properties.png
+ - from: /_static/lecture_specific/networks/weighted.png
+ to: _static/networks/weighted.png
+ - from: /_static/lecture_specific/networks/mc.png
+ to: _static/networks/mc.png
lectures/odu.md:
canonical: dp-test
interim: true
@@ -764,6 +799,16 @@ lectures/perm_income_cons.md:
dp-test: 8d26b9349f55f6f185395c82279518bd3ba6a51c
assets: []
rewrites: []
+lectures/phillips_drifts_volatilities.md:
+ canonical: intermediate
+ sources:
+ - series: intermediate
+ path: phillips_drifts_volatilities.md
+ digest: fe0f8eb5238437d02a201dd2894a20941ce227f8e378e11046835e69bcc3ade5
+ promoted_at:
+ intermediate: 0abf4031c9db8f7a50abb2b21af78929ca39e179
+ assets: []
+ rewrites: []
lectures/rs_inventory_q.md:
canonical: intermediate
interim: true
@@ -1069,6 +1114,36 @@ assets:
path: _static/lecture_specific/lqramsey/firenze.pdf
- series: dp-test
path: _static/lecture_specific/lqramsey/firenze.pdf
+ _static/networks/mc.png:
+ canonical: intro
+ digest: f60ac44c43cf138f7d3483e77681a7bacd99fff90ffbdac4cd49584df1f62020
+ sources:
+ - series: intro
+ path: _static/lecture_specific/networks/mc.png
+ _static/networks/poverty_trap_1.png:
+ canonical: intro
+ digest: 38d83489c1c406c89df60dc745e0dd26a49b3c322c9eaf22d98906034bc0ff11
+ sources:
+ - series: intro
+ path: _static/lecture_specific/networks/poverty_trap_1.png
+ _static/networks/poverty_trap_2.png:
+ canonical: intro
+ digest: 827efb64c396387fe0a3ef5f2545d4e11e20347dafcb445575587d59a913c233
+ sources:
+ - series: intro
+ path: _static/lecture_specific/networks/poverty_trap_2.png
+ _static/networks/properties.png:
+ canonical: intro
+ digest: 29c1499ad1708f38173d26e1d5080b234dc778bbe7e0a42bef2cb78ea1028422
+ sources:
+ - series: intro
+ path: _static/lecture_specific/networks/properties.png
+ _static/networks/weighted.png:
+ canonical: intro
+ digest: f9cdae92a8aafa07b8dfe436c589274058a317ae09cdb306f2a6e6398c9dba5d
+ sources:
+ - series: intro
+ path: _static/lecture_specific/networks/weighted.png
_static/opt_tax_recur/recursive_allocation.py:
canonical: advanced
digest: 0031ff1af560c329dc9011269d8a9b403c4a1afb8004a183a6a24fe3f0a41e59
+
+