打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
数量经济学的Python实现(八)|RCK模型下的计划经济问题

译者:邱任翔,中国人民大学,qrx_math@ruc.edu.cn

「本部分的原文参考为:25. Cass-Koopmans Planning Problem[1]

目录

  • 概览
  • 模型
  • 计划经济问题
  • 打靶法(Shooting Algorithm)
  • 将初始资本存量设置为稳态
  • 大道理论(Turnpike Theory)
  • 无限期问题
  • 总结

概览

本篇讲稿以及下一篇讲稿描述了Tjalling Koopmans和David Cass分析最优经济增长时采用的模型。

这个模型可以看作是索罗模型(Solow model)的拓展,使得储蓄率不再是外生给定,而是最优选择下的结果。

本系列讲稿的这篇和下一篇中,我们介绍这个模型的两个版本。

这两篇文章在一起说明了「计划经济」和去中心化的经济——「竞争性均衡」之间更一般的联系。

这一篇主要介绍计划经济的版本。

本文主要涉及到的观点有:

  • 解决计划经济问题的极小化极大(a min-max problem)问题
  • 解决有初值和终止条件的差分方程的打靶法(shooting algorithm)
  • 大道理论(turnpike property),证明长期经济发展存在着最有效率的资本积累路径
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (115)  #set default figure size
from numba import njit, float64
from numba.experimental import jitclass
import numpy as np

模型

假设时间是离散且有限的,取值为

(在最后我们会探讨的极限情况)

一个商品既可以被消费也可以被投资进实体资本中。

消费品并不是耐用品,在一次消费后就会完全折旧。

被投资的商品是耐用品,每一期都会折旧一部分。

期,我们让表示非耐用的消费品,表示投资额。

= =

代表性的家庭在每一期会被赋予一单位的劳动力,并且他们喜欢消费。

代表性的家庭在每一期会无弹性地提供一单位劳动力,所以对于所有的都有

代表性家庭对于消费束的偏好如下:

其中是折现因子,则反应了单期效用函数的曲率(curvature),越大,则曲率也会越大。

注意:

满足0,u''<0 '="" data-formula-type='inline-equation'>

0' data-formula-type='inline-equation'>意味着消费者消费越多效用越高。

<0 '="" data-formula-type='inline-equation'>意味着边际效用递减。

我们假设是一个外生给定的初始资本存量。

下面我们给出一个全经济范围(economy-wide)的生产函数:

其中,

一个关于的可行性(feasible)约束为:

其中是资本的折旧率。

计划经济问题

一个计划经济者会选择一个配置在可行性约束式(4)的条件下最大化家庭的效用流式(1)。

为非负的拉格朗日乘子序列。

为了找到最优的配置,我们有如下的拉格朗日函数:

则变成如下的极小化极大问题(min-max problem):

  • 极值化(Extremization)意味着在最大化的条件下最大化
  • 我们的问题所满足的条件保证了所需要的二阶条件被满足我们将要计算得到一阶条件的配置所满足。(这段翻译过于复杂,原文是:Our problem satisfies conditions that assure that required second-order conditions are satisfied at an allocation that satisfies the first-order conditions that we are about to compute.)

在计算一阶条件前,我们提供了一些方便的公式。

线性齐次生产函数的有用性质

下面的一些技巧将会帮助我们。

注意到:

定义如下的「人均生产函数」

它的参数是「人均资本拥有量」

这在计算资本的边际产出时是相当有用的:

同理,劳动的边际产出是:

一阶必要条件

我们现在写出拉格朗日函数极值化的一阶必要条件:

(7)式

(8)式

(9)式

(10)式

在计算(9)式时,我们注意到将会同时出现在期和期的可行性约束当中。

(10)式主要来自于「卡鲁什-库恩-塔克条件」(KKT,Karush-Kuhn-Tucker condition)

结合式(7)和式(8),我们有(原文在这里漏了一个beta符号,译者已经修正):

整理一下可以得到:

(12)式

对等式的两边都取的反函数:

最后应用效用函数的具体形式,我们得到了如下的「欧拉方程」(Euler equation):

下面我们定义一个jitclass来存储描述我们经济的参数和方程:

