From 7f70be1976c41f802de5c0f4b404c7ec0bd36b35 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 17:07:26 +1000 Subject: [PATCH 1/3] =?UTF-8?q?=F0=9F=94=84=20resync=20aiyagari.md?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .translate/state/aiyagari.md.yml | 6 + lectures/aiyagari.md | 867 +++++++++++++++++++++++-------- 2 files changed, 662 insertions(+), 211 deletions(-) create mode 100644 .translate/state/aiyagari.md.yml diff --git a/.translate/state/aiyagari.md.yml b/.translate/state/aiyagari.md.yml new file mode 100644 index 00000000..3112f1be --- /dev/null +++ b/.translate/state/aiyagari.md.yml @@ -0,0 +1,6 @@ +source-sha: 0bfcac8105ac8c11f06797579a9f4d565db7b53d +synced-at: "2026-07-18" +model: claude-sonnet-5 +mode: RESYNC +section-count: 4 +tool-version: 0.17.0 diff --git a/lectures/aiyagari.md b/lectures/aiyagari.md index bbe5db70..75ec7d02 100644 --- a/lectures/aiyagari.md +++ b/lectures/aiyagari.md @@ -7,6 +7,22 @@ kernelspec: display_name: Python 3 language: python name: python3 +translation: + title: 艾亚加里模型 + headings: + Overview: 概述 + Overview::Preliminaries: 预备知识 + Overview::References: 参考文献 + The Economy: 经济模型 + The Economy::Households: 家庭 + The Economy::Firms: 企业 + The Economy::Equilibrium: 均衡 + Implementation: 代码实现 + Implementation::Primitives and operators: 原语与算子 + Implementation::Capital supply: 资本供给 + Implementation::Equilibrium: 均衡 + Implementation::Supply and demand curves: 供给和需求曲线 + Exercises: 练习 --- (aiyagari)= @@ -20,17 +36,19 @@ kernelspec: # 艾亚加里模型 +```{include} _admonition/gpu.md +``` + ```{contents} 目录 :depth: 2 ``` -除了Anaconda中包含的包之外，本讲座还需要以下库： +除了Anaconda中包含的包之外，我们还需要安装JAX -```{code-cell} ipython ---- -tags: [hide-output] ---- -!pip install quantecon +```{code-cell} ipython3 +:tags: [hide-output] + +!pip install quantecon jax ``` ## 概述 @@ -42,7 +60,7 @@ tags: [hide-output] 该模型具有以下特点： * 异质性主体 -* 单一的借贷工具 +* 单一的外生借贷工具 * 对个人主体借款额度的限制 艾亚加里模型已被用于研究多个主题，包括： @@ -52,19 +70,35 @@ tags: [hide-output] * 财富分布的形状 {cite}`benhabib2015` * 等等 -让我们从导入必要的包开始： +### 预备知识 + +我们使用以下导入： -```{code-cell} ipython +```{code-cell} ipython3 +import quantecon as qe import matplotlib.pyplot as plt -import matplotlib as mpl -FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf" -mpl.font_manager.fontManager.addfont(FONTPATH) -plt.rcParams['font.family'] = ['Source Han Serif SC'] +import jax +import jax.numpy as jnp +from typing import NamedTuple +from scipy.optimize import bisect +``` -plt.rcParams["figure.figsize"] = (11, 5) #设置默认图形大小 -import numpy as np -from quantecon.markov import DiscreteDP -from numba import jit +我们将在JAX中使用64位浮点数以提高精度。 + +```{code-cell} ipython3 +jax.config.update("jax_enable_x64", True) +``` + +我们将使用以下函数来计算随机矩阵的平稳分布（关于该算法的参考资料，请参见[Economic Dynamics](https://johnstachurski.net/edtc)第88页）。 + +```{code-cell} ipython3 +@jax.jit +def compute_stationary(P): + n = P.shape[0] + I = jnp.identity(n) + O = jnp.ones((n, n)) + A = I - jnp.transpose(P) + O + return jnp.linalg.solve(A, jnp.ones(n)) ``` ### 参考文献 @@ -103,7 +137,7 @@ $$ * $c_t$ 是当前消费 * $a_t$ 是资产 -* $z_t$ 是劳动收入的随机组成部分，捕捉了失业风险等 +* $z_t$ 是劳动收入的外生组成部分，捕捉了随机失业风险等 * $w$ 是工资率 * $r$ 是净利率 * $B$ 是主体允许借入的最大金额 @@ -114,7 +148,7 @@ $$ 在这个模型的简单版本中，家庭无弹性地供给劳动，因为他们不重视闲暇。 -## 企业 +### 企业 企业通过雇佣资本和劳动来生产产出。 @@ -127,23 +161,33 @@ $$ 企业的产出为 $$ -Y_t = A K_t^{\alpha} N^{1 - \alpha} +Y = A K^{\alpha} N^{1 - \alpha} $$ 其中： * $A$ 和 $\alpha$ 是参数，$A > 0$ 且 $\alpha \in (0, 1)$ -* $K_t$ 是总资本 +* $K$ 是总资本 * $N$ 是总劳动供给（在这个简单版本的模型中保持不变） 企业的问题是 $$ -max_{K, N} \left\{ A K_t^{\alpha} N^{1 - \alpha} - (r + \delta) K - w N \right\} +\max_{K, N} \left\{ A K^{\alpha} N^{1 - \alpha} - (r + \delta) K - w N \right\} $$ 参数 $\delta$ 是折旧率。 +这些参数存储在以下namedtuple中： + +```{code-cell} ipython3 +class Firm(NamedTuple): + A: float = 1.0 # 全要素生产率 + N: float = 1.0 # 总劳动供给 + α: float = 0.33 # 资本份额 + δ: float = 0.05 # 折旧率 +``` + 从关于资本的一阶条件，企业的资本需求反函数为 ```{math} @@ -152,6 +196,15 @@ $$ r = A \alpha \left( \frac{N}{K} \right)^{1 - \alpha} - \delta ``` +```{code-cell} ipython3 +def r_given_k(K, firm): + """ + 资本的需求反函数。与给定资本需求K相关的利率。 + """ + A, N, α, δ = firm + return A * α * (N / K)**(1 - α) - δ +``` + 使用这个表达式和企业关于劳动的一阶条件，我们可以将均衡工资率表示为 $r$ 的函数： ```{math} @@ -160,9 +213,18 @@ r = A \alpha \left( \frac{N}{K} \right)^{1 - \alpha} - \delta w(r) = A (1 - \alpha) (A \alpha / (r + \delta))^{\alpha / (1 - \alpha)} ``` +```{code-cell} ipython3 +def r_to_w(r, firm): + """ + 与给定利率r相关的均衡工资。 + """ + A, N, α, δ = firm + return A * (1 - α) * (A * α / (r + δ))**(α / (1 - α)) +``` + ### 均衡 -我们构建一个*平稳理性预期均衡*（SREE）。 +我们构建一个**平稳理性预期均衡**（SREE）。 在这样的均衡中： @@ -171,267 +233,650 @@ w(r) = A (1 - \alpha) (A \alpha / (r + \delta))^{\alpha / (1 - \alpha)} 更详细地说，SREE列出了价格、储蓄和生产策略的集合，使得： -* 家庭在给定价格下选择指定的储蓄策略 +* 家庭在给定价格下希望选择指定的储蓄策略 * 企业在相同价格下最大化利润 * 产生的总量与价格一致；特别是，资本需求等于供给 * 总量（定义为横截面平均值）保持不变 -在实践中，一旦参数值设定，我们可以通过以下步骤检查SREE： +## 代码实现 -1. 选择一个提议的总资本量 $K$ -2. 确定相应的价格，其中利率 $r$ 由 {eq}`aiy_rgk` 决定，工资率 $w(r)$ 由 {eq}`aiy_wgr` 给出 -3. 确定给定这些价格下家庭的共同最优储蓄策略 -4. 计算给定这个储蓄策略下的稳态资本平均值 +让我们看看如何在实践中计算这样的均衡。 -如果最终数量与 $K$ 一致，那么我们就得到了一个SREE。 +下面我们提供代码来求解家庭问题，将 $r$ 和 $w$ 视为固定值。 -## 代码 +### 原语与算子 -让我们看看如何在实践中计算这样的均衡。 +我们将使用价值函数迭代来求解家庭问题。 -为了求解家庭的动态规划问题，我们将使用来自 [QuantEcon.py](https://quantecon.org/quantecon-py/) 的 [DiscreteDP](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/markov/ddp.py) 类。 +首先我们设置一个 `NamedTuple` 来存储定义家庭资产积累问题的参数，以及用于求解的网格 -我们的第一个任务是最不令人兴奋的：编写代码，将家庭问题的参数映射到生成 `DiscreteDP` 实例所需的 `R` 和 `Q` 矩阵。 +```{code-cell} ipython3 +class Household(NamedTuple): + β: float # 贴现因子 + a_grid: jnp.ndarray # 资产网格 + z_grid: jnp.ndarray # 外生状态 + Π: jnp.ndarray # 转移矩阵 -下面是一段样板代码，它完成了这个任务。 +def create_household(β=0.96, # 贴现因子 + Π=[[0.9, 0.1], [0.1, 0.9]], # 马尔可夫链 + z_grid=[0.1, 1.0], # 外生状态 + a_min=1e-10, a_max=12.5, # 资产网格 + a_size=100): + """ + 使用自定义网格创建Household namedtuple。 + """ + a_grid = jnp.linspace(a_min, a_max, a_size) + z_grid, Π = map(jnp.array, (z_grid, Π)) + return Household(β=β, a_grid=a_grid, z_grid=z_grid, Π=Π) +``` + +现在我们假设 $u(c) = \log(c)$ + +```{code-cell} ipython3 +u = jnp.log +``` + +这是一个存储工资率和利率（带默认值）的namedtuple -在阅读代码时，以下信息会很有帮助： +```{code-cell} ipython3 +class Prices(NamedTuple): + r: float = 0.01 # 利率 + w: float = 1.0 # 工资 +``` -* `R` 需要是一个矩阵，其中 `R[s, a]` 是在状态 `s` 下采取行动 `a` 的回报。 -* `Q` 需要是一个三维数组，其中 `Q[s, a, s']` 是在当前状态为 `s` 且当前行动为 `a` 时转移到状态 `s'` 的概率。 +现在我们建立贝尔曼方程右侧（最大化之前）的向量化版本，它是一个表示以下内容的三维数组 -（关于 `DiscreteDP` 的更详细讨论可以在 [Advanced Quantitative Economics with Python](https://python-advanced.quantecon.org) 讲座系列的 [Discrete State Dynamic Programming](https://python-advanced.quantecon.org/discrete_dp.html) 讲座中找到。） +$$ +B(a, z, a') = u(wz + (1+r)a - a') + \beta \sum_{z'} v(a', z') \Pi(z, z') +$$ -这里我们将状态设为 $s_t := (a_t, z_t)$，其中 $a_t$ 是资产，$z_t$ 是冲击。 +对所有 $(a, z, a')$ 成立。 -行动是选择下一期的资产水平 $a_{t+1}$。 +```{code-cell} ipython3 +def B(v, household, prices): + # 解包 + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) + r, w = prices -我们使用Numba来加速循环，这样当参数改变时我们可以高效地更新矩阵。 + # 将当前消费计算为数组 c[i, j, ip] + a = jnp.reshape(a_grid, (a_size, 1, 1)) # a[i] -> a[i, j, ip] + z = jnp.reshape(z_grid, (1, z_size, 1)) # z[j] -> z[i, j, ip] + ap = jnp.reshape(a_grid, (1, 1, a_size)) # ap[ip] -> ap[i, j, ip] + c = w * z + (1 + r) * a - ap -该类还包括一组默认参数，除非另有说明，否则我们将采用这些参数。 + # 计算(a, z, ap)所有组合的延续回报 + v = jnp.reshape(v, (1, 1, a_size, z_size)) # v[ip, jp] -> v[i, j, ip, jp] + Π = jnp.reshape(Π, (1, z_size, 1, z_size)) # Π[j, jp] -> Π[i, j, ip, jp] + EV = jnp.sum(v * Π, axis=-1) # 对最后一个索引jp求和 -```{code-cell} python3 -class Household: + # 计算贝尔曼方程的右侧 + return jnp.where(c > 0, u(c) + β * EV, -jnp.inf) +``` + +下一个函数计算贪婪策略 + +```{code-cell} ipython3 +def get_greedy(v, household, prices): + """ + 计算v-贪婪策略σ，以一组索引的形式返回。如果 + σ[i, j]等于ip，则a_grid[ip]是i, j处的最大化元素。 """ - 这个类接收定义家庭资产积累问题的参数，并计算相应的回报和转移矩阵R - 和Q，这些矩阵是生成DiscreteDP实例所必需的，从而求解最优策略。 + # 对ap求argmax + return jnp.argmax(B(v, household, prices), axis=-1) +``` - 关于索引的说明：我们需要将状态空间S枚举为序列S = {0, ..., n}。 - 为此，(a_i, z_i)索引对根据以下规则映射到s_i索引： +我们定义贝尔曼算子 $T$，它接受一个价值函数 $v$ 并返回贝尔曼方程给出的 $Tv$ - s_i = a_i * z_size + z_i +```{code-cell} ipython3 +def T(v, household, prices): + """ + 贝尔曼算子。接受一个价值函数v并返回Tv。 + """ + return jnp.max(B(v, household, prices), axis=-1) +``` - 要反转这个映射，使用： +这是价值函数迭代，它反复应用贝尔曼算子直至收敛 - a_i = s_i // z_size (整数除法) - z_i = s_i % z_size +```{code-cell} ipython3 +@jax.jit +def value_function_iteration(household, prices, tol=1e-4, max_iter=10_000): """ + 使用编译的JAX循环实现价值函数迭代。 + """ + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) + def condition_function(loop_state): + i, v, error = loop_state + return jnp.logical_and(error > tol, i < max_iter) - def __init__(self, - r=0.01, # 利率 - w=1.0, # 工资 - β=0.96, # 贴现因子 - a_min=1e-10, - Π=[[0.9, 0.1], [0.1, 0.9]], # 马尔可夫链 - z_vals=[0.1, 1.0], # 外生状态 - a_max=18, - a_size=200): - - # 存储值，设置a和z的网格 - self.r, self.w, self.β = r, w, β - self.a_min, self.a_max, self.a_size = a_min, a_max, a_size - - self.Π = np.asarray(Π) - self.z_vals = np.asarray(z_vals) - self.z_size = len(z_vals) - - self.a_vals = np.linspace(a_min, a_max, a_size) - self.n = a_size * self.z_size - - # 构建数组Q - self.Q = np.zeros((self.n, a_size, self.n)) - self.build_Q() - - # 构建数组R - self.R = np.empty((self.n, a_size)) - self.build_R() - - def set_prices(self, r, w): - """ - 使用此方法重置价格。调用此方法将触发R的重新构建。 - """ - self.r, self.w = r, w - self.build_R() - - def build_Q(self): - populate_Q(self.Q, self.a_size, self.z_size, self.Π) - - def build_R(self): - self.R.fill(-np.inf) - populate_R(self.R, - self.a_size, - self.z_size, - self.a_vals, - self.z_vals, - self.r, - self.w) - - -# 使用JIT编译的函数进行繁重工作 - -@jit -def populate_R(R, a_size, z_size, a_vals, z_vals, r, w): - n = a_size * z_size - for s_i in range(n): - a_i = s_i // z_size - z_i = s_i % z_size - a = a_vals[a_i] - z = z_vals[z_i] - for new_a_i in range(a_size): - a_new = a_vals[new_a_i] - c = w * z + (1 + r) * a - a_new - if c > 0: - R[s_i, new_a_i] = np.log(c) # 效用 - -@jit -def populate_Q(Q, a_size, z_size, Π): - n = a_size * z_size - for s_i in range(n): - z_i = s_i % z_size - for a_i in range(a_size): - for next_z_i in range(z_size): - Q[s_i, a_i, a_i*z_size + next_z_i] = Π[z_i, next_z_i] + def update(loop_state): + i, v, error = loop_state + v_new = T(v, household, prices) + error = jnp.max(jnp.abs(v_new - v)) + return i + 1, v_new, error + # 初始循环状态 + v_init = jnp.zeros((a_size, z_size)) + loop_state_init = (0, v_init, tol + 1) -@jit -def asset_marginal(s_probs, a_size, z_size): - a_probs = np.zeros(a_size) - for a_i in range(a_size): - for z_i in range(z_size): - a_probs[a_i] += s_probs[a_i*z_size + z_i] - return a_probs -``` + # 运行不动点迭代 + i, v, error = jax.lax.while_loop(condition_function, update, loop_state_init) -作为第一个例子，让我们计算并绘制在固定价格下的最优积累策略。 + return get_greedy(v, household, prices) +``` -```{code-cell} python3 -# 示例价格 -r = 0.03 -w = 0.956 +作为我们能做的第一个例子，让我们计算并绘制在固定价格下的最优积累策略 +```{code-cell} ipython3 # 创建Household实例 -am = Household(a_max=20, r=r, w=w) +household = create_household() +prices = Prices() -# 使用实例构建离散动态规划 -am_ddp = DiscreteDP(am.R, am.Q, am.β) +r, w = prices +print(f"Interest rate: {r}, Wage: {w}") +``` -# 使用策略函数迭代求解 -results = am_ddp.solve(method='policy_iteration') +```{code-cell} ipython3 +with qe.Timer(): + σ_star = value_function_iteration(household, prices).block_until_ready() +``` -# 简化名称 -z_size, a_size = am.z_size, am.a_size -z_vals, a_vals = am.z_vals, am.a_vals -n = a_size * z_size +下图显示了在不同外生状态值下的资产积累策略 -# 获取所有最优行动，在每行中固定z值 -a_star = np.empty((z_size, a_size)) -for s_i in range(n): - a_i = s_i // z_size - z_i = s_i % z_size - a_star[z_i, a_i] = a_vals[results.sigma[s_i]] +```{code-cell} ipython3 +β, a_grid, z_grid, Π = household fig, ax = plt.subplots(figsize=(9, 9)) -ax.plot(a_vals, a_vals, 'k--') # 45度线 -for i in range(z_size): - lb = f'$z = {z_vals[i]:.2}$' - ax.plot(a_vals, a_star[i, :], lw=2, alpha=0.6, label=lb) +ax.plot(a_grid, a_grid, 'k--', label="45度线") +for j, z in enumerate(z_grid): + lb = f'$z = {z:.2}$' + policy_vals = a_grid[σ_star[:, j]] + ax.plot(a_grid, policy_vals, lw=2, alpha=0.6, label=lb) ax.set_xlabel('当前资产') ax.set_ylabel('下一期资产') ax.legend(loc='upper left') - plt.show() ``` 该图显示了在不同外生状态值下的资产积累策略。 -现在我们要计算均衡。 +### 资本供给 -让我们先通过可视化来做到这一点。 +要开始考虑均衡，我们需要知道在给定利率 $r$ 下家庭供给多少资本。 -以下代码绘制了总供给和需求曲线。 +该数量可以通过取最优策略下资产的平稳分布并计算其均值来计算。 -交点给出了均衡利率和资本。 +下一个函数通过以下步骤计算给定策略 $\sigma$ 的平稳分布： -```{code-cell} python3 -A = 1.0 -N = 1.0 -α = 0.33 -β = 0.96 -δ = 0.05 +* 计算 $P_{\sigma}$ 的平稳分布 $\psi = (\psi(a, z))$，它定义了策略 $\sigma$ 下状态 $(a_t, z_t)$ 的马尔可夫链。 +* 对 $z_t$ 求和以得到 $a_t$ 的边际分布。 +```{code-cell} ipython3 +@jax.jit +def compute_asset_stationary(σ, household): + # 解包 + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) -def r_to_w(r): - """ - 与给定利率r相关的均衡工资。 - """ - return A * (1 - α) * (A * α / (r + δ))**(α / (1 - α)) + # 将P_σ构建为形式为P_σ[i, j, ip, jp]的数组 + ap_idx = jnp.arange(a_size) + ap_idx = jnp.reshape(ap_idx, (1, 1, a_size, 1)) + σ = jnp.reshape(σ, (a_size, z_size, 1, 1)) + A = jnp.where(σ == ap_idx, 1, 0) + Π = jnp.reshape(Π, (1, z_size, 1, z_size)) + P_σ = A * Π + + # 将P_σ重塑为矩阵 + n = a_size * z_size + P_σ = jnp.reshape(P_σ, (n, n)) + + # 获取平稳分布并重塑回[i, j]网格 + ψ = compute_stationary(P_σ) + ψ = jnp.reshape(ψ, (a_size, z_size)) + + # 沿行求和以得到资产的边际分布 + ψ_a = jnp.sum(ψ, axis=1) + return ψ_a +``` + +让我们试运行一下。 + +```{code-cell} ipython3 +ψ_a = compute_asset_stationary(σ_star, household) + +fig, ax = plt.subplots() +ax.bar(household.a_grid, ψ_a) +ax.set_xlabel("资产水平") +ax.set_ylabel("概率质量") +plt.show() +``` + +该分布应该总和为一： -def rd(K): +```{code-cell} ipython3 +ψ_a.sum() +``` + +下一个函数计算给定工资和利率下，策略 $\sigma$ 下家庭的总资本供给 + +```{code-cell} ipython3 +def capital_supply(σ, household): """ - 资本的需求反函数。与给定资本需求K相关的利率。 + 在给定r和w的情况下，策略下诱致的资本存量水平。 """ - return A * α * (N / K)**(1 - α) - δ + β, a_grid, z_grid, Π = household + ψ_a = compute_asset_stationary(σ, household) + return float(jnp.sum(ψ_a * a_grid)) +``` +### 均衡 -def prices_to_capital_stock(am, r): - """ - 将价格映射到诱致的资本存量水平。 +我们通过以下方式计算SREE： - 参数: - ---------- +1. 设 $n=0$ 并以总资本的初始猜测值 $K_0$ 开始。 +1. 给定 $K_n$，从企业的决策问题中确定价格 $r, w$。 +1. 计算给定这些价格下家庭的最优储蓄策略。 +1. 将总资本 $K_{n+1}$ 计算为给定该储蓄策略下稳态资本的均值。 +1. 如果 $K_{n+1} \approx K_n$，停止；否则转到步骤2。 - am : Household - aiyagari_household.Household的实例 - r : float - 利率 - """ - w = r_to_w(r) - am.set_prices(r, w) - aiyagari_ddp = DiscreteDP(am.R, am.Q, β) - # 计算最优策略 - results = aiyagari_ddp.solve(method='policy_iteration') - # 计算稳态分布 - stationary_probs = results.mc.stationary_distributions[0] - # 提取资产的边际分布 - asset_probs = asset_marginal(stationary_probs, am.a_size, am.z_size) - # 返回K - return np.sum(asset_probs * am.a_vals) +我们可以将步骤2-4中的操作序列写为 +$$ +K_{n + 1} = G(K_n) +$$ -# 创建Household实例 -am = Household(a_max=20) +如果 $K_{n+1}$ 与 $K_n$ 一致，那么我们就得到了一个SREE。 -# 使用实例构建离散动态规划 -am_ddp = DiscreteDP(am.R, am.Q, am.β) +换句话说，我们的问题是找到一维映射 $G$ 的不动点。 + +以下是用Python函数表示的 $G$ + +```{code-cell} ipython3 +def G(K, firm, household): + # 获取与K相关的价格r, w + r = r_given_k(K, firm) + w = r_to_w(r, firm) + + # 用这些价格生成一个household对象，计算 + # 总资本。 + prices = Prices(r=r, w=w) + σ_star = value_function_iteration(household, prices) + return capital_supply(σ_star, household) +``` + +作为第一步，让我们直观地检查一下 + +```{code-cell} ipython3 +num_points = 50 +firm = Firm() +household = create_household() +k_vals = jnp.linspace(4, 12, num_points) +out = [G(k, firm, household) for k in k_vals] + +fig, ax = plt.subplots(figsize=(11, 8)) +ax.plot(k_vals, out, lw=2, alpha=0.6, label='$G$') +ax.plot(k_vals, k_vals, 'k--', label="45度线") +ax.set_xlabel('资本') +ax.legend() +plt.show() +``` + +现在让我们来计算均衡。 + +看上图，我们发现简单的迭代方案 $K_{n+1} = G(K_n)$ 会在高低值之间循环，导致收敛缓慢。 + +因此，我们使用如下形式的阻尼迭代方案 + +$$ +K_{n+1} = \alpha K_n + (1-\alpha) G(K_n) +$$ + +```{code-cell} ipython3 +def compute_equilibrium(firm, household, + K0=6, α=0.99, max_iter=1_000, tol=1e-4, + print_skip=10, verbose=False): + n = 0 + K = K0 + error = tol + 1 + while error > tol and n < max_iter: + new_K = α * K + (1 - α) * G(K, firm, household) + error = abs(new_K - K) + K = new_K + n += 1 + if verbose and n % print_skip == 0: + print(f"At iteration {n} with error {error}") + return K, n +``` + +```{code-cell} ipython3 +firm = Firm() +household = create_household() +print("\nComputing equilibrium capital stock") +with qe.Timer(): + K_star, n = compute_equilibrium(firm, household, K0=6.0) +print(f"Computed equilibrium {K_star:.5} in {n} iterations") +``` + +考虑到我们可以多快地求解家庭问题，这种收敛速度并不算快。 + +你可以尝试改变 $\alpha$，但通常这个参数很难事先设定。 + +在下面的练习中，你将被要求改用二分法，这种方法通常表现更好。 + +### 供给和需求曲线 + +我们可以使用供给和需求曲线来可视化均衡。 + +以下代码绘制了总供给和需求曲线。 + +交点给出了均衡利率和资本 + +```{code-cell} ipython3 +def prices_to_capital_stock(household, r, firm): + """ + 将价格映射到诱致的资本存量水平。 + """ + w = r_to_w(r, firm) + prices = Prices(r=r, w=w) + + # 计算最优策略 + σ_star = value_function_iteration(household, prices) + + # 计算资本供给 + return capital_supply(σ_star, household) # 创建计算资本需求和供给的r值网格 num_points = 20 -r_vals = np.linspace(0.005, 0.04, num_points) +r_vals = jnp.linspace(0.005, 0.04, num_points) # 计算资本供给 -k_vals = np.empty(num_points) -for i, r in enumerate(r_vals): - k_vals[i] = prices_to_capital_stock(am, r) +k_vals = [] +for r in r_vals: + k_vals.append(prices_to_capital_stock(household, r, firm)) # 绘制与企业的资本需求相对 fig, ax = plt.subplots(figsize=(11, 8)) -ax.plot(k_vals, r_vals, lw=2, alpha=0.6, label='资本供给') -ax.plot(k_vals, rd(k_vals), lw=2, alpha=0.6, label='资本需求') -ax.grid() +ax.plot(k_vals, r_vals, lw=2, alpha=0.6, + label='资本供给') +ax.plot(k_vals, r_given_k( + jnp.array(k_vals), firm), lw=2, alpha=0.6, + label='资本需求') + +# 在均衡点添加标记 +r_star = r_given_k(K_star, firm) +ax.plot(K_star, r_star, 'o', markersize=10, label='均衡') + ax.set_xlabel('资本') ax.set_ylabel('利率') ax.legend(loc='upper right') plt.show() -``` \ No newline at end of file +``` + +## 练习 + +```{exercise} +:label: aiyagari_ex1 + +编写一个新版本的 `compute_equilibrium`，使用 `scipy.optimize` 中的 `bisect` 而不是阻尼迭代。 + +看看你能否使它比之前的版本更快。 + +在 `bisect` 中， + +* 你应该设置 `xtol=1e-4`，以获得与之前版本相同的误差容限。 +* 对于二分法程序的上下界，尝试使用 `a = 1.0` 和 `b = 20.0`。 +``` + +```{solution-start} aiyagari_ex1 +:class: dropdown +``` + +我们使用二分法找到函数 $h(k) = k - G(k)$ 的零点 + +```{code-cell} ipython3 +def compute_equilibrium_bisect(firm, household, a=1.0, b=20.0): + K = bisect(lambda k: k - G(k, firm, household), a, b, xtol=1e-4) + return K + +firm = Firm() +household = create_household() +print("\nComputing equilibrium capital stock using bisection") +with qe.Timer(): + K_star = compute_equilibrium_bisect(firm, household) +print(f"Computed equilibrium capital stock {K_star:.5}") +``` + +二分法比阻尼迭代方案更快。 + +```{solution-end} +``` + +```{exercise-start} +:label: aiyagari_ex2 +``` + +展示均衡资本存量如何随 $\beta$ 变化。 + +使用以下 $\beta$ 值并绘制你发现的关系。 + +```{code-cell} ipython3 +:tags: [hide-output] + +β_vals = jnp.linspace(0.94, 0.98, 20) +``` + +```{exercise-end} +``` + +```{solution-start} aiyagari_ex2 +:class: dropdown +``` + +```{code-cell} ipython3 +K_vals = [] +K = 6.0 # 初始猜测值 + +for β in β_vals: + household = create_household(β=β) + K = compute_equilibrium_bisect(firm, household, 0.5 * K, 1.5 * K) + print(f"Computed equilibrium {K:.4} at β = {β}") + K_vals.append(K) + +fig, ax = plt.subplots() +ax.plot(β_vals, K_vals, ms=2) +ax.set_xlabel(r'$\beta$') +ax.set_ylabel('资本') +plt.show() +``` + +```{solution-end} +``` + +```{exercise-start} +:label: aiyagari_ex3 +``` + +在本讲座中，我们使用价值函数迭代来求解家庭问题。 + +另一种方法是霍华德策略迭代（HPI），在[Dynamic Programming](https://dp.quantecon.org/)中有详细讨论。 + +对于某些问题，HPI可以比VFI更快，因为它使用更少但计算量更大的迭代。 + +你的任务是实现霍华德策略迭代，并将结果与价值函数迭代进行比较。 + +**你需要的关键概念：** + +霍华德策略迭代需要计算策略 $\sigma$ 的价值 $v_{\sigma}$，定义为： + +$$ +v_{\sigma} = (I - \beta P_{\sigma})^{-1} r_{\sigma} +$$ + +其中 $r_{\sigma}$ 是策略 $\sigma$ 下的回报向量，$P_{\sigma}$ 是由 $\sigma$ 诱导的转移矩阵。 + +要解决这个问题，你需要： +1. 计算当前回报 $r_{\sigma}(a, z) = u((1 + r)a + wz - \sigma(a, z))$ +2. 建立线性算子 $R_{\sigma}$，其中 $(R_{\sigma} v)(a, z) = v(a, z) - \beta \sum_{z'} v(\sigma(a, z), z') \Pi(z, z')$ +3. 使用 `jax.scipy.sparse.linalg.bicgstab` 求解 $v_{\sigma} = R_{\sigma}^{-1} r_{\sigma}$ + +你可以使用本讲座中已经定义的 `get_greedy` 函数。 + +实现以下霍华德策略迭代程序： + +```python +def howard_policy_iteration(household, prices, + tol=1e-4, max_iter=10_000, verbose=False): + """ + 霍华德策略迭代程序。 + """ + # 你的代码在这里 + pass +``` + +实现后，使用HPI计算均衡资本存量，并验证在默认参数值下它是否产生与VFI大致相同的结果。 + +```{exercise-end} +``` + +```{solution-start} aiyagari_ex3 +:class: dropdown +``` + +首先，我们需要为霍华德策略迭代实现辅助函数。 + +以下函数计算数组 $r_{\sigma}$，它给出策略 $\sigma$ 下的当前回报： + +```{code-cell} ipython3 +def compute_r_σ(σ, household, prices): + """ + 计算策略σ下每个i, j处的当前回报。特别地， + + r_σ[i, j] = u((1 + r)a[i] + wz[j] - a'[ip]) + + 当 ip = σ[i, j] 时。 + """ + # 解包 + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) + r, w = prices + + # 计算 r_σ[i, j] + a = jnp.reshape(a_grid, (a_size, 1)) + z = jnp.reshape(z_grid, (1, z_size)) + ap = a_grid[σ] + c = (1 + r) * a + w * z - ap + r_σ = u(c) + + return r_σ +``` + +线性算子 $R_{\sigma}$ 定义为： + +```{code-cell} ipython3 +def R_σ(v, σ, household): + # 解包 + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) + + # 建立数组 v[σ[i, j], jp] + zp_idx = jnp.arange(z_size) + zp_idx = jnp.reshape(zp_idx, (1, 1, z_size)) + σ = jnp.reshape(σ, (a_size, z_size, 1)) + V = v[σ, zp_idx] + + # 将 Π[j, jp] 扩展为 Π[i, j, jp] + Π = jnp.reshape(Π, (1, z_size, z_size)) + + # 计算并返回 v[i, j] - β Σ_jp v[σ[i, j], jp] * Π[j, jp] + return v - β * jnp.sum(V * Π, axis=-1) +``` + +下一个函数计算给定策略的终身价值： + +```{code-cell} ipython3 +def get_value(σ, household, prices): + """ + 通过计算以下内容获得策略σ的终身价值 + + v_σ = R_σ^{-1} r_σ + """ + r_σ = compute_r_σ(σ, household, prices) + + # 将 R_σ 简化为关于v的函数 + _R_σ = lambda v: R_σ(v, σ, household) + + # 使用迭代程序计算 v_σ = R_σ^{-1} r_σ。 + return jax.scipy.sparse.linalg.bicgstab(_R_σ, r_σ)[0] +``` + +现在我们可以实现霍华德策略迭代： + +```{code-cell} ipython3 +@jax.jit +def howard_policy_iteration(household, prices, tol=1e-4, max_iter=10_000): + """ + 使用编译的JAX循环实现霍华德策略迭代程序。 + """ + β, a_grid, z_grid, Π = household + a_size, z_size = len(a_grid), len(z_grid) + + def condition_function(loop_state): + i, σ, v_σ, error = loop_state + return jnp.logical_and(error > tol, i < max_iter) + + def update(loop_state): + i, σ, v_σ, error = loop_state + σ_new = get_greedy(v_σ, household, prices) + v_σ_new = get_value(σ_new, household, prices) + error = jnp.max(jnp.abs(v_σ_new - v_σ)) + return i + 1, σ_new, v_σ_new, error + + # 初始循环状态 + σ_init = jnp.zeros((a_size, z_size), dtype=int) + v_σ_init = get_value(σ_init, household, prices) + loop_state_init = (0, σ_init, v_σ_init, tol + 1) + + # 运行不动点迭代 + i, σ, v_σ, error = jax.lax.while_loop(condition_function, update, loop_state_init) + + return σ +``` + +现在让我们创建一个使用HPI的G函数的修改版本： + +```{code-cell} ipython3 +def G_hpi(K, firm, household): + # 获取与K相关的价格r, w + r = r_given_k(K, firm) + w = r_to_w(r, firm) + + # 生成价格并使用HPI计算总资本。 + prices = Prices(r=r, w=w) + σ_star = howard_policy_iteration(household, prices) + return capital_supply(σ_star, household) +``` + +并使用HPI计算均衡： + +```{code-cell} ipython3 +def compute_equilibrium_bisect_hpi(firm, household, a=1.0, b=20.0): + K = bisect(lambda k: k - G_hpi(k, firm, household), a, b, xtol=1e-4) + return K + +firm = Firm() +household = create_household() +print("\nComputing equilibrium capital stock using HPI") +with qe.Timer(): + K_star_hpi = compute_equilibrium_bisect_hpi(firm, household) +print(f"Computed equilibrium capital stock with HPI: {K_star_hpi:.5}") +print(f"Previous equilibrium capital stock with VFI: {K_star:.5}") +print(f"Difference: {abs(K_star_hpi - K_star):.6}") +``` + +结果显示两种方法产生了大致相同的均衡，证实了HPI是VFI的有效替代方案。 + +```{solution-end} +``` \ No newline at end of file From 87c90987f4286bfc651222c6ced525a09a9de71a Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 21:20:39 +1000 Subject: [PATCH 2/3] Restore CJK font configuration dropped by the resync The resynced file keeps Chinese plot labels but reverted the import cell to the source's exact form, losing the Source Han Serif setup - Chinese in figures would render as missing glyphs. Residual #107 class found in the merge review; applied wave-wide. Co-Authored-By: Claude Fable 5 --- lectures/aiyagari.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lectures/aiyagari.md b/lectures/aiyagari.md index 75ec7d02..25bf0cc9 100644 --- a/lectures/aiyagari.md +++ b/lectures/aiyagari.md @@ -77,6 +77,10 @@ translation: ```{code-cell} ipython3 import quantecon as qe import matplotlib.pyplot as plt +import matplotlib as mpl +FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf" +mpl.font_manager.fontManager.addfont(FONTPATH) +plt.rcParams['font.family'] = ['Source Han Serif SC'] import jax import jax.numpy as jnp from typing import NamedTuple From 4edc693b245342f047600c4be043ffad19548c00 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sun, 19 Jul 2026 09:58:10 +1000 Subject: [PATCH 3/3] Restore trailing newline (engine issue tracked in QuantEcon/action-translation#116) Co-Authored-By: Claude Fable 5 --- lectures/aiyagari.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/aiyagari.md b/lectures/aiyagari.md index 25bf0cc9..f44e5c91 100644 --- a/lectures/aiyagari.md +++ b/lectures/aiyagari.md @@ -883,4 +883,4 @@ print(f"Difference: {abs(K_star_hpi - K_star):.6}") 结果显示两种方法产生了大致相同的均衡，证实了HPI是VFI的有效替代方案。 ```{solution-end} -``` \ No newline at end of file +```