28. 世代交叠模型

在这节课中，我们将学习著名的世代交叠（OLG）模型。政策制定者和研究人员使用这个模型来研究：

• 财政政策

• 货币政策

• 长期增长

以及许多其他主题。

世代交叠模型的第一个严谨版本由保罗·萨缪尔森(Paul Samuelson)提出[Samuelson, 1958]。

我们的目标是深入理解世代交叠模型的一个简单版本。

28.1. 概述

OLG模型的动态与索洛-斯旺增长模型非常相似。

与此同时，OLG模型增加了一个重要的新特征：储蓄多少的选择是内生的。

这一特征为什么重要呢？举例来说，假设我们想要预测一项新税收政策对长期增长的影响。

我们可以在索洛-斯旺模型中引入税收，然后观察稳态的变化。

但这种方法忽略了一个事实：当家庭面对新的税率时，他们会改变自己的储蓄和消费行为。

这样的变化可能会显著改变模型的预测结果。

因此，如果我们关心预测的准确性，就应该对经济主体的决策问题进行建模。

具体而言，模型中的家庭应该根据他们所面临的环境（技术、税收、 价格等）来决定储蓄多少、消费多少。

OLG模型正是为解决这类问题而设计的。

我们将介绍一个简单版本的OLG模型，阐明家庭的决策问题，并研究其对长期增长的影响。

首先，让我们导入需要的库。

import numpy as np
from scipy import optimize
from collections import namedtuple
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']

28.2. 环境

我们假设时间是离散的，因此 \(t=0, 1, 2, \ldots\)。

在时间 \(t\) 出生的个体存活两个时期，\(t\) 和 \(t + 1\)。

我们称一个个体

• 在生命的第一个时期为”年轻人”

• 在生命的第二个时期为”老年人”。

年轻人工作，提供劳动并赚取劳动收入。

他们还决定存多少钱。

老年人不工作，所以所有收入都是金融收入。

他们此前工资的储蓄转化为资本， 这些资本会在与 \(t+1\) 时期新的年轻一代的劳动相结合时获得利息收入。

工资率和利率由市场供求关系在均衡中决定。

为了使代数计算稍微简单一些，我们假设人口规模恒定。

我们将人口规模标准化为1。

我们还假设每个个体提供一个”单位”的劳动时间，因此总劳动供给为1。

28.3. 资本供给

现在让我们从单个家庭的角度开始分析模型。

28.3.1. 消费者问题

假设在时间 \(t\) 出生的个体的效用函数形式为

(28.1)\[ U_t = u(c_t) + \beta u(c_{t+1})\]

其中

• \(u: \mathbb R_+ \to \mathbb R\) 是”瞬时”效用函数（或”流”效用函数）

• \(\beta \in (0, 1)\) 是贴现因子

• \(c_t\) 表示在时间 \(t\) 出生的个体在时间 \(t\) 的消费

• \(c_{t+1}\) 表示同一个体在时间 \(t+1\) 的消费

我们假设效用函数 \(u\) 是严格递增的。

个体的储蓄决策来自于解决以下最优化问题

(28.2)\[ \max_{c_t, c_{t+1}} \, \left \{ u(c_t) + \beta u(c_{t+1}) \right \} \]

受制于

\[ c_t + s_t \le w_t \quad \text{and} \quad c_{t+1} \le R_{t+1} s_t \]

这里

• \(s_t\) 是出生于时间 \(t\) 的个人的储蓄

• \(w_t\) 是时间 \(t\) 的工资率

• \(R_{t+1}\) 是时间 \(t\) 投资的储蓄在时间 \(t+1\) 支付的总利率（即 \(R_{t+1} = 1 + r_{t+1}\)，其中 \(r_{t+1}\) 是净利率）

由于效用函数 \(u\) 是严格递增的，在最大化时这两个约束都会以等式成立。

利用这一事实，将第一个约束中的 \(s_t\) 代入第二个约束，我们得到 \(c_{t+1} = R_{t+1}(w_t - c_t)\)。

最大化的一阶条件可以通过将 \(c_{t+1}\) 代入目标函数， 对 \(c_t\) 求导并令其等于零得到。

这给出了OLG模型的欧拉方程，它刻画了最优的跨期消费动态：

(28.3)\[ u'(c_t) = \beta R_{t+1} u'( R_{t+1} (w_t - c_t))\]

从第一个约束我们得到 \(c_t = w_t - s_t\)，所以欧拉方程 也可以表示为

(28.4)\[ u'(w_t - s_t) = \beta R_{t+1} u'( R_{t+1} s_t)\]

