動的計画法で解く効用最大化問題

代表的家計の最適行動

• 離散時間モデル
• 消費者（代表的個人）は無限期間生存するとし，単一の家計を構成¹
• 家計は$t$期において，資産$w_t$を消費$c_t$と貯蓄$s_t$に配分する
• 家計は労働の対価として$p$の所得を得る
  • よって，来期の資産は $w_{t+1} = s_t + p$
• $u(c_t)$を瞬時的な効用関数 (instantenous utility function) とする
• 家計は0期において，次式で表される生涯効用$U_0$を最大化する消費（貯蓄）計画を立てる

U_0 = \sum_{t=0}^{\infty} \beta^t u(c_t)

上記の設定の下で，家計は次式で表現される制約付きの生涯効用最大化問題を解く．

\begin{equation*}
\begin{aligned}
\max_{\{c_t,\,s_t,\,w_{t+1}\}_{t=0}^{\infty}} & \sum_{t=0}^{\infty} \beta^t u(c_t) \\
& \text{s.t.} \quad w_0 \text{ given}, \\
& \quad c_t + s_t = w_t, \\
& \quad w_{t+1} = s_t + p, \\
& \quad c_t, s_t \geq 0.
\end{aligned}
\end{equation*}

ベルマン方程式

$t$期の価値関数を次式で定義する．

\begin{align}
V_t(w_t) = & \max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=0}^{\infty}}
\sum_{\tau=0}^\infty\beta^\tau u(c_{t+\tau}) \\
& \hspace{3em} \begin{aligned}
\text{s.t.} \quad & c_{t+\tau} + s_{t+\tau} = w_{t+\tau}, \\
& w_{t+\tau+1} = s_{t+\tau} + p, \\
& c_{t+\tau}, s_{t+\tau} \geq 0.
\end{aligned}
\end{align}

これを書き下すと，$t$期のベルマン方程式

\begin{align}
V_t(w_t) = \max_{w_{t+1}} &\{u(c_t) + \beta V_{t+1}(w_{t+1})\} \\
& \begin{aligned}
\text{s.t.} \quad & c_t + s_t = w_t, \\
                  & w_{t+1} = s_t + p, \\
                  & c_t, s_t \geq 0.
\end{aligned}
\end{align}

を得る．

途中式が知りたい方向け

導出

\begin{align}
V_t(w_t) = & \max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=0}^{\infty}}
\sum_{\tau=0}^\infty \beta^\tau u(c_{t+\tau}) \quad \text{s.t.} \cdots \\
&= \max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=0}^{\infty}} 
\left\{ u(c_t) + \sum_{\tau=1}^\infty \beta^\tau u(c_{t+\tau}) \right\} \quad \text{s.t.} \cdots \\
&= \max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=0}^{\infty}} 
\left\{ u(c_t) + \beta\sum_{\tau=0}^\infty \beta^\tau u(c_{t+\tau+1}) \right\} \quad \text{s.t.} \cdots \\
&= \max_{\{c_t,\,s_t,\,w_{t+1}\}}
\max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=1}^{\infty}}
\left\{ u(c_t) + \beta\sum_{\tau=0}^\infty \beta^\tau u(c_{t+\tau+1}) \right\} \quad \text{s.t.} \cdots \\
&= \max_{\{c_t,\,s_t,\,w_{t+1}\}}
\left\{ u(c_t) + 
\max_{\{c_{t+\tau},\,s_{t+\tau},\,w_{t+\tau+1}\}_{\tau=1}^{\infty}}
\beta\sum_{\tau=0}^\infty \beta^\tau u(c_{t+\tau+1}) \right\} \quad \text{s.t.} \cdots \\
&= \max_{\{c_t,\,s_t,\,w_{t+1}\}}
\left\{ u(c_t) + 
\beta \max_{\{c_{t+\tau+1},\,s_{t+\tau+1},\,w_{t+\tau+2}\}_{\tau=0}^{\infty}}
\sum_{\tau=0}^\infty \beta^\tau u(c_{t+\tau+1}) \right\} \quad \text{s.t.} \cdots \\
&= \max_{\{c_t,\,s_t,\,w_{t+1}\}}
\left\{ u(c_t) + \beta V_{t+1}(w_{t+1}) \right\} \\
& \hspace{5em}
\begin{aligned}
\text{s.t.} \quad & c_t + s_t = w_t, \\
                  & w_{t+1} = s_t + p, \\
                  & c_t, s_t \geq 0.
\end{aligned}
\end{align}

無限期間の場合はどの期から問題を解いても同じなので，$V_t = V_{t+1} = V$とすると，ベルマン方程式は

\begin{align}
V(w_t) = \max_{w_{t+1}} &\{u(c_t) + \beta V(w_{t+1})\} \\
& \begin{aligned}
\text{s.t.} \quad & c_t + s_t = w_t, \\
                  & w_{t+1} = s_t + p, \\
                  & c_t, s_t \geq 0.
\end{aligned}
\end{align}

と表現することができる．動的計画法の嬉しいところは，このように無限期間の問題を今期と来期の関係性の問題に落とし込むことができる点にある（と筆者は理解している）．

実際にベルマン方程式の解$V$を得るには，価値反復法 (value iteration)という方法を用いれば良い．本記事の問題設定の下では，以下の反復を繰り返す．

V_{j+1} = \max_{\{c_t,\,s_t,\,w_{t+1}\}} \{u(c_t) + \beta V_{j}(w_{t+1})\} \quad \text{s.t.} \cdots

1. $V_0 = (0,\ldots,0)$を価値関数の初期値とする
2. 上の式に従って，$V_1, V_2, \ldots$を繰り返し計算する
3. これを$V_{j+1} \approx V_j$と収束するまで行う