planning_data = [
    ('γ', float64),    # Coefficient of relative risk aversion
    ('β', float64),    # Discount factor
    ('δ', float64),    # Depreciation rate on capital
    ('α', float64),    # Return to capital per capita
    ('A', float64)     # Technology
]
@jitclass(planning_data)
class PlanningProblem():

    def __init__(self, γ=2, β=0.95, δ=0.02, α=0.33, A=1):

        self.γ, self.β = γ, β
        self.δ, self.α, self.A = δ, α, A

    def u(self, c):
        '''
        Utility function
        ASIDE: If you have a utility function that is hard to solve by hand
        you can use automatic or symbolic differentiation
        See https://github.com/HIPS/autograd
        '''

        γ = self.γ

        return c ** (1 - γ) / (1 - γ) if γ!= 1 else np.log(c)

    def u_prime(self, c):
        'Derivative of utility'
        γ = self.γ

        return c ** (-γ)

    def u_prime_inv(self, c):
        'Inverse of derivative of utility'
        γ = self.γ

        return c ** (-1 / γ)

    def f(self, k):
        'Production function'
        α, A = self.α, self.A

        return A * k ** α

    def f_prime(self, k):
        'Derivative of production function'
        α, A = self.α, self.A

        return α * A * k ** (α - 1)

    def f_prime_inv(self, k):
        'Inverse of derivative of production function'
        α, A = self.α, self.A

        return (k / (A * α)) ** (1 / (α - 1))

    def next_k_c(self, k, c):
        ''''
        Given the current capital Kt and an arbitrary feasible
        consumption choice Ct, computes Kt+1 by state transition law
        and optimal Ct+1 by Euler equation.
        '''

        β, δ = self.β, self.δ
        u_prime, u_prime_inv = self.u_prime, self.u_prime_inv
        f, f_prime = self.f, self.f_prime

        k_next = f(k) + (1 - δ) * k - c
        c_next = u_prime_inv(u_prime(c) / (β * (f_prime(k_next) + (1 - δ))))

        return k_next, c_next

下面的代码就建立了一个经济模型:

pp = PlanningProblem()

打靶法(Shooting Algorithm)

我们使用打靶法来计算最优的配置以及其相关的拉格朗日乘子序列

计划经济下的一阶必要条件式(7)(8)(9)形成了一个具有两个边界条件的差分方程系统。

  • 是资本的初始条件
  • 是我们从一阶必要条件的KKT定理中推导得出的资本的终止条件

对于拉格朗日乘子我们并没有初始条件。

如果我们有的话,那么问题就会变得相当简单:-给定,我们就可以分别从式(7)(9)(8)中计算出

  • 继续这种方法我们就可以计算出中的剩余的元素

但是我们并没有的初值,所以这个方法并不奏效。

事实上,我们工作的一部分就是计算出的最优值。

为了计算以及我们需要的其他参数,只需要对上述过程进行一些简单的修改。

这被称为「打靶法」(shooting algorithm)

它是一种猜测与验证反复结合的算法,主要由以下几个步骤组成:

  • 假设一个拉格朗日乘子的初值
  • 应用上述的简单算法
  • 计算并检查其是否等于0
  • 如果,那么问题就解决了。
  • 如果,减小的取值,再试一次。
  • 如果,增加的取值,再试一次。

下面的python代码给出了打靶法对该问题的实现过程。

在下面代码中,我们稍微修改了一点在于我们直接猜测而不是.

@njit
def shooting(pp, c0, k0, T=10):
    '''
    Given the initial condition of capital k0 and an initial guess
    of consumption c0, computes the whole paths of c and k
    using the state transition law and Euler equation for T periods.
    '''

    if c0 > pp.f(k0):
        print('initial consumption is not feasible')

        return None

    # initialize vectors of c and k
    c_vec = np.empty(T+1)
    k_vec = np.empty(T+2)

    c_vec[0] = c0
    k_vec[0] = k0

    for t in range(T):
        k_vec[t+1], c_vec[t+1] = pp.next_k_c(k_vec[t], c_vec[t])

    k_vec[T+1] = pp.f(k_vec[T]) + (1 - pp.δ) * k_vec[T] - c_vec[T]

    return c_vec, k_vec
