From 470ef8135720433b87b99cf2e8ecef2df2e94891 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 17:38:31 +1000 Subject: [PATCH 1/4] =?UTF-8?q?=F0=9F=94=84=20resync=20ifp=5Fadvanced.md?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .translate/state/ifp_advanced.md.yml | 6 + lectures/ifp_advanced.md | 797 +++++++++++++++++---------- 2 files changed, 511 insertions(+), 292 deletions(-) create mode 100644 .translate/state/ifp_advanced.md.yml diff --git a/.translate/state/ifp_advanced.md.yml b/.translate/state/ifp_advanced.md.yml new file mode 100644 index 00000000..ab4ef609 --- /dev/null +++ b/.translate/state/ifp_advanced.md.yml @@ -0,0 +1,6 @@ +source-sha: 60235ebbd73f767403ac39907ac68560fbb444ba +synced-at: "2026-07-18" +model: claude-sonnet-5 +mode: RESYNC +section-count: 7 +tool-version: 0.17.0 diff --git a/lectures/ifp_advanced.md b/lectures/ifp_advanced.md index 0d1cd3b9..866a062f 100644 --- a/lectures/ifp_advanced.md +++ b/lectures/ifp_advanced.md @@ -7,9 +7,39 @@ kernelspec: display_name: Python 3 language: python name: python3 +translation: + title: 收入波动问题 V：资产随机收益 + headings: + Overview: 概述 + The Model: 模型 + The Model::Set Up: 设定 + The Model::Assumptions: 假设 + The Model::Optimality: 最优性 + Solution Algorithm: 求解算法 + Solution Algorithm::A Time Iteration Operator: 时间迭代算子 + Solution Algorithm::Convergence Properties: 收敛性质 + Solution Algorithm::Using an Endogenous Grid: 使用内生网格 + Solution Algorithm::Using an Endogenous Grid::Finding Optimal Consumption: 寻找最优消费 + Solution Algorithm::Using an Endogenous Grid::Iterating: 迭代 + Implementation: 实现 + Simulation: 模拟 + Wealth Inequality: 财富不平等 + Wealth Inequality::Measuring Inequality: 度量不平等 + Exercises: 练习 --- -# 收入波动问题 II：资产随机收益 +```{raw} jupyter +
+ + QuantEcon + +
+``` + +# {index}`收入波动问题 V：资产随机收益 ` + +```{include} _admonition/gpu.md +``` ```{contents} 目录 :depth: 2 @@ -17,7 +47,7 @@ kernelspec: 除了 Anaconda 中的内容外，本讲座还需要以下库： -```{code-cell} ipython +```{code-cell} ipython3 --- tags: [hide-output] --- @@ -26,7 +56,7 @@ tags: [hide-output] ## 概述 -在本讲座中，我们继续研究 {doc}`收入波动问题 `。 +在本讲座中，我们继续研究 {doc}`ifp_egm` 中描述的收入波动问题。 之前假设利率是固定的，但现在我们允许资产收益随状态变化。 @@ -40,20 +70,20 @@ tags: [hide-output] 我们需要以下导入： -```{code-cell} ipython +```{code-cell} ipython3 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 numpy as np -from numba import jit, float64 -from numba.experimental import jitclass -from quantecon import MarkovChain +import quantecon as qe +import jax +import jax.numpy as jnp +from jax import vmap +from typing import NamedTuple +from functools import partial ``` -## 储蓄问题 + + +## 模型 在本节中，我们回顾家庭问题及其最优性结果。 @@ -79,41 +109,37 @@ a_{t+1} = R_{t+1} (a_t - c_t) + Y_{t+1} 初始条件 $(a_0, Z_0)=(a,z)$ 视为给定。 -注意，财富的总收益率序列 ${R_t}_{t \geq 1}$ 允许是随机的。 +与 {doc}`ifp_egm_transient_shocks` 唯一的不同之处在于，财富的总收益率 $\{R_t\}_{t \geq 1}$ 现在允许是随机的。 -序列 $\{Y_t \}_{t \geq 1}$ 是非金融收入。 - -问题的随机成分服从 +具体而言，我们假设 ```{math} :label: eq:RY_func -R_t = R(Z_t, \zeta_t) - \quad \text{且} \quad -Y_t = Y(Z_t, \eta_t), + R_t = R(Z_t, \zeta_t) + \quad \text{且} \quad + Y_t = Y(Z_t, \eta_t), ``` 其中 -* 映射 $R$ 和 $Y$ 是时不变的非负函数， +* $R$ 和 $Y$ 是时不变的非负函数， * 创新过程 $\{\zeta_t\}$ 和 $\{\eta_t\}$ 独立同分布且相互独立， -* $\{Z_t\}_{t \geq 0}$ 是有限集 $\mathsf Z$ 上的不可约齐次马尔可夫链 +* $\{Z_t\}_{t \geq 0}$ 是有限集 $\mathsf Z$ 上的马尔可夫链 令 $P$ 表示链 $\{Z_t\}_{t \geq 0}$ 的马尔可夫矩阵。 -我们对偏好的假设与 {doc}`之前的讲座 ` 中关于收入波动问题的假设相同。 - -如前所述，$\mathbb E_z \hat X$ 表示给定当前值 $Z = z$ 时下一期值 $\hat X$ 的期望。 +在下文中，$\mathbb E_z \hat X$ 表示给定当前值 $Z = z$ 时下一期值 $\hat X$ 的期望。 ### 假设 我们需要一些限制条件来确保目标 {eq}`trans_at` 是有限的，并且下面描述的解法能够收敛。 -我们还需要确保财富的现值值不会增长得太快。 +我们还需要确保财富的现值不会增长得太快。 当 $\{R_t\}$ 是常数时，我们要求 $\beta R < 1$。 -现在它是随机的，我们要求 +现在它是随机的，我们要求（参见 {cite}`ma2020income`） ```{math} :label: fpbc2 @@ -124,25 +150,26 @@ G_R := \lim_{n \to \infty} \left(\mathbb E \prod_{t=1}^n R_t \right)^{1/n} ``` -注意，当 $\{R_t\}$ 取某个常数值 $R$ 时，这一条件简化为之前的限制 $\beta R < 1$ - 值 $G_R$ 可以理解为长期（几何）平均总收益率。 -{cite}`ma2020income`提供了{eq}`fpbc2` 背后的更多直觉。 +为了简化本讲座，我们将*假设利率过程是独立同分布的*。 -我们在下面讨论如何检验该条件。 +在这种情况下，从 $G_R$ 的定义可以清楚地看出 $G_R$ 就是 $\mathbb E R_t$。 + +我们在下面的代码中检验条件 $\beta \mathbb E R_t < 1$。 -最后，我们对非金融收入施加一些常规的技术性限制： +最后，我们对非金融收入施加一些常规的技术性限制。 $$ \mathbb E \, Y_t < \infty \text{ 且 } \mathbb E \, u'(Y_t) < \infty +\label{a:y0} $$ 一个相对简单且满足所有这些限制的环境是 {cite}`benhabib2015` 的独立同分布和 CRRA 环境。 ### 最优性 -令候选消费政策类 $\mathscr C$ 的定义 {doc}`如前 `。 +令候选消费政策类 $\mathscr C$ 的定义如 {doc}`ifp_egm` 中所述。 在 {cite}`ma2020income` 中证明，在所述假设下， @@ -162,15 +189,18 @@ $$ \right\} ``` -（直觉和推导与我们在 {doc}`早期讲座 ` 中关于收入波动问题的内容类似。） +（直觉和推导与 {doc}`ifp_egm` 中的内容类似。） -我们再次使用时间迭代来求解欧拉方程，使用 Coleman--Reffett 算子 $K$ 来匹配欧拉方程 {eq}`ifpa_euler`。 +我们再次使用时间迭代来求解欧拉方程，使用与欧拉方程 {eq}`ifpa_euler` 相匹配的 Coleman--Reffett 算子 $K$ 进行迭代。 ## 求解算法 +```{index} single: Optimal Savings; Computation +``` + ### 时间迭代算子 -我们对候选类 $\sigma \in \mathscr C$ 消费政策的定义与{doc}`之前关于收入波动问题的讲座 `中的定义相同。 +我们对候选类 $\sigma \in \mathscr C$ 消费政策的定义与 {doc}`ifp_egm` 中的定义相同。 对于固定的 $\sigma \in \mathscr C$ 和 $(a,z) \in \mathbf S$，函数 $K\sigma$ 在 $(a,z)$ 处的值 $K\sigma(a,z)$ 定义为满足以下方程的 $\xi \in (0,a]$ @@ -191,7 +221,7 @@ $K$ 背后的思想是，从定义可以看出，$\sigma \in \mathscr C$ 满足 ### 收敛性质 -如前所述，我们在 $\mathscr C$ 上定义如下度量 +如前所述，我们在 $\mathscr C$ 上配以如下度量 $$ \rho(c,d) @@ -212,7 +242,7 @@ $$ ### 使用内生网格 -在研究该模型时，我们发现可以通过 {doc}`内生网格方法 ` 进一步加速时间迭代。 +在研究该模型时，我们发现可以通过 {doc}`内生网格方法 ` 进一步加速时间迭代。 我们将在这里使用相同的方法。 @@ -250,357 +280,540 @@ c_i = #### 迭代 -一旦我们得到 ${s_i, c_i}$ 对，内生资产网格通过 $a_i = c_i + s_i$ 获得。 +一旦我们得到 $\{s_i, c_i\}$ 对，内生资产网格通过 $a_i = c_i + s_i$ 获得。 另外，在上面的讨论中我们固定了 $z \in \mathsf Z$，所以可以将其与 $a_i$ 配对。 -通过在每个 $z$ 上对 ${a_i, c_i}$ 插值，就可以得到政策 $(a,z) \mapsto \sigma(a,z)$ 的近似。 +通过在每个 $z$ 上对 $\{a_i, c_i\}$ 插值，就可以得到政策 $(a,z) \mapsto \sigma(a,z)$ 的近似。 在下面的内容中，我们使用线性插值。 -### 检验假设 - -时间迭代的收敛性依赖于条件 $\beta G_R < 1$ 的满足。 -我们可以利用 $G_R$ 等于矩阵 $L$ 的谱半径这一事实来检验。矩阵 $L$ 定义为 +## 实现 -$$ -L(z, \hat z) := P(z, \hat z) \int R(\hat z, x) \phi(x) dx -$$ +以下是以 `NamedTuple` 表示的模型。 -这个恒等式在 {cite}`ma2020income` 中得到证明，其中 $\phi$ 是资产收益创新 $\zeta_t$ 的密度函数。 +```{code-cell} ipython3 +class IFP(NamedTuple): + """ + 一个 NamedTuple，使用 JAX 存储收入波动问题的基本参数。 + """ + γ: float + β: float + P: jnp.ndarray + a_r: float + b_r: float + a_y: float + b_y: float + s_grid: jnp.ndarray + η_draws: jnp.ndarray + ζ_draws: jnp.ndarray + + +def create_ifp( + γ=1.5, # 效用参数 + β=0.96, # 折现因子 + P=jnp.array([(0.9, 0.1), # Z 的默认马尔可夫链 + (0.1, 0.9)]), + a_r=0.16, # R 冲击中的波动率项 + b_r=0.0, # R 冲击的均值偏移 + a_y=0.2, # Y 冲击中的波动率项 + b_y=0.5, # Y 冲击的均值偏移 + shock_draw_size=100, # 用于蒙特卡洛 + grid_max=100, # 外生网格最大值 + grid_size=100, # 外生网格大小 + seed=1234 # 随机种子 + ): + """ + 使用给定参数创建一个 IFP 实例。 -（注意，$\mathsf Z$ 是一个有限集，所以这个表达式定义了一个矩阵。） + """ + # 假设 {R_t} 独立同分布且 ln R ~ N(b_r, a_r)，检验稳定性 + ER = np.exp(b_r + a_r**2 / 2) + assert β * ER < 1, "稳定性条件不成立。" -当 $\{R_t\}$ 是独立同分布时，检查这一条件甚更容易。 + # 使用 JAX 生成随机抽取 + key = jax.random.PRNGKey(seed) + subkey1, subkey2 = jax.random.split(key) + η_draws = jax.random.normal(subkey1, (shock_draw_size,)) + ζ_draws = jax.random.normal(subkey2, (shock_draw_size,)) + s_grid = jnp.linspace(0, grid_max, grid_size) -在这种情况下，从 $G_R$ 的定义可以清楚地看出 $G_R$ 就是 $\mathbb E R_t$。 + return IFP( + γ, β, P, a_r, b_r, a_y, b_y, s_grid, η_draws, ζ_draws + ) -我们在下面的代码中检验条件 $\beta \mathbb E R_t < 1$。 -## 实现 +def u_prime(c, γ): + """边际效用""" + return c**(-γ) -我们将假设 $R_t = \exp(a_r \zeta_t + b_r)$，其中 $a_r, b_r$ 是常数，$\{\zeta_t\}$ 是独立同分布的标准正态。 +def u_prime_inv(c, γ): + """边际效用的逆函数""" + return c**(-1/γ) -我们允许劳动收入相关，即 +def R(z, ζ, a_r, b_r): + """资产的总收益率""" + return jnp.exp(a_r * ζ + b_r) -$$ -Y_t = \exp(a_y \eta_t + Z_t b_y) -$$ +def Y(z, η, a_y, b_y): + """劳动收入""" + return jnp.exp(a_y * η + (z * b_y)) +``` -其中 $\{\eta_t\}$ 也是独立同分布的标准正态，$\{ Z_t\}$ 是取值于 $\{0, 1\}$ 的马尔可夫链。 - -```{code-cell} ipython -ifp_data = [ - ('γ', float64), # 效用参数 - ('β', float64), # 折现因子 - ('P', float64[:, :]), # z_t 的转移概率 - ('a_r', float64), # R_t 的尺度参数 - ('b_r', float64), # R_t 的加性参数 - ('a_y', float64), # Y_t 的尺度参数 - ('b_y', float64), # Y_t 的加性参数 - ('s_grid', float64[:]), # 储蓄网格 - ('η_draws', float64[:]), # 用于 MC 的创新 η 的抽取 - ('ζ_draws', float64[:]) # 用于 MC 的创新 ζ 的抽取 -] -``` - -```{code-cell} ipython -@jitclass(ifp_data) -class IFP: - """ - 用于收入波动问题的基本类 - """ +这是使用 JAX 的 Coleman-Reffett 算子： - def __init__(self, - γ=1.5, - β=0.96, - P=np.array([(0.9, 0.1), - (0.1, 0.9)]), - a_r=0.1, - b_r=0.0, - a_y=0.2, - b_y=0.5, - shock_draw_size=50, - grid_max=10, - grid_size=100, - seed=1234): - - np.random.seed(seed) # 任意随机种子 - - self.P, self.γ, self.β = P, γ, β - self.a_r, self.b_r, self.a_y, self.b_y = a_r, b_r, a_y, b_y - self.η_draws = np.random.randn(shock_draw_size) - self.ζ_draws = np.random.randn(shock_draw_size) - self.s_grid = np.linspace(0, grid_max, grid_size) - - # 假设 ${R_t}$ 独立同分布并服从下面给定的对数正态分布设定，进行稳定性检验。 - # 检验 β E R_t < 1。 - ER = np.exp(b_r + a_r**2 / 2) - assert β * ER < 1, "稳定性条件不成立。" - - # 边际效用 - def u_prime(self, c): - return c**(-self.γ) - - # 边际效用的逆函数 - def u_prime_inv(self, c): - return c**(-1/self.γ) - - def R(self, z, ζ): - return np.exp(self.a_r * ζ + self.b_r) - - def Y(self, z, η): - return np.exp(self.a_y * η + (z * self.b_y)) -``` - -这是基于 EGM 的 Coleman-Reffett 算子： - -```{code-cell} ipython -@jit -def K(a_in, σ_in, ifp): +```{code-cell} ipython3 +def K( + a_in: jnp.array, # a_in[i, z] 是资产网格 + c_in: jnp.array, # c_in[i, z] = a_in[i, z] 处的消费 + ifp: IFP + ): """ - 收入波动问题的 Coleman--Reffett 算子， - 使用内生网格方法。 + 使用 JAX 结合内生网格方法的收入波动问题的 + Coleman--Reffett 算子。 - * ifp 是 IFP 的实例 - * a_in[i, z] 是资产网格 - * σ_in[i, z] 是在 a_in[i, z] 处的消费 """ - # 简化名称 - u_prime, u_prime_inv = ifp.u_prime, ifp.u_prime_inv - R, Y, P, β = ifp.R, ifp.Y, ifp.P, ifp.β - s_grid, η_draws, ζ_draws = ifp.s_grid, ifp.η_draws, ifp.ζ_draws + # 从 ifp 中提取参数 + γ, β, P, a_r, b_r, a_y, b_y, s_grid, η_draws, ζ_draws = ifp n = len(P) - # 通过线性插值创建消费函数 - σ = lambda a, z: np.interp(a, a_in[:, z], σ_in[:, z]) - - # 分配内存 - σ_out = np.empty_like(σ_in) - - # 在每个 s_i, z 处获得 c_i，存储在 σ_out[i, z] 中， - # 通过蒙特卡洛计算期望项 - for i, s in enumerate(s_grid): - for z in range(n): - # 计算期望 - Ez = 0.0 - for z_hat in range(n): - for η in ifp.η_draws: - for ζ in ifp.ζ_draws: - R_hat = R(z_hat, ζ) - Y_hat = Y(z_hat, η) - U = u_prime(σ(R_hat * s + Y_hat, z_hat)) - Ez += R_hat * U * P[z, z_hat] - Ez = Ez / (len(η_draws) * len(ζ_draws)) - σ_out[i, z] = u_prime_inv(β * Ez) - + def compute_expectation(s, z): + def inner_expectation(z_hat): + def compute_term(η, ζ): + R_hat = R(z_hat, ζ, a_r, b_r) + Y_hat = Y(z_hat, η, a_y, b_y) + a_val = R_hat * s + Y_hat + # 对消费进行插值 + c_interp = jnp.interp(a_val, a_in[:, z_hat], c_in[:, z_hat]) + mu = u_prime(c_interp, γ) + return R_hat * mu + # 对所有冲击组合进行向量化 + η_grid, ζ_grid = jnp.meshgrid(η_draws, ζ_draws, indexing='ij') + terms = vmap(vmap(compute_term))(η_grid, ζ_grid) + return P[z, z_hat] * jnp.mean(terms) + # 对 z_hat 状态求和 + Ez = jnp.sum(vmap(inner_expectation)(jnp.arange(n))) + return u_prime_inv(β * Ez, γ) + + # 对 s_grid 和 z 进行向量化 + compute_exp_v1 = vmap(compute_expectation, in_axes=(None, 0)) + compute_exp_v2 = vmap(compute_exp_v1, in_axes=(0, None)) + c_out = compute_exp_v2(s_grid, jnp.arange(n)) # 计算内生资产网格 - a_out = np.empty_like(σ_out) - for z in range(n): - a_out[:, z] = s_grid + σ_out[:, z] + a_out = s_grid[:, None] + c_out + # 在 (0, 0) 处固定消费-资产对 + c_out = c_out.at[0, :].set(0) + a_out = a_out.at[0, :].set(0) - # 在 (0, 0) 处固定消费-资产对有助于插值 - σ_out[0, :] = 0 - a_out[0, :] = 0 - - return a_out, σ_out + return a_out, c_out ``` -下一个函数通过时间迭代求解最优消费政策的近似。 - -```{code-cell} ipython -def solve_model_time_iter(model, # 包含模型信息的类 - a_vec, # 资产的初始条件 - σ_vec, # 消费的初始条件 - tol=1e-4, - max_iter=1000, - verbose=True, - print_skip=25): - - # 设置循环 - i = 0 - error = tol + 1 - - while i < max_iter and error > tol: - a_new, σ_new = K(a_vec, σ_vec, model) - error = np.max(np.abs(σ_vec - σ_new)) +下一个函数使用 JAX 通过时间迭代求解最优消费政策的近似： + +```{code-cell} ipython3 +@jax.jit +def solve_model( + ifp: IFP, + c_init: jnp.ndarray, # 内生网格上 σ 的初始猜测 + a_init: jnp.ndarray, # 初始内生网格 + tol: float = 1e-5, + max_iter: int = 1000 + ) -> jnp.ndarray: + " 使用 EGM 的时间迭代求解模型。 " + + def condition(loop_state): + c_in, a_in, i, error = loop_state + return (error > tol) & (i < max_iter) + + def body(loop_state): + c_in, a_in, i, error = loop_state + c_out, a_out = K(c_in, a_in, ifp) + error = jnp.max(jnp.abs(c_out - c_in)) i += 1 - if verbose and i % print_skip == 0: - print(f"第{i}迭代的误差是 {error}。") - a_vec, σ_vec = np.copy(a_new), np.copy(σ_new) + return c_out, a_out, i, error - if error > tol: - print("未能收敛！") - elif verbose: - print(f"\n在第{i}次迭代中收敛。") + i, error = 0, tol + 1 + initial_state = (c_init, a_init, i, error) + final_loop_state = jax.lax.while_loop(condition, body, initial_state) + c_out, a_out, i, error = final_loop_state - return a_new, σ_new + return c_out, a_out ``` -现在我们可以用默认参数创建一个实例。 +现在我们可以创建一个实例并使用 JAX 求解模型： -```{code-cell} ipython -ifp = IFP() +```{code-cell} ipython3 +ifp = create_ifp() ``` -接下来我们设置一个初始条件，对应“消费掉所有资产”。 +设置初始条件： -```{code-cell} ipython +```{code-cell} ipython3 # 初始猜测 σ = 消费所有资产 k = len(ifp.s_grid) n = len(ifp.P) -σ_init = np.empty((k, n)) +σ_init = jnp.empty((k, n)) for z in range(n): - σ_init[:, z] = ifp.s_grid -a_init = np.copy(σ_init) + σ_init = σ_init.at[:, z].set(ifp.s_grid) +a_init = σ_init.copy() ``` -让我们生成一个近似解。 +让我们用 JAX 生成一个近似解： -```{code-cell} ipython -a_star, σ_star = solve_model_time_iter(ifp, a_init, σ_init, print_skip=5) +```{code-cell} ipython3 +a_star, σ_star = solve_model(ifp, a_init, σ_init) ``` -这是结果消费政策的图： +让我们再用计时器试一次。 -```{code-cell} ipython -fig, ax = plt.subplots() -for z in range(len(ifp.P)): - ax.plot(a_star[:, z], σ_star[:, z], label=f"当 $z={z}$ 时的消费") +```{code-cell} python3 +with qe.Timer(precision=8): + a_star, σ_star = solve_model(ifp, a_init, σ_init) + a_star.block_until_ready() +``` -plt.legend() -plt.show() +## 模拟 + +让我们回到默认模型，研究资产的平稳分布。 + +我们的计划是让大量家庭向前推进 $T$ 期，然后绘制资产横截面分布的直方图。 + +设置 `num_households=50_000, T=500`。 + +首先我们编写一个函数，将单个家庭向前模拟，并记录资产的最终值。 + +该函数接受一对解 `c_vec` 和 `a_vec`，将其理解为与给定模型 `ifp` 相关联的最优政策。 + +```{code-cell} ipython3 +def simulate_household( + key, a_0, z_idx_0, c_vec, a_vec, ifp, T + ): + """ + 模拟单个家庭 T 期，以逼近资产的平稳分布。 + + - key 是随机数生成器的状态 + - ifp 是 IFP 的一个实例 + - c_vec, a_vec 是 ifp 的最优消费政策和内生网格 + + """ + # 从 ifp 中提取参数 + γ, β, P, a_r, b_r, a_y, b_y, s_grid, η_draws, ζ_draws = ifp + n_z = len(P) + + # 为消费政策创建插值函数 + σ = lambda a, z_idx: jnp.interp(a, a_vec[:, z_idx], c_vec[:, z_idx]) + + # 向前模拟 T 期 + def update(t, state): + a, z_idx = state + # 从 P[z, z'] 中抽取下一期冲击 z' + current_key = jax.random.fold_in(key, 3*t) + z_next_idx = jax.random.choice(current_key, n_z, p=P[z_idx]).astype(jnp.int32) + # 为收入抽取 η 冲击 + η_key = jax.random.fold_in(key, 3*t + 1) + η = jax.random.normal(η_key) + # 为收益率抽取 ζ 冲击 + ζ_key = jax.random.fold_in(key, 3*t + 2) + ζ = jax.random.normal(ζ_key) + # 计算随机收益率 + R_next = R(z_next_idx, ζ, a_r, b_r) + # 计算收入 + Y_next = Y(z_next_idx, η, a_y, b_y) + # 更新资产：a' = R' * (a - c) + Y' + a_next = R_next * (a - σ(a, z_idx)) + Y_next + # 返回更新后的状态 + return a_next, z_next_idx + + initial_state = a_0, z_idx_0 + final_state = jax.lax.fori_loop(0, T, update, initial_state) + a_final, _ = final_state + return a_final ``` -注意，在资产空间的较低区间，我们会消费掉所有资产。 +现在我们编写一个函数，并行模拟许多家庭。 -这是因为我们预期下一期会有收入 $Y_{t+1}$，因此储蓄的紧迫性较低。 +```{code-cell} ipython3 +@partial(jax.jit, static_argnums=(3, 4, 5)) +def compute_asset_stationary( + c_vec, a_vec, ifp, num_households=50_000, T=500, seed=1234 + ): + """ + 模拟 num_households 个家庭 T 期，以逼近资产的平稳分布。 -你能解释为什么在 $z=0$ 时，消费掉所有资产会更早结束（即在较低的资产水平就停止）吗？ + 返回资产持有量的最终横截面。 -### 运动规律 + - ifp 是 IFP 的一个实例 + - c_vec, a_vec 是最优消费政策和内生网格。 -让我们试着了解，在这种消费政策下，从长期来看资产会如何变化。 + """ + # 从 ifp 中提取参数 + γ, β, P, a_r, b_r, a_y, b_y, s_grid, η_draws, ζ_draws = ifp + + # 从 资产 = 储蓄网格最大值 / 2 开始 + a_0_vector = jnp.full(num_households, s_grid[-1] / 2) + # 初始化每个家庭的外生状态 + z_idx_0_vector = jnp.zeros(num_households).astype(jnp.int32) + + # 对许多家庭进行向量化 + key = jax.random.PRNGKey(seed) + keys = jax.random.split(key, num_households) + # 在 (key, a_0, z_idx_0) 上向量化 simulate_household + sim_all_households = jax.vmap( + simulate_household, in_axes=(0, 0, 0, None, None, None, None) + ) + assets = sim_all_households(keys, a_0_vector, z_idx_0_vector, c_vec, a_vec, ifp, T) + + return jnp.array(assets) +``` -与我们在 {doc}`之前关于收入波动问题的讲座` 中一样，我们首先制作一个 45 度图，展示资产的运动规律： +我们需要一些不平等度量来进行可视化，所以让我们先定义它们： -```{code-cell} python3 -# 好状态和坏状态的平均劳动收入 -Y_mean = [np.mean(ifp.Y(z, ifp.η_draws)) for z in (0, 1)] -# 平均收益 -R_mean = np.mean(ifp.R(z, ifp.ζ_draws)) +```{code-cell} ipython3 +def gini_coefficient(x): + """ + 计算数组 x 的基尼系数。 -a = a_star -fig, ax = plt.subplots() -for z, lb in zip((0, 1), ('坏状态', '好状态')): - ax.plot(a[:, z], R_mean * (a[:, z] - σ_star[:, z]) + Y_mean[z] , label=lb) + """ + x = jnp.asarray(x) + n = len(x) + x_sorted = jnp.sort(x) + # 计算基尼系数 + cumsum = jnp.cumsum(x_sorted) + a = (2 * jnp.sum((jnp.arange(1, n+1)) * x_sorted)) / (n * cumsum[-1]) + return a - (n + 1) / n + + +def top_share( + x: jnp.array, # 财富值数组 + p: float=0.01 # 头部家庭的比例（默认 0.01 表示前 1%） + ): + """ + 计算前 p 比例家庭所持有的总财富份额。 -ax.plot(a[:, 0], a[:, 0], 'k--') -ax.set(xlabel='当前资产', ylabel='下一期资产') + """ + x = jnp.asarray(x) + x_sorted = jnp.sort(x) + # 前 p% 中的家庭数量 + n_top = int(jnp.ceil(len(x) * p)) + # 前 p% 持有的财富 + wealth_top = jnp.sum(x_sorted[-n_top:]) + # 总财富 + wealth_total = jnp.sum(x_sorted) + return wealth_top / wealth_total +``` -ax.legend() +现在我们调用该函数，生成资产分布并将其可视化： + +```{code-cell} ipython3 +ifp = create_ifp() +# 提取用于初始化的参数 +s_grid = ifp.s_grid +n_z = len(ifp.P) +a_init = s_grid[:, None] * jnp.ones(n_z) +c_init = a_init +a_vec, c_vec = solve_model(ifp, a_init, c_init) +assets = compute_asset_stationary(c_vec, a_vec, ifp, num_households=200_000) + +# 为图形计算基尼系数 +gini_plot = gini_coefficient(assets) + +# 绘制对数财富直方图 +fig, ax = plt.subplots(figsize=(10, 6)) +ax.hist(jnp.log(assets), bins=40, alpha=0.5, density=True) +ax.set(xlabel='对数资产', ylabel='密度', title="财富分布") +plt.tight_layout() plt.show() ``` -图中实线表示对于每个 $z$，资产的平均更新函数，由下式给出： +直方图显示了对数财富的分布。 -$$ -a \mapsto \bar R (a - \sigma^*(a, z)) + \bar Y(z) -$$ +请记住我们看的是对数值，直方图表明分布有较长的右尾。 -其中 +下面我们更详细地研究这一点。 -* $\bar R = \mathbb E R_t$，即平均收益率，且 -* $\bar Y(z) = \mathbb E_z Y(z, \eta_t)$，即状态 $z$ 下的平均劳动收入。 -虚线是 45 度线。 -从图中可以看出，动态是稳定的——即使在最高的状态下，资产也不会发散。 +## 财富不平等 -## 练习 +让我们通过计算这一现象的一些标准度量来考察财富不平等。 -```{exercise} -:label: ifpa_ex1 +我们还将考察不平等程度如何随利率变化。 + + +### 度量不平等 + +让我们打印出模拟结果中的基尼系数和前 1% 财富份额： + +```{code-cell} ipython3 +gini = gini_coefficient(assets) +top1 = top_share(assets, p=0.01) + +print(f"基尼系数：{gini:.4f}") +print(f"前 1% 财富份额：{top1:.4f}") +``` + +最近的数据表明 + +* 美国财富的基尼系数约为 0.8 +* 前 1% 的财富份额超过 0.3 + +我们具有随机收益的模型生成的基尼系数接近经验值，这表明资本收入风险是财富不平等的一个重要因素。 + +然而，前 1% 的财富份额过大。 + +我们的模型需要适当的校准和进一步的工作——我们暂时搁置这些任务。 -让我们重复 {ref}`之前的练习 `，研究资产的长期横截面分布。 +## 练习 -在那个练习中，我们使用了一个相对简单的收入波动模型。 +```{exercise} +:label: ifp_advanced_ex1 -在解答中，我们发现资产分布的形状不切实际。 +绘制基尼系数如何随资产收益的波动性变化。 -特别是，我们未能匹配财富分布的长右尾。 +具体而言，计算 `a_r` 从 0.10 到 0.16 变化时的基尼系数（至少使用 5 个不同的值），并绘制结果图。 -你的任务是再次尝试这个练习，但这一次使用我们更复杂的模型。 +这告诉我们资本收入风险与财富不平等之间的关系是什么？ -使用默认参数。 ``` -```{solution-start} ifpa_ex1 +```{solution-start} ifp_advanced_ex1 :class: dropdown ``` -首先我们编写一个函数来生成一个较长的资产序列。 +我们对不同的 `a_r` 值进行循环，为每个值求解模型，模拟财富分布，并计算基尼系数。 + +```{code-cell} ipython3 +# 需要探索的 a_r 值范围 +a_r_vals = np.linspace(0.10, 0.16, 5) +gini_vals = [] + +print("正在计算不同收益波动性下的基尼系数...\n") + +for a_r in a_r_vals: + print(f"a_r = {a_r:.3f}...", end=" ") + + # 用这个 a_r 值创建模型 + ifp_temp = create_ifp(a_r=a_r, grid_max=100) + + # 求解模型 + s_grid_temp = ifp_temp.s_grid + n_z_temp = len(ifp_temp.P) + a_init_temp = s_grid_temp[:, None] * jnp.ones(n_z_temp) + c_init_temp = a_init_temp + a_vec_temp, c_vec_temp = solve_model( + ifp_temp, a_init_temp, c_init_temp + ) + + # 模拟家庭 + assets_temp = compute_asset_stationary( + c_vec_temp, a_vec_temp, ifp_temp, num_households=200_000 + ) + + # 计算基尼系数 + gini_temp = gini_coefficient(assets_temp) + gini_vals.append(gini_temp) + print(f"基尼系数 = {gini_temp:.4f}") + +# 绘制结果图 +fig, ax = plt.subplots(figsize=(10, 6)) +ax.plot(a_r_vals, gini_vals, 'o-', linewidth=2, markersize=8) +ax.set(xlabel='收益波动性 (a_r)', + ylabel='基尼系数', + title='财富不平等与收益波动性') +ax.axhline(y=0.8, color='k', linestyle='--', linewidth=1, + label='美国经验基尼系数 (~0.8)') +ax.legend() +plt.tight_layout() +plt.show() +``` -因为我们希望用 JIT 来编译函数，所以在写代码时不得不打破一些良好的编程风格。 +图中显示，财富不平等（用基尼系数衡量）随收益波动性的增大而增加。 -例如，我们会把解 `a_star, σ_star` 以及 `ifp` 一起传入，尽管更自然的方式是只传入 `ifp` 然后在函数内求解。 +这表明资本收入风险是财富不平等的一个关键驱动因素。 -我们这样做的原因是 `solve_model_time_iter` 不是 JIT 编译的。 +当收益波动性更大时，经历了一系列高收益的幸运家庭会积累远多于不幸家庭的财富，从而导致财富分布的不平等程度加剧。 -```{code-cell} python3 -@jit -def compute_asset_series(ifp, a_star, σ_star, z_seq, T=500_000): - """ - 在最优储蓄行为下，模拟长度为 T 的资产时间序列 - * ifp 是 IFP 的实例 - * a_star 是内生网格解 - * σ_star 是网格上的最优消费 - * z_seq 是 {Z_t} 的时间路径 +```{solution-end} +``` - """ +```{exercise} +:label: ifp_advanced_ex2 - # 通过线性插值创建消费函数 - σ = lambda a, z: np.interp(a, a_star[:, z], σ_star[:, z]) +绘制基尼系数如何随劳动收入的波动性变化。 - # 模拟资产路径 - a = np.zeros(T+1) - for t in range(T): - z = z_seq[t] - ζ, η = np.random.randn(), np.random.randn() - R = ifp.R(z, ζ) - Y = ifp.Y(z, η) - a[t+1] = R * (a[t] - σ(a[t], z)) + Y - return a -``` +具体而言，计算 `a_y` 从 0.125 到 0.20 变化时的基尼系数，并绘制结果图。在本练习中设置 `a_r=0.10`。 -接下来，我们调用该函数，生成资产序列，并利用上面的解来绘制直方图。 +这告诉我们劳动收入风险与财富不平等之间的关系是什么？通过改变劳动收入波动性，我们能否达到与改变收益波动性同样程度的不平等上升？ -```{code-cell} python3 -T = 1_000_000 -mc = MarkovChain(ifp.P) -z_seq = mc.simulate(T, random_state=1234) +``` -a = compute_asset_series(ifp, a_star, σ_star, z_seq, T=T) +```{solution-start} ifp_advanced_ex2 +:class: dropdown +``` -fig, ax = plt.subplots() -ax.hist(a, bins=40, alpha=0.5, density=True) -ax.set(xlabel='资产') +我们对不同的 `a_y` 值进行循环，为每个值求解模型，模拟财富分布，并计算基尼系数。 + +```{code-cell} ipython3 +# 需要探索的 a_y 值范围 +a_y_vals = np.linspace(0.125, 0.20, 5) +gini_vals_y = [] + +print("正在计算不同劳动收入波动性下的基尼系数...\n") + +for a_y in a_y_vals: + print(f"a_y = {a_y:.3f}...", end=" ") + + # 用这个 a_y 值和 a_r=0.10 创建模型 + ifp_temp = create_ifp(a_y=a_y, a_r=0.10, grid_max=100) + + # 求解模型 + s_grid_temp = ifp_temp.s_grid + n_z_temp = len(ifp_temp.P) + a_init_temp = s_grid_temp[:, None] * jnp.ones(n_z_temp) + c_init_temp = a_init_temp + a_vec_temp, c_vec_temp = solve_model( + ifp_temp, a_init_temp, c_init_temp + ) + + # 模拟家庭 + assets_temp = compute_asset_stationary( + c_vec_temp, a_vec_temp, ifp_temp, num_households=200_000 + ) + + # 计算基尼系数 + gini_temp = gini_coefficient(assets_temp) + gini_vals_y.append(gini_temp) + print(f"基尼系数 = {gini_temp:.4f}") + +# 绘制结果图 +fig, ax = plt.subplots(figsize=(10, 6)) +ax.plot(a_y_vals, gini_vals_y, 'o-', linewidth=2, markersize=8, color='green') +ax.set(xlabel='劳动收入波动性 (a_y)', + ylabel='基尼系数', + title='财富不平等与劳动收入波动性') +ax.axhline(y=0.8, color='k', linestyle='--', linewidth=1, + label='美国经验基尼系数 (~0.8)') +ax.legend() +plt.tight_layout() plt.show() ``` -我们现在已经成功再现了财富分布的长右尾。 +图中显示，财富不平等随劳动收入波动性的增大而增加，但这一效应比收益波动性的效应要弱得多。 -下面用一张水平小提琴图来展示这一结果的另一种视角。 +比较这两个练习： + +- 当收益波动性（`a_r`）从 0.10 变化到 0.16 时，基尼系数从约 0.20 急剧上升到 0.79 +- 当劳动收入波动性（`a_y`）以类似的百分比幅度从 0.125 变化到 0.20 时，基尼系数虽有增加，但增幅小得多 + +这表明，相较于劳动收入风险，资本收入风险是财富不平等更重要的驱动因素。 + +其直觉在于，财富积累会随时间产生复利效应：经历了有利资产收益的家庭可以将这些收益进行再投资，从而实现指数级增长。 + +相比之下，劳动收入冲击虽然会影响当期消费和储蓄，但对财富积累不具有同样的复利效应。 -```{code-cell} python3 -fig, ax = plt.subplots() -ax.violinplot(a, vert=False, showmedians=True) -ax.set(xlabel='资产') -plt.show() -``` ```{solution-end} -``` +``` \ No newline at end of file From 5c1e1ef8385cb958a6782b6757b941db1e2c1083 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 21:16:16 +1000 Subject: [PATCH 2/4] Restore trailing newline (review nit) Co-Authored-By: Claude Fable 5 --- lectures/ifp_advanced.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/ifp_advanced.md b/lectures/ifp_advanced.md index 866a062f..9d818b17 100644 --- a/lectures/ifp_advanced.md +++ b/lectures/ifp_advanced.md @@ -816,4 +816,4 @@ plt.show() ```{solution-end} -``` \ No newline at end of file +``` From efb9cb47da113e97116d75a2a4e00d5bdf50ae0c Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 21:20:46 +1000 Subject: [PATCH 3/4] 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/ifp_advanced.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lectures/ifp_advanced.md b/lectures/ifp_advanced.md index 9d818b17..c9c073a6 100644 --- a/lectures/ifp_advanced.md +++ b/lectures/ifp_advanced.md @@ -72,6 +72,10 @@ tags: [hide-output] ```{code-cell} ipython3 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 numpy as np import quantecon as qe import jax From ed36bad9077f6f5752902651fa5ed9a7bad7b09c Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Sat, 18 Jul 2026 21:21:53 +1000 Subject: [PATCH 4/4] Point references to untranslated lectures at the English site These {doc} targets exist only in lecture-python.myst until Phase 2 translates them; qualifying with the intermediate: intersphinx prefix gives working links now, and a future resync restores local refs once the targets exist. Program decision recorded 2026-07-18 (Matt). Co-Authored-By: Claude Fable 5 --- lectures/ifp_advanced.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/ifp_advanced.md b/lectures/ifp_advanced.md index c9c073a6..7ca0b874 100644 --- a/lectures/ifp_advanced.md +++ b/lectures/ifp_advanced.md @@ -113,7 +113,7 @@ a_{t+1} = R_{t+1} (a_t - c_t) + Y_{t+1} 初始条件 $(a_0, Z_0)=(a,z)$ 视为给定。 -与 {doc}`ifp_egm_transient_shocks` 唯一的不同之处在于，财富的总收益率 $\{R_t\}_{t \geq 1}$ 现在允许是随机的。 +与 {doc}`intermediate:ifp_egm_transient_shocks` 唯一的不同之处在于，财富的总收益率 $\{R_t\}_{t \geq 1}$ 现在允许是随机的。 具体而言，我们假设 @@ -246,7 +246,7 @@ $$ ### 使用内生网格 -在研究该模型时，我们发现可以通过 {doc}`内生网格方法 ` 进一步加速时间迭代。 +在研究该模型时，我们发现可以通过 {doc}`内生网格方法 ` 进一步加速时间迭代。 我们将在这里使用相同的方法。