假设对于每个 \(w_t\) 和 \(R_{t+1}\)，恰好有一个 \(s_t\) 可以满足 (28.4)。

那么储蓄可以被写成 \(w_t\) 和 \(R_{t+1}\) 的固定函数。

我们将其表示为

(28.5)\[ s_t = s(w_t, R_{t+1})\]

函数 \(s\) 的具体形式将取决于流效用函数 \(u\) 的选择。

更高的工资会增加可用于储蓄的资源，而利率则影响今天消费与为明天储蓄之间的权衡。

\(w_t\) 和 \(R_{t+1}\) 共同代表经济中的价格（劳动力价格和资本租赁率）。

因此，(28.5) 表示给定价格下的储蓄量。

28.3.2. 特殊情形：对数偏好

在特殊情况 \(u(c) = \log c\) 下，欧拉方程简化为 \(s_t= \beta (w_t - s_t)\)。

求解储蓄，我们得到

(28.6)\[ s_t = s(w_t, R_{t+1}) = \frac{\beta}{1+\beta} w_t\]

在这种特殊情况下，储蓄不依赖于利率。

28.3.3. 储蓄与投资的关系

由于人口规模标准化为1，\(s_t\)也是\(t\)时刻整个经济的总储蓄。

在我们考虑的封闭经济中，不存在外国投资，因此净储蓄等于总投资，这可以理解为对企业的资本供给。

在下一节，我们将研究资本需求。

使供给与需求相等，将使我们能够确定OLG经济中的均衡状态。

28.4. 资本需求

下面我们先描述企业的问题，然后写出一个方程，用于描述给定价格条件下的资本需求。

28.4.1. 企业问题

对于每个整数 \(t \geq 0\)，t期的产出 \(y_t\) 由**柯布-道格拉斯生产函数**给出

(28.7)\[ y_t = k_t^{\alpha} \ell_t^{1-\alpha}\]

在这里，\(k_t\) 表示资本，\(\ell_t\) 表示劳动，而 \(\alpha\) 是一个参数 （有时被称为”资本的产出弹性”）。

企业面临的利润最大化问题是

(28.8)\[ \max_{k_t, \ell_t} \{ k^{\alpha}_t \ell_t^{1-\alpha} - R_t k_t -w_t \ell_t \}\]

一阶条件是通过对资本和劳动分别求导，并令导数等于零得到的：

\[ (1-\alpha)(k_t / \ell_t)^{\alpha} = w_t \quad \text{and} \quad \alpha (k_t / \ell_t)^{\alpha - 1} = R_t\]

这些条件将工资 \(w_t\) 等同于劳动的边际产出，将利率 \(R_t\) 等同于资本的边际产出。

28.4.2. 需求

根据我们的假设\(\ell_t = 1\)，我们可以得到

(28.9)\[ w_t = (1-\alpha)k_t^\alpha \]

并且

(28.10)\[ R_t = \alpha k_t^{\alpha - 1} \]

重新整理 (28.10) 得出在时间 \(t+1\) 的总资本需求

(28.11)\[ k^d (R_{t+1}) := \left (\frac{\alpha}{R_{t+1}} \right )^{1/(1-\alpha)}\]

这里 \(k^d(R_{t+1})\) 是当总利率为 \(R_{t+1}\) 时，企业希望使用的资本量。

由于资本的边际产出是递减的，资本需求是其成本 \(R_{t+1}\) 的递减函数。

在Python中这可以表示为

def capital_demand(R, α):
    return (α/R)**(1/(1-α)) 

def capital_supply(R, β, w):
    R = np.ones_like(R)
    return R * (β / (1 + β)) * w

28.5. 均衡

在本节中，我们将推导均衡条件并研究一个例子。

28.5.1. 均衡条件

在均衡状态下，\(t\)时刻的储蓄等于\(t\)时刻的投资，也等于\(t+1\)时刻的资本供给。

通过将这些数量设为相等来计算均衡，即

(28.12)\[ s(w_t, R_{t+1}) = k^d(R_{t+1}) = \left (\frac{\alpha}{R_{t+1}} \right )^{1/(1-\alpha)}\]

原则上，我们现在可以给定\(w_t\)来求解均衡价格\(R_{t+1}\)。 （实际上，我们首先需要指定函数\(u\)，从而得到\(s\)。）

当我们求解这个关于\(t+1\)时刻结果的方程时，\(t\)时刻的量已经确定，所以我们可以将\(w_t\)视为常数。

从均衡\(R_{t+1}\)和(28.11)，我们可以得到 均衡数量\(k_{t+1}\)。

28.5.2. 示例：对数效用

在对数效用的情况下，我们可以使用(28.12)和(28.6)得到