paths = shooting(pp, 0.20.3, T=10)
fig, axs = plt.subplots(12, figsize=(145))

colors = ['blue''red']
titles = ['Consumption''Capital']
ylabels = ['$c_t$''$k_t$']

T = paths[0].size - 1
for i in range(2):
    axs[i].plot(paths[i], c=colors[i])
    axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])

axs[1].scatter(T+10, s=80)
axs[1].axvline(T+1, color='k', ls='--', lw=1)

plt.show()

显然,我们最初猜测的太高,所以初始的消费偏低。

我们知道这个结论因为在终止的时候并不能使得

现在,我们要让能够进行自动修正,使得最终可以满足的终止条件。

我们使用「二分法」

我们首先猜测的初值。(我们可以不使用,因为就是关于的函数)

我们知道一开始最小可以取到0,最大不超过.

所以我们随便假设一个,最终根据得到的分情况讨论:

如果,这个值就变成的新下界。

如果,则变成新上界。

然后取上下界的中值作为新的进行计算,重复进行直到最后结果收敛。

最后当足够接近0的时候我们就停止计算。

@njit
def bisection(pp, c0, k0, T=10, tol=1e-4, max_iter=500, k_ter=0, verbose=True):

    # initial boundaries for guess c0
    c0_upper = pp.f(k0)
    c0_lower = 0

    i = 0
    while True:
        c_vec, k_vec = shooting(pp, c0, k0, T)
        error = k_vec[-1] - k_ter

        # check if the terminal condition is satisfied
        if np.abs(error) < tol:
            if verbose:
                print('Converged successfully on iteration ', i+1)
            return c_vec, k_vec

        i += 1
        if i == max_iter:
            if verbose:
                print('Convergence failed.')
            return c_vec, k_vec

        # if iteration continues, updates boundaries and guess of c0
        if error > 0:
            c0_lower = c0
        else:
            c0_upper = c0

        c0 = (c0_lower + c0_upper) / 2
def plot_paths(pp, c0, k0, T_arr, k_ter=0, k_ss=None, axs=None):

    if axs is None:
        fix, axs = plt.subplots(13, figsize=(164))
    ylabels = ['$c_t$''$k_t$''$\mu_t$']
    titles = ['Consumption''Capital''Lagrange Multiplier']

    c_paths = []
    k_paths = []
    for T in T_arr:
        c_vec, k_vec = bisection(pp, c0, k0, T, k_ter=k_ter, verbose=False)
        c_paths.append(c_vec)
        k_paths.append(k_vec)

        μ_vec = pp.u_prime(c_vec)
        paths = [c_vec, k_vec, μ_vec]

        for i in range(3):
            axs[i].plot(paths[i])
            axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])

        # Plot steady state value of capital
        if k_ss is not None:
            axs[1].axhline(k_ss, c='k', ls='--', lw=1)

        axs[1].axvline(T+1, c='k', ls='--', lw=1)
        axs[1].scatter(T+1, paths[1][-1], s=80)

    return c_paths, k_paths
plot_paths(pp, 0.30.3, [10]);

将初始资本存量设置为稳态的资本存量

时,最优的资源配置会收敛到稳态

等于是相当具有指导意义的,我们将其称为稳态资本存量

在稳态的情况下对于一个较大的都有.

在稳态的情形下,可行性约束式(4)就变成了:

,带入式(12)得到:

定义,化简得到:

带入生产函数的形式的得到:

我们取, ,得到:

让我们使用python进行验证——把这个稳态的值作为初值带入,

ρ = 1 / pp.β - 1
k_ss = pp.f_prime_inv(ρ+pp.δ)

print(f'steady state for capital is: {k_ss}')
steady state for capital is: 9.57583816331462
plot_paths(pp, 0.3, k_ss, [150], k_ss=k_ss);

显然的,当比较大时,一直在的附近,直到接近为止。

让我们看看当低于时计划者会怎么做:

plot_paths(pp, 0.3, k_ss/3, [150], k_ss=k_ss);