Pythonによる実装

これまで一般形で数式を記載していたが，数値計算ではコンピュータに色々と指示する必要がある：

• 効用関数の具体的な形
• 割引因子$\beta$の値
• 各期における状態数（資産$w_t$が取りうる値の数）
  • 資産$w_t$の最小値と最大値をそれぞれ w_min, w_max として ns で状態数を指定
  • 資産$w_t$を numpy.linspace(w_min, w_max, ns) で離散化
• 価値関数の収束を判定するための評価値 epsilon
• 二つの価値関数$V_{j+1}, V_j$の距離を計算する関数 distance()
  • ここではチェビシェフ距離を採用

使用する関数を定義

効用関数

u(c_t) = c_t^\alpha, \quad \alpha \in (0,1)

効用関数

def u(c, alpha):
    return c ** alpha

試しに可視化してみると，alpha=0.3 の時の効用関数のグラフは次のとおり．

チェビシェフ距離については，NumPyを使ってスマートに実装できるが，ここではベクトルの要素を一つずつ取り出して計算する方式を採用．

チェビシェフ距離

def distance(vector1, vector2):
    # Raise an error if the dimensions of the vectors are different
    if len(vector1) != len(vector2):
        raise ValueError("Vector dimensions are different")
    
    e = 0 # evaluation value
    for i in range(len(vector1)):
        dist = abs(vector1[i] - vector2[i])
        if e == 0 or e < dist:
            e = dist
    return e

数値計算の結果を可視化する関数．plt.plot() でも可視化できるが，状態が離散値なので散布図の方が好ましいと思われる．

可視化用の関数

def plot_result(w, lst, xlabel="", ylabel="", title="", plot_type="scatter"):
    if plot_type == "scatter":
        plt.scatter(w, lst, facecolors="none", edgecolors="tab:blue")
    elif plot_type == "line":
        plt.plot(w, lst, lw=2)
    plt.xlabel(xlabel, fontsize=15)
    plt.ylabel(ylabel, fontsize=15)
    plt.title(title, fontsize=20)
    plt.show()

動的計画法

実装の核となる部分は，無限ループの中で価値反復を行って，価値関数が収束した際にループを終了させるというもの．

dynamic_programming.py

import numpy as np
import matplotlib.pyplot as plt

ns = 100                  # the number of states for each period
w_min, w_max = 10, 40     # the minimum and maximum asset value
alpha = 0.3               # the exponent of c 
beta = 0.9                # discount factor
p = 10                    # payment for labor
accurary = 10 ** 6        # how accurate the convergence criteria is
epsilon = u(w_min, alpha) / accurary # evaluation value for the proximity of two vectors

# discretize the given asset interval [w_min, w_max]
w = np.linspace(w_min, w_max, ns)
# lists that contain corresponding quantity (value, index, consumption, saving) at the period in question
v_lst = [0] * ns
g_lst = [0] * ns
c_lst = [0] * ns
s_lst = [0] * ns

count = 0 # the number of steps for the value function to converge
while True:
    count += 1
    v_tmp = v_lst.copy() # v_tmp is the values in the next period
    for i in range(ns): # i is the index of state at period t
        e_val = 0 # evaluation value
        for j in range(ns): # path (i,j) is considered here
            s = w[j] - p # saving at period t on the path (i,j)
            c = w[i] - s # consumption at period t on the path (i,j)
            if c >= 0:
                vf_t = u(c, alpha) + beta * v_tmp[j] # value function at period t on the path (i,j)
                if vf_t > e_val:
                    e_val = vf_t
                    c_opt = c
                    s_opt = s
                    g = j
        v_lst[i] = e_val # the maximized value at state i
        c_lst[i] = c_opt # the optimal consumption at state i
        s_lst[i] = s_opt # the optimal saving at state i
        g_lst[i] = g  # the optimal direction to go at state i
        delta = distance(v_tmp, v_lst) # the distance btw. the previous and current value functions
    if delta < epsilon:
        print("Converged to the optimal value function.")
        print(f"delta: {delta}, epsilon: {epsilon}")
        break
    print(f"==================== step: {count} ====================") 
    # check if the simulation is running well
    print(f"delta: {delta}, epsilon: {epsilon}")
    print()

print(f"The number of steps required for the value function to converge: {count}")

実行結果

次のように，各ステップにおける価値関数の距離 delta と評価値 epsilon（固定）が出力され，収束すればメッセージ "Converged to the optimal value function." が表示される．

terminal

==================== step: 130 ====================
delta: 2.496224951187287e-06, epsilon: 1.9952623149688796e-06

==================== step: 131 ====================
delta: 2.2466024560685582e-06, epsilon: 1.9952623149688796e-06

==================== step: 132 ====================
delta: 2.0219422047773605e-06, epsilon: 1.9952623149688796e-06

Converged to the optimal value function.
delta: 1.8197479860759813e-06, epsilon: 1.9952623149688796e-06

The number of steps required for the value function to converge: 133

動的計画法で求めた価値関数や消費関数を可視化する．

dynamic_programming.py

# value function
plot_result(w, v_lst, r"$w_t$ (state)", r"$V(w_t)$", "Value function")
# optimal consumption
plot_result(w, c_lst, r"$w_t$ (state)", "Consumption", "Consumption function")
# optimal saving
plot_result(w, s_lst, r"$w_t$ (state)", "Saving", "Saving function")

1. 人間が無限期間生きるというのはおかしな話だが，0期の消費者の生涯効用$U_0$を見ると将来世代の効用関数が入っている．つまり，遺産を残すことが自分の効用の増大にもつながるので，分析者からすると0期の消費者はあたかも無限期間生きているような振る舞いを見せる． ↩