(28.13)\[ \frac{\beta}{1+\beta} w_t = \left( \frac{\alpha}{R_{t+1}} \right)^{1/(1-\alpha)}\]

求解均衡利率得到

(28.14)\[ R_{t+1} = \alpha \left( \frac{\beta}{1+\beta} w_t \right)^{\alpha-1}\]

在Python中，我们可以通过以下方式计算这个

def equilibrium_R_log_utility(α, β, w):
    R = α * ( (β * w) / (1 + β))**(α - 1)
    return R

在对数效用的情况下，由于资本供给不依赖于利率，均衡数量是由供给固定的。

也就是说，

(28.15)\[ k_{t+1} = s(w_t, R_{t+1}) = \frac{\beta }{1+\beta} w_t\]

下图绘制了资本供给和需求作为利率 \(R_{t+1}\) 的函数，同时标出了均衡数量和均衡价格。

R_vals = np.linspace(0.3, 1)
α, β = 0.5, 0.9
w = 2.0

fig, ax = plt.subplots()

ax.plot(R_vals, capital_demand(R_vals, α), 
        label="总需求")
ax.plot(R_vals, capital_supply(R_vals, β, w), 
        label="总供给")

R_e = equilibrium_R_log_utility(α, β, w)
k_e = (β / (1 + β)) * w

ax.plot(R_e, k_e, 'o',label='均衡')

ax.set_xlabel("$R_{t+1}$")
ax.set_ylabel("$k_{t+1}$")
ax.legend()
plt.show()

findfont: Failed to find font weight normal, now using 600.
findfont: Failed to find font weight normal, now using 600.

28.6. 动态分析

在本节中，我们将讨论动态问题。

目前我们将专注于对数效用的情况，此时均衡由(28.15)确定。

28.6.1. 资本积累过程

前面我们已经展示了如何根据\(w_t\)得到均衡的\(k_{t+1}\)。

根据(28.9)，我们可以将其转化为\(k_{t+1}\)作为\(k_t\)的函数。

具体来说，由于\(w_t = (1-\alpha)k_t^\alpha\)，我们有

(28.16)\[ k_{t+1} = \frac{\beta}{1+\beta} (1-\alpha)(k_t)^{\alpha}\]

如果我们对这个方程进行迭代，我们将得到一个资本存量序列。

让我们绘制这些动态的45度图，我们将其表示为

\[ k_{t+1} = g(k_t) \quad \text{其中 } g(k) := \frac{\beta}{1+\beta} (1-\alpha)(k)^{\alpha} \]

def k_update(k, α, β):
    return β * (1 - α) * k**α /  (1 + β)

α, β = 0.5, 0.9
kmin, kmax = 0, 0.1
n = 1000
k_grid = np.linspace(kmin, kmax, n)
k_grid_next = k_update(k_grid,α,β)

fig, ax = plt.subplots(figsize=(6, 6))

ax.plot(k_grid, k_grid_next,  lw=2, alpha=0.6, label='$g$')
ax.plot(k_grid, k_grid, 'k-', lw=1, alpha=0.7, label=r'$45^{\circ}$')


ax.legend(loc='upper left', frameon=False, fontsize=12)
ax.set_xlabel('$k_t$', fontsize=12)
ax.set_ylabel('$k_{t+1}$', fontsize=12)

plt.show()

findfont: Failed to find font weight normal, now using 600.
findfont: Failed to find font weight normal, now using 600.

28.6.2. 稳态（对数情况）

图表显示，该模型具有唯一的正稳态，我们将其表示为\(k^*\)。

我们可以通过设置\(k^* = g(k^*)\)来求解\(k^*\)，即

(28.17)\[ k^* = \frac{\beta (1-\alpha) (k^*)^{\alpha}}{(1+\beta)}\]

求解方程得到

(28.18)\[ k^* = \left (\frac{\beta (1-\alpha)}{1+\beta} \right )^{1/(1-\alpha)}\]

我们可以从(28.10)得到稳态利率，即

\[ R^* = \alpha (k^*)^{\alpha - 1} = \frac{\alpha}{1 - \alpha} \frac{1 + \beta}{\beta} \]

在Python中，我们有

k_star = ((β * (1 - α))/(1 + β))**(1/(1-α))
R_star = (α/(1 - α)) * ((1 + β) / β)

28.6.3. 时间序列

上面的45度图显示，具有正初始条件的资本时间序列会收敛到这个稳态。

让我们绘制一些时间序列来可视化这一点。

ts_length = 25
k_series = np.empty(ts_length)
k_series[0] = 0.02
for t in range(ts_length - 1):
    k_series[t+1] = k_update(k_series[t], α, β)