注意到计划者是如何将资本存量增加到稳定状态,在那里停留一段时间并且在接近时逐渐让资本存量变为0。

下图比较了不同的情况下的最优产出:

plot_paths(pp, 0.3, k_ss/3, [150755025], k_ss=k_ss);

大道理论(A Turnpike Property)

下面的计算指出:当非常大时,最优的资本存量将在大多数情况下接近他稳态时的资本存量。

plot_paths(pp, 0.3, k_ss/3, [2501505025], k_ss=k_ss);

注意到当时期增加时,计划者使得接近的时间也更长了。

这个特性反映了稳定状态的大道理论。

一个计划者的经验法则是:

  • 开始,让接近稳态并且在接近前尽可能地保持在这一状态下。

计划者通过调整储蓄率来实现这一点。

让我们计算并画出储蓄率变化图。

@njit
def saving_rate(pp, c_path, k_path):
    'Given paths of c and k, computes the path of saving rate.'
    production = pp.f(k_path[:-1])

    return (production - c_path) / production
def plot_saving_rate(pp, c0, k0, T_arr, k_ter=0, k_ss=None, s_ss=None):

    fix, axs = plt.subplots(22, figsize=(129))

    c_paths, k_paths = plot_paths(pp, c0, k0, T_arr, k_ter=k_ter, k_ss=k_ss, axs=axs.flatten())

    for i, T in enumerate(T_arr):
        s_path = saving_rate(pp, c_paths[i], k_paths[i])
        axs[11].plot(s_path)

    axs[11].set(xlabel='t', ylabel='$s_t$', title='Saving rate')

    if s_ss is not None:
        axs[11].hlines(s_ss, 0, np.max(T_arr), linestyle='--')
plot_saving_rate(pp, 0.3, k_ss/3, [2501507550], k_ss=k_ss)

无限期的经济

我们想设置为

合适的做法是将式(10)的终止条件设置为:

这是一个满足路径最后会收敛到稳态的条件。

我们会估计从任意一个初值开始到达稳态的最优路径。

在下面的代码中我们在一个较大的情况下绘制出相应的消费、资本存量和储蓄率的图像。

我们知道在稳态的情况下,储蓄率是恒定的并且等于

从式(13)可以得到稳定的储蓄率等于:

所以每期的储蓄额是用于抵消每期资本折旧所需。

我们首先研究初值低于稳态的最优资本变化路径:

# steady state of saving rate
s_ss = pp.δ * k_ss / pp.f(k_ss)

plot_saving_rate(pp, 0.3, k_ss/3, [130], k_ter=k_ss, k_ss=k_ss, s_ss=s_ss)

因为, \rho +\delta ' data-formula-type='inline-equation'>

所以计划者采取了一个积极的储蓄策略,其高于稳态情况的储蓄率。

注意到当增加时,<0' data-formula-type='inline-equation'>,所以是一直下降的。

所以计划者缓慢地降低储蓄率知道到达稳态:

练习

  • 绘制初始资本存量是稳态水平1.5倍的最优消费、资本和储蓄率的变化路径,
  • 为什么储蓄率会这样变化?

答案

plot_saving_rate(pp, 0.3, k_ss*1.5, [130], k_ter=k_ss, k_ss=k_ss, s_ss=s_ss)

总结

在下一篇讲稿中,我们要研究一个去中心化版本的经济体,其中用到的技术和偏好结构与本文是完全相同的。

在下一篇中,我们使用亚当·斯密的「看不见的手」代替本文中的计划者。

而代替计划者做出数量决策的是看不见的手所导致的市场价格。

市场的价格必须调整以配合代表性家庭和厂商所独立做出的决策。

本文中计划经济和下一篇中的市场经济之间的关系便是一般均衡理论和福利经济学研究中的基本主题。

Reference

[1]

25. Cass-Koopmans Planning Problem: https://python.quantecon.org/cass_koopmans_1.html

本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
万字长文:预训练词向量模型的方法、应用场景、变体延伸与实践总结
(ok)python3 将print输出的内容保存到txt文件中
42
python怎么实现canopy聚类
零基础入门深度学习(6)
ML之HierarchicalClustering:自定义HierarchicalClustering层次聚类算法
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服