fig, ax = plt.subplots()
ax.plot(k_series, label="资本序列")
ax.plot(range(ts_length), np.full(ts_length, k_star), 'k--', label=fr'$k^* = {k_star:.4f}$')
ax.set_ylim(0, 0.1)
ax.set_ylabel("资本")
ax.set_xlabel("$t$")
ax.legend()
plt.show()

如果你尝试不同的正初始值，你会发现这个序列总是收敛于\(k^*\)。

下面我们还绘制了随时间变化的总利率。

R_series = α * k_series**(α - 1)

fig, ax = plt.subplots()
ax.plot(R_series, label="毛利率")
ax.plot(range(ts_length), np.full(ts_length, R_star), 'k--', label="$R^*$")
ax.set_ylim(0, 4)
ax.set_ylabel("毛利率")
ax.set_xlabel("$t$")
ax.legend()
plt.show()

利率反映了资本的边际产出，当资本存量较低时，利率较高。

28.7. CRRA偏好

在本节中，我们将对数效用推广为CRRA（常相对风险厌恶）效用函数

\[ u(c) = \frac{ c^{1-\gamma}-1}{1-\gamma} \]

其中\(\gamma >0, \gamma\neq 1\)。

CRRA效用在宏观经济学中很受欢迎，因为参数\(\gamma\)控制着相对风险厌恶系数。

它还控制着跨期替代弹性，后者决定了个体在利率变化时愿意在不同时期之间调整消费的程度。

上面研究的对数效用情形是一个特例，当\(\gamma \to 1\)时得到。

在其他方面，模型是相同的。

下面我们在Python中定义效用函数，并构造一个namedtuple来存储参数。

def crra(c, γ):
    return c**(1 - γ) / (1 - γ)

Model = namedtuple('Model', ['α',        # 柯布-道格拉斯参数
                             'β',        # 折现因子
                             'γ']        # CRRA效用函数中的参数
                   )

def create_olg_model(α=0.4, β=0.9, γ=0.5):
    return Model(α=α, β=β, γ=γ)

让我们也重新定义资本需求函数，使其能够与这个namedtuple一起工作。

def capital_demand(R, model):
    return (model.α/R)**(1/(1-model.α)) 

28.7.1. 供给

对于家庭来说，欧拉方程变为

(28.19)\[ (w_t - s_t)^{-\gamma} = \beta R^{1-\gamma}_{t+1} (s_t)^{-\gamma}\]

求解储蓄，我们得到

(28.20)\[ s_t = s(w_t, R_{t+1}) = w_t \left [ 1 + \beta^{-1/\gamma} R_{t+1}^{(\gamma-1)/\gamma} \right ]^{-1}\]

请注意，与对数情况不同，储蓄现在取决于利率。

def savings_crra(w, R, model):
    α, β, γ = model
    return w / (1 + β**(-1/γ) * R**((γ-1)/γ)) 

model = create_olg_model()
w = 2.0

fig, ax = plt.subplots()

ax.plot(R_vals, capital_demand(R_vals, model), 
        label="总需求")
ax.plot(R_vals, savings_crra(w, R_vals, model), 
        label="总供给")

ax.set_xlabel("$R_{t+1}$")
ax.set_ylabel("$k_{t+1}$")
ax.legend()
plt.show()

28.7.2. 均衡

将资本的总需求（参见(28.11)）与我们新的总供给函数相等，即可得出均衡资本。

因此，我们设定

(28.21)\[ w_t \left [ 1 + \beta^{-1/\gamma} R_{t+1}^{(\gamma-1)/\gamma} \right ]^{-1} = \left (\frac{R_{t+1}}{\alpha} \right )^{1/(\alpha - 1)}\]

这个表达式相当复杂，我们无法解析地求解 \(R_{t+1}\)。

将(28.10)和(28.21)结合起来，得到

(28.22)\[ k_{t+1} = \left [ 1 + \beta^{-1/\gamma} (\alpha k^{\alpha - 1}_{t+1})^{(\gamma-1)/\gamma} \right ]^{-1} (1-\alpha)(k_t)^{\alpha}\]

同样，有了这个等式和给定的 \(k_t\)，我们无法用纸笔求解 \(k_{t+1}\)。

在下面的练习中，你将被要求用数值方法解这些方程。

28.8. 练习

Exercise 28.1

使用(28.22)在 CRRA 情况下数值求解均衡资本存量的动态。

使用 45 度图来可视化这个动态过程。

Solution to Exercise 28.1

为了在给定 \(k_t\) 的情况下求解 \(k_{t+1}\)，我们使用牛顿法。

设

(28.23)\[ f(k_{t+1}, k_t) = k_{t+1} \left[ 1 + \beta^{-1/\gamma} \left ( \alpha k^{\alpha-1}_{t+1} \right )^{(\gamma-1)/\gamma} \right] - (1-\alpha) k^{\alpha}_t =0\]

如果给定了 \(k_t\)，那么 \(f\) 就是未知数 \(k_{t+1}\) 的函数。

然后我们可以使用 scipy.optimize.newton 求解 \(f(k_{t+1}, k_t)=0\) 来得到 \(k_{t+1}\)。

首先让我们定义 \(f\)。

def f(k_prime, k, model):
    α, β, γ = model.α, model.β, model.γ
    z = (1 - α) * k**α
    a = α**(1-1/γ)
    b = k_prime**((α * γ - α + 1) / γ)
    p = k_prime + k_prime * β**(-1/γ) * a * b
    return p - z

现在让我们定义一个函数来找出 \(k_{t+1}\) 的值。

def k_update(k, model):
    return optimize.newton(lambda k_prime: f(k_prime, k, model), 0.1)

最后，这是 45 度图。

kmin, kmax = 0, 0.5
n = 1000
k_grid = np.linspace(kmin, kmax, n)
k_grid_next = np.empty_like(k_grid)

for i in range(n):
    k_grid_next[i] = k_update(k_grid[i], model)

fig, ax = plt.subplots(figsize=(6, 6))

ax.plot(k_grid, k_grid_next,  lw=2, alpha=0.6, label='$g$')
ax.plot(k_grid, k_grid, 'k-', lw=1, alpha=0.7, label=r'$45^{\circ}$')


ax.legend(loc='upper left', frameon=False, fontsize=12)
ax.set_xlabel('$k_t$', fontsize=12)
ax.set_ylabel('$k_{t+1}$', fontsize=12)

plt.show()

Exercise 28.2

上一个练习中的45度图表明存在唯一的正稳态。

通过在(28.22)中设置 \(k_{t+1} = k_t = k^*\)，可以得到正稳态，即：

\[ k^* = \frac{(1-\alpha)(k^*)^{\alpha}} {1 + \beta^{-1/\gamma} (\alpha (k^*)^{\alpha-1})^{(\gamma-1)/\gamma}} \]

与对数效用情况不同，CRRA效用的稳态 \(k^*\) 无法通过解析方法获得。

相反，我们使用牛顿法求解 \(k^*\)。

Solution to Exercise 28.2

我们引入一个函数 \(h\)，使得 正稳态是 \(h\) 的根。

(28.24)\[ h(k^*) = k^* \left [ 1 + \beta^{-1/\gamma} (\alpha (k^*)^{\alpha-1})^{(\gamma-1)/\gamma} \right ] - (1-\alpha)(k^*)^{\alpha}\]

下面是Python代码实现

def h(k_star, model):
    α, β, γ = model.α, model.β, model.γ
    z = (1 - α) * k_star**α
    R1 = α ** (1-1/γ)
    R2 = k_star**((α * γ - α + 1) / γ)
    p = k_star + k_star * β**(-1/γ) * R1 * R2
    return p - z

让我们运用牛顿法来求根：

k_star = optimize.newton(h, 0.2, args=(model,))
print(f"k_star = {k_star}")

k_star = 0.25788950250843484

Exercise 28.3

根据上面列出的参数化，生成三个资本的时间路径，来自三个不同的初始条件。

使用 \(k_0\) 的初始条件为 \(0.001, 1.2, 2.6\)，并且时间序列长度为 10。

Solution to Exercise 28.3

让我们定义常数和三个不同的初始条件

ts_length = 10
k0 = np.array([0.001, 1.2, 2.6])

def simulate_ts(model, k0_values, ts_length):

    fig, ax = plt.subplots()

    ts = np.zeros(ts_length)

    # 模拟并绘制时间序列
    for k_init in k0_values:
        ts[0] = k_init
        for t in range(1, ts_length):
            ts[t] = k_update(ts[t-1], model)
        ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6,
                label=r'$k_0=%g$' %k_init)
    ax.plot(np.arange(ts_length), np.full(ts_length, k_star),
            alpha=0.6, color='red', label=r'$k^*$')
    ax.legend(fontsize=10)

    ax.set_xlabel(r'$t$', fontsize=14)
    ax.set_ylabel(r'$k_t$', fontsize=14)

    plt.show()

simulate_ts(model, k0, ts_length)

findfont: Failed to find font weight normal, now using 600.
findfont: Failed to find font weight normal, now using 600.