Appendix B — 連続時間モデル

using LinearAlgebra
using SparseArrays
using Plots
using LaTeXStrings
using QuantEcon
using PrettyTables

本文の章はすべて離散時間でモデルを組んできましたが, 現代のマクロ経済学では連続時間の定式化も広く使われています. この章の前半では, 連続時間の確率過程 (ブラウン運動, 伊藤積分, OU 過程) と最大値原理 (HJB 方程式) という理論の道具を整えます. 後半では, その道具を使って Aiyagari モデルの連続時間版を有限差分法 (finite difference method) で解きます (Achdou et al. 2022).

B.1 確率過程

連続時間の確率過程は, 時間 $t$ に対してランダムな値をとる関数 $x_t$ の族 $(x_t)_{t \geq 0}$ として定義されます. 連続時間の確率過程として基礎になるのは, 連続的なノイズを表すブラウン運動 (Brownian motion) です.

B.1.1 ブラウン運動

定義 B.1 (ブラウン運動) $b_0 = 0$ を満たす確率過程 $(b_t)_{t \geq 0}$ は次の性質を満たすとき, ブラウン運動と呼ばれる.

独立増分: 任意の $0 \leq t_1 < t_2 < \cdots < t_n$ に対して, $b_{t_2} - b_{t_1}, b_{t_3} - b_{t_2}, \ldots, b_{t_n} - b_{t_{n-1}}$ は独立である.
正規性: 任意の $t$ に対して, $b_t \sim \mathcal{N}(\mu t, \sigma^2 t)$
連続パス: $b_t$ は $t$ の関数として連続である.

特に, $\mu = 0, \sigma = 1$ の場合を標準ブラウン運動と呼ぶ.

なお標準ブラウン運動 $B_t$ は明らかに $B_{t+1} - B_t \sim \mathcal{N}(0, 1)$ を満たすので, 離散時間のホワイトノイズの連続時間版とみなすことができます.

後で伊藤積分や伊藤の補題で繰り返し使う $(dB_t)^2 = dt$ という計算規則は, ブラウン運動の二次変分 (quadratic variation) が時間 $t$ に等しいという性質に由来します.

命題 B.1 (ブラウン運動の二次変分) 標準ブラウン運動 $B_t$ に対して, 区間 $[0, t]$ の分割 $0 = s_0 < s_1 < \cdots < s_n = t$ を考え, その幅 (メッシュ) を $\|P\| := \max_i (s_{i+1} - s_i)$ とする. メッシュを $0$ に近づけると, 増分の二乗和は $t$ に二乗平均収束する. \[ \sum_{i=0}^{n-1} \left(B_{s_{i+1}} - B_{s_i}\right)^2 \;\xrightarrow{\ L^2\ }\; t \qquad (\|P\| \to 0). \] これを微分の記法で略記したものが $(dB_t)^2 = dt$ である.

証明. 二乗和を $S := \sum_{i=0}^{n-1} \Delta_i^2$ (ただし $\Delta_i := B_{s_{i+1}} - B_{s_i}$, 区間幅 $\delta_i := s_{i+1} - s_i$) と書く. $S$ が $t$ に二乗平均収束 ($L^2$ 収束) するとは, 定義により \[ \mathbb{E}\left[(S - t)^2\right] \longrightarrow 0 \qquad (\|P\| \to 0) \] が成り立つことである. これを示す. ブラウン運動の正規性と独立増分性から $\Delta_i \sim \mathcal{N}(0, \delta_i)$ で互いに独立である.

Step1: $\mathbb{E}[\Delta_i^2] = \text{Var}(\Delta_i) = \delta_i$ なので, 区間幅の総和が $t$ であることから \[ \mathbb{E}[S] = \sum_i \delta_i = t. \]

Step2: $\Delta_i^2$ は独立なので分散は和になる. 正規分布 $\mathcal{N}(0, \delta_i)$ は $\mathbb{E}[\Delta_i^4] = 3\delta_i^2$ を満たすから, $\text{Var}(\Delta_i^2) = 3\delta_i^2 - \delta_i^2 = 2\delta_i^2$ であり, \[ \text{Var}(S) = \sum_i 2\delta_i^2 \le 2 \left(\max_i \delta_i\right) \sum_i \delta_i = 2 t\, \|P\|. \]

Step3: $\mathbb{E}[S] = t$ なので二乗平均誤差は分散に一致し, \[ \mathbb{E}\left[(S - t)^2\right] = \text{Var}(S) \le 2 t\, \|P\| \xrightarrow{\ \|P\| \to 0\ } 0. \] よって $S \xrightarrow{\ L^2\ } t$ が成り立つ.

B.1.2 伊藤積分

ここからブラウン運動を含む確率微分方程式を扱います. このとき, 伊藤積分を用います.

定義 B.2 (伊藤積分) 過去の情報のみに依存する関数 $f_s$ と標準ブラウン運動 $B_t$ に対して, 伊藤積分 $\int_0^t f_s\, dB_s$ を, 区間分割 $0 = s_0 < s_1 < \cdots < s_n = t$ に沿った和の二乗平均極限として定義する. \[ \int_0^t f_s\, dB_s := \lim_{\max_i |s_{i+1} - s_i| \to 0} \sum_{i=0}^{n-1} f_{s_i}\,\left(B_{s_{i+1}} - B_{s_i}\right). \] このとき次の2つの性質が成り立つ.

平均ゼロ: $\displaystyle \mathbb{E}\left[\int_0^t f_s\, dB_s\right] = 0$
伊藤の等長性 (Itô isometry): $\displaystyle \mathbb{E}\left[\left(\int_0^t f_s\, dB_s\right)^2\right] = \int_0^t \mathbb{E}\left[f_s^2\right]\, ds$

各増分 $B_{s_{i+1}} - B_{s_i}$ を区間の左端 $s_i$ の値 $f_{s_i}$ で評価する点が重要です. これにより, $f_{s_i}$ は区間 $[s_i, s_{i+1}]$ の将来の増分 $B_{s_{i+1}} - B_{s_i}$ に依存しない (独立である) という性質が使えます. この独立性を踏まえると, 平均ゼロは直感的です. \[ \mathbb{E}\left[f_{s_i} \Delta_i\right] = \mathbb{E}[f_{s_i}]\, \mathbb{E}[\Delta_i] = 0. \]

等長性に関しては, 和の二乗を展開すると, 異なる区間の積 ($i \neq j$) は独立かつ平均ゼロより消え, 対角項だけが残ります. 各対角項は $\mathbb{E}[f_{s_i}^2\, \Delta_i^2] = \mathbb{E}[f_{s_i}^2]\, \mathbb{E}[\Delta_i^2]$ となるので, 命題 B.1 の二次変分を使うと, \[ \sum_i \mathbb{E}[f_{s_i}^2\, \Delta_i^2] = \sum_i \mathbb{E}[f_{s_i}^2]\, (s_{i+1} - s_i) \xrightarrow{\ n \to \infty\ } \int_0^t \mathbb{E}[f_s^2]\, ds. \] 微分の記法でいえば, $\mathbb{E}[dB_t] = 0$ が平均ゼロを, $(dB_t)^2 = dt$ が等長性を生んでいる, と見ることができます.

命題 B.2 (伊藤積分の正規性) 被積分関数 $f_s$ が決定論的な ($B_s$ に依存しない) とき, 伊藤積分は正規分布に従う. \[ \int_0^t f_s\, dB_s \;\sim\; \mathcal{N}\left(0,\ \int_0^t f_s^2\, ds\right). \]

証明. 伊藤積分は近似和 $\sum_i f_{s_i}\, \Delta_i$ ($\Delta_i := B_{s_{i+1}} - B_{s_i}$) の $L^2$ 極限である. $f_s$ が決定論的なら $f_{s_i}$ は定数なので, 各項は $f_{s_i}\, \Delta_i \sim \mathcal{N}\left(0,\, f_{s_i}^2 (s_{i+1} - s_i)\right)$ で正規である. さらにブラウン運動の独立増分性から $\Delta_i$ は互いに独立なので, 独立な正規確率変数の和である近似和もまた正規分布に従う. \[ \sum_i f_{s_i}\, \Delta_i \;\sim\; \mathcal{N}\left(0,\ \sum_i f_{s_i}^2 (s_{i+1} - s_i)\right). \] この近似和は $\|P\| \to 0$ で伊藤積分に $L^2$ 収束し, その分散は $\sum_i f_{s_i}^2 (s_{i+1} - s_i) \to \int_0^t f_s^2\, ds$ に収束する. 正規分布の族は分布収束に関して閉じている¹ ので, 平均ゼロ性と合わせて, 極限である伊藤積分も正規分布 $\mathcal{N}\left(0,\, \int_0^t f_s^2\, ds\right)$ に従う.

B.1.3 Orstein-Uhlenbeck 過程

AR(1) 過程の連続時間版として, Orstein-Uhlenbeck (OU) 過程を考えます. OU過程は以下の確率微分方程式で定義されます.

定義 B.3 (OU過程) OU過程 $x_t$ は以下の確率微分方程式で定義される. \[ dx_t = -\gamma (x_t - \mu) dt + \sigma dB_t. \] ここで $B_t$ は標準ブラウン運動.

命題 B.3 (OU過程の解) OU過程 $x_t$ は次のように陽に解ける. \[ x_t = \mu + (x_0 - \mu) e^{-\gamma t} + \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s. \tag{B.1}\]

証明. 確率微分方程式を積分因子 (integrating factor) $e^{\gamma t}$ を用いて解く. 積分因子は決定論的なので, 通常の積の微分法則が使えて \[ d\left(e^{\gamma t}(x_t - \mu)\right) = \gamma e^{\gamma t}(x_t - \mu)\, dt + e^{\gamma t}\, dx_t. \] ここに $dx_t = -\gamma (x_t - \mu)\, dt + \sigma\, dB_t$ を代入するとドリフト項が打ち消し合い, \[ d\left(e^{\gamma t}(x_t - \mu)\right) = e^{\gamma t}\, \sigma\, dB_t. \] 両辺を $0$ から $t$ まで積分する. 左辺は微分を積分した増分 $\left[e^{\gamma s}(x_s - \mu)\right]_0^t$ なので, \[ e^{\gamma t}(x_t - \mu) - (x_0 - \mu) = \sigma \int_0^t e^{\gamma s}\, dB_s. \] 両辺を $e^{\gamma t}$ で割ると, 解が陽に求まる. \[ x_t = \mu + (x_0 - \mu)\, e^{-\gamma t} + \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s. \]

命題 B.4 (OU過程の平均と分散) \[ \mathbb{E}[x_t] = \mu + (x_0 - \mu)e^{-\gamma t}, \quad \text{Var}(x_t) = \frac{\sigma^2}{2\gamma} \left(1 - e^{-2\gamma t}\right). \] また, $t \to \infty$ のとき, $x_t$ は平均 $\mu$, 分散 $\frac{\sigma^2}{2\gamma}$ の正規分布に収束する.

証明. Step1: 式 B.1 の期待値をとる. 右辺第3項は伊藤積分なので平均ゼロであり, \[ \mathbb{E}[x_t] = \mu + (x_0 - \mu)\, e^{-\gamma t}. \]

Step2: 式 B.1 より $x_t - \mathbb{E}[x_t] = \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s$ である. 被積分関数 $e^{-\gamma(t-s)}$ は決定論的なので, 伊藤の等長性を用いると \[ \text{Var}(x_t) = \sigma^2\, \mathbb{E}\left[\left(\int_0^t e^{-\gamma (t - s)}\, dB_s\right)^{2}\right] = \sigma^2 \int_0^t e^{-2\gamma (t - s)}\, ds. \] $\tau = t - s$ と置換すると $\displaystyle \int_0^t e^{-2\gamma (t - s)}\, ds = \int_0^t e^{-2\gamma \tau}\, d\tau = \frac{1 - e^{-2\gamma t}}{2\gamma}$ なので, \[ \text{Var}(x_t) = \frac{\sigma^2}{2\gamma}\left(1 - e^{-2\gamma t}\right). \]

Step 3: 決定論的な被積分関数に対する伊藤積分は正規分布に従うので, 式 B.1 より $x_t$ も正規分布である. $t \to \infty$ で平均は $\mu$, 分散は $\dfrac{\sigma^2}{2\gamma}$ に収束するから, $x_t$ は定常分布 $\mathcal{N}\left(\mu,\, \dfrac{\sigma^2}{2\gamma}\right)$ に収束する.

B.2 最大値原理

B.2.1 HJB方程式 (Deterministic)

復習として, ノイズのない決定論的な最適制御問題の Hamilton-Jacobi-Bellman (HJB) 方程式を導きます.

命題 B.5 (HJB 方程式 (決定論的)) 状態 $x_t$ が制御 $u_t$ のもとで \[ \dot{x}_t = g(x_t, u_t) \] に従い, 意思決定者は割引率 $\rho$ で \[ V(x_0) = \max_{\{u_t\}} \int_0^{\infty} e^{-\rho t} f(x_t, u_t)\, dt \] を最大化するとする. このとき, 最適価値関数 $V(x)$ は次のHJB方程式を満たす. \[ \rho V(x) = \max_{u}\left\{ f(x, u) + V'(x)\, g(x, u) \right\}. \tag{B.2}\]

証明. 微小区間 $[t, t + dt]$ に最適性原理 (principle of optimality) を適用する. いま $dt$ だけ制御 $u$ を選び, 以降を最適に続けると \[ V(x) = \max_{u}\left\{ f(x, u)\, dt + e^{-\rho\, dt}\, V(x_{t+dt}) \right\}. \] 決定論的なので $x_{t+dt} = x + g(x, u)\, dt + o(dt)$ であり, テイラー展開より \[ V(x_{t+dt}) = V(x) + V'(x)\, g(x, u)\, dt + o(dt). \] これに $e^{-\rho\, dt} \approx 1 - \rho\, dt$ を代入し, $dt$ の1次まで残すと \[ V(x) = \max_{u}\left\{ V(x) + \left(f(x, u) + V'(x)\, g(x, u) - \rho V(x)\right) dt \right\} + o(dt). \] 両辺から $V(x)$ を引いて $dt$ で割り, $dt \to 0$ とすると式 B.2 を得る.

ここに確率的なノイズを加えた場合, HJB方程式を導くには伊藤の補題が必要になります.

B.2.2 伊藤の補題

補題 B.1 (伊藤の補題 (Itô’s lemma)) $x_t$ が確率微分方程式 $dx_t = g_t\, dt + \sigma_t\, dB_t$ ($g_t$ はドリフト, $\sigma_t$ は拡散係数で, ともに時刻 $t$ までの情報に依存してよい) に従うとき, 2階連続微分可能な関数 $\phi$ に対して \[ \begin{aligned} d\phi(x_t) &= \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, (dx_t)^2 \\ &= \left[\phi'(x_t)\, g_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\right] dt + \phi'(x_t)\, \sigma_t\, dB_t. \end{aligned} \]

証明. $\phi(x_t)$ の微小時間 $[t, t + dt]$ での変化をテイラー展開し, 2次の項まで残すと \[ \phi(x_{t+dt}) - \phi(x_t) = \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, (dx_t)^2 + o\left((dx_t)^2\right). \] ここに $dx_t = g_t\, dt + \sigma_t\, dB_t$ を代入すると \[ (dx_t)^2 = g_t^2 (dt)^2 + 2 g_t \sigma_t\, dt\, dB_t + \sigma_t^2 (dB_t)^2. \] 命題 B.1 の二次変分 $(dB_t)^2 = dt$ を用い, より高次の $(dt)^2$ と $dt\, dB_t$ を $o(dt)$ として落とすと $(dx_t)^2 = \sigma_t^2\, dt$ となる. これを代入すると \[ \begin{aligned} d\phi(x_t) &= \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\, dt \\ &= \left[\phi'(x_t)\, g_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\right] dt + \phi'(x_t)\, \sigma_t\, dB_t. \end{aligned} \]

B.2.3 HJB方程式 (Stochastic)

状態が確率微分方程式 $dx_t = g(x_t, u_t)\, dt + \sigma(x_t, u_t)\, dB_t$ に従い, 割引率 $\rho$ で \[ V(x_0) = \max_{\{u_t\}} \mathbb{E}_0 \int_0^{\infty} e^{-\rho t} f(x_t, u_t)\, dt \] を最大化する確率的な問題を考えます. $V(x)$ を最適価値関数と呼びます. 伊藤の補題を使うと, 決定論的な式 B.2 が次のように拡張されます.

命題 B.6 (HJB 方程式 (確率的)) 状態が $dx_t = g(x_t, u_t)\, dt + \sigma(x_t, u_t)\, dB_t$ に従うとき, 最適価値関数 $V$ は次の Hamilton-Jacobi-Bellman (HJB) 方程式を満たす. \[ \rho V(x) = \max_{u}\left\{ f(x, u) + V'(x)\, g(x, u) + \frac{1}{2}\,\sigma(x, u)^2\, V''(x) \right\}. \]

証明. 微小区間 $[t, t + dt]$ に最適性原理 (principle of optimality) を適用する. いま $dt$ だけ制御 $u$ を選び以降を最適に続けると, 状態が確率的なので将来の価値は期待値で評価される. \[ V(x) = \max_{u}\left\{ f(x, u)\, dt + e^{-\rho\, dt}\, \mathbb{E}\left[V(x_{t+dt}) \mid x_t = x\right] \right\}. \] 伊藤の補題より $dV = \left[V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x)\right] dt + V'(x)\, \sigma(x, u)\, dB_t$ であり, 伊藤積分の項は平均ゼロなので, \[ \mathbb{E}\left[V(x_{t+dt}) \mid x_t = x\right] = V(x) + \left[V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x)\right] dt. \] これと $e^{-\rho\, dt} \approx 1 - \rho\, dt$ を代入し, $dt$ の1次まで残すと \[ V(x) = \max_{u}\left\{ V(x) + \left(f(x, u) + V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x) - \rho V(x)\right) dt \right\} + o(dt). \] 両辺から $V(x)$ を引いて $dt$ で割り, $dt \to 0$ とすると命題の HJB 方程式を得る.

B.3 連続時間の Aiyagari モデル

ここからは, 上で整えた道具を使って, 連続時間の異質エージェントモデルを実際に解きます. 題材は Aiyagari モデル (3.3節) の連続時間版です. 生産性の対数を OU 過程でモデル化するので, HJB 方程式には伊藤の補題が生む2階微分の項がそのまま現れます. 数値解法は Achdou et al. (2022) の有限差分法です. 離散時間版と比べながら, 連続時間の数値計算がなぜ速く, なぜ頑健なのかを見ていきます.

B.3.1 設定

家計は割引率 $\rho$ で期待生涯効用 $\mathbb{E}_0 \int_0^\infty e^{-\rho t} u\left(c_t\right) dt$ を最大化します. 資産 $a$ は資本で, 予算制約と借入制約は \[ \dot{a}_t = w z_t + r a_t - c_t, \qquad a_t \geq 0 \] です. 労働生産性 $z_t = e^{x_t}$ の対数は OU 過程 \[ dx_t = -\theta x_t\, dt + \sigma\, dB_t \] に従います (命題 B.3 の記法で $\gamma = \theta$, $\mu = 0$ とし, 平均回帰の速さを CRRA 係数と区別して $\theta$ と書きます). OU 過程は AR(1) の連続時間版なので, 離散時間の章との対応は明快です. 1年間隔で観測すると自己相関 $e^{-\theta}$, 定常分散 $\sigma^2 / \left(2\theta\right)$ の AR(1) になるので (命題 B.4), 動的計画法の章のカリブレーション $\left(\rho_z, \sigma_\varepsilon\right) = \left(0.9, 0.2\right)$ をそのまま \[ \theta = -\log 0.9, \qquad \frac{\sigma^2}{2\theta} = \frac{\sigma_\varepsilon^2}{1 - \rho_z^2} \] で写し取れます. 割引率も同様に $\rho = -\log \beta = -\log 0.96$ とします. 企業は Cobb-Douglas 生産関数 $K^\alpha L^{1-\alpha}$ を持ち, $k = K / L$ として $r = \alpha k^{\alpha - 1} - \delta$, $w = \left(1 - \alpha\right) k^{\alpha}$ が成り立ちます. 労働供給は非弾力的で, $\mathbb{E}\left[z\right] = 1$ に基準化すれば $L = 1$ です.

定義 B.4 (定常均衡の方程式系) 定常均衡は, 価値関数 $v\left(a, x\right)$, 密度 $g\left(a, x\right)$, 価格 $\left(r, w\right)$ の組であって次を満たすものである.

HJB 方程式: \[ \rho v = \max_{c}\ u\left(c\right) + \partial_a v \left(w e^x + ra - c\right) - \theta x\, \partial_x v + \frac{\sigma^2}{2} \partial_{xx} v. \tag{B.3}\]

Kolmogorov Forward (KF) 方程式 (定常): \[ 0 = -\partial_a \left[s\left(a, x\right) g\right] + \partial_x \left[\theta x\, g\right] + \frac{\sigma^2}{2} \partial_{xx} g, \tag{B.4}\] ここで $s\left(a, x\right) = w e^x + ra - c\left(a, x\right)$ は貯蓄関数, $c\left(a, x\right) = \left(u'\right)^{-1}\left(\partial_a v\right)$.

資本市場の清算: \[ \iint a\, g\left(a, x\right)\, da\, dx = K, \qquad r = \alpha \left(K / L\right)^{\alpha - 1} - \delta. \tag{B.5}\]

HJB 方程式が価値関数を「後ろ向き」に決め, その政策関数のもとで KF 方程式が分布を「前向き」に決め, 分布が価格 $r$ を決めて HJB に戻る, というループ構造です.² アルゴリズムの全体像を先に図で示しておきます (図 B.1). 以下の数値解法の節は, この図の各ブロックを上から順に実装していきます.

図 B.1: Equilibrium loop of the finite difference method

借入制約の入り方が離散時間と決定的に違います. 連続時間では制約は内点で決して binding にならないため, 一階条件 $u'\left(c\left(a, x\right)\right) = \partial_a v\left(a, x\right)$ が $a > 0$ のすべての点で等号で成り立ちます. 制約は HJB 方程式そのものではなく, 境界条件 \[ \partial_a v\left(0, x\right) \geq u'\left(w e^x\right) \tag{B.6}\] として入ります. これは一階条件を通じて $s\left(0, x\right) \geq 0$ を含意するので, 制約が破られることはありません. KF 方程式には境界条件すら不要です. 最適政策が制約を尊重している限り, 分布は自動的に $\left[0, \infty\right)$ に収まるからです. $x$ 方向はグリッドを定常標準偏差の $\pm 3$ 倍で打ち切り, 端では反射壁 (reflecting barrier) を仮定します. OU 過程のドリフトは端では常に内側を向くので, この打ち切りの影響はわずかです.

B.4 有限差分法

B.4.1 Upwind スキーム

資産グリッド $a_1 = 0 < a_2 < \cdots < a_{N_a}$ (等間隔, 幅 $\Delta a$) を取り, $x$ を固定して資産方向の微分を前方差分と後方差分で近似します. \[ \partial_a v_{i, F} = \frac{v_{i+1} - v_i}{\Delta a}, \qquad \partial_a v_{i, B} = \frac{v_i - v_{i-1}}{\Delta a}. \]

どちらを使うべきでしょうか. 答えは「状態が動く向きに合わせる」です. 貯蓄 (ドリフト) が正なら前方差分, 負なら後方差分を使う. これが upwind スキームで, $s_i^{+} = \max\left\{s_i, 0\right\}$, $s_i^{-} = \min\left\{s_i, 0\right\}$ として HJB 方程式の資産方向を \[ \frac{v_{i+1} - v_i}{\Delta a}\, s_i^{+} + \frac{v_i - v_{i-1}}{\Delta a}\, s_i^{-} \tag{B.7}\] と離散化します (図 B.2).

図 B.2: Upwind scheme: the difference follows the drift

有限差分スキームが真の解 (粘性解) に収束するための Barles–Souganidis の十分条件は単調性・整合性・安定性の3つで, このうち非自明なのが単調性です. Upwind は, まさに単調性を保証するための差分の選び方になっています (Achdou et al. 2022, sec. 5.3).

実装上は, まず前方差分で仮の貯蓄 $s_{i,F} = w e^{x} + r a_i - \left(u'\right)^{-1}\left(\partial_a v_{i,F}\right)$ を, 後方差分で $s_{i,B}$ を計算し,

$s_{i,F} > 0$ なら前方差分を採用する
$s_{i,B} < 0$ なら後方差分を採用する
どちらでもない場合 ($s_{i,F} \leq 0 \leq s_{i,B}$) はその点を定常点とみなして $s_i = 0$, $c_i = w e^{x} + r a_i$ とする

という規則にします.³

境界条件の扱いは驚くほど簡単です. 下端 $a_1 = 0$ では後方差分だけに状態制約の境界条件式 B.6 を等号で課します: $\partial_a v_{1, B} = u'\left(w e^{x}\right)$. 前方差分は普通に計算し, どちらを使うかは upwind に選ばせます. 貯蓄が正なら境界条件は使われず, 負になりそうなときだけ発動して $s = 0$ で止める, という仕掛けです. 上端では貯蓄が負になるので前方差分は使われず, 境界条件は実質不要です.

B.4.2 生産性方向の離散化

OU 過程の項 $-\theta x\, \partial_x v + \frac{\sigma^2}{2} \partial_{xx} v$ も同じ発想で離散化します. ドリフト $-\theta x$ は upwind (符号に応じて片側差分), 拡散項は中心差分です. グリッド $x_1 < \cdots < x_{N_x}$ (幅 $\Delta x$) 上の点 $j$ について, 隣へ移る「強度」をまとめると \[ \lambda_j^{\downarrow} = \frac{\sigma^2}{2 \Delta x^2} + \frac{\max\left\{\theta x_j, 0\right\}}{\Delta x}, \qquad \lambda_j^{\uparrow} = \frac{\sigma^2}{2 \Delta x^2} + \frac{\max\left\{-\theta x_j, 0\right\}}{\Delta x} \tag{B.8}\] となります ($x_j > 0$ なら平均回帰のドリフトは下向きに働く, という形です). 端のグリッドでは外向きの強度を落とします. これが反射壁にあたります. できあがる三重対角行列 $B_x$ は各行の和がゼロの強度行列 (連続時間マルコフ連鎖の生成作用素) です. つまりこの離散化は, AR(1) を有限状態マルコフ連鎖で近似した Tauchen 法・Rouwenhorst 法 (2.4節) の連続時間版です.⁴

B.4.3 遷移行列 $A$

以上を合わせると, 離散化された HJB 方程式は \[ \rho v = u + A\left(v\right) v \tag{B.9}\] という形になります. $A$ は $N_a N_x \times N_a N_x$ の疎行列で, 資産方向の upwind 部分 (三重対角) と生産性方向の $B_x$ (帯2本) の和です. 重要な性質が2つあります.

各行の和はゼロ: $A$ は全体として連続時間マルコフ連鎖の強度行列である
有限差分法は, 富と生産性の連続な動きを「隣のグリッド点へのポアソンジャンプ」で近似していると解釈できる. 離散時間で遷移確率行列を作ったことの, 連続時間版にあたる

図 B.3 がこのジャンプの解釈を図解しています. 状態 $\left(a_i, x_j\right)$ から確率的に移れる先は, 上下左右の隣接グリッドの4点だけです. 資産方向は upwind が選んだ片側 (実際に正の強度を持つのは矢印2本のうち1本), 生産性方向は式 B.8 の2本です.

図 B.3: Poisson-jump interpretation of the generator

この「隣にしか跳べない」構造が, そのまま $A$ の疎パターンになります. 図 B.4 は小さいグリッド ($N_a = 12$, $N_x = 4$) で $A$ の非ゼロ要素の位置を描いたものです. 資産方向のジャンプは各ブロック内の三重対角 (対角とその上下1本ずつ), 生産性方向のジャンプはブロックを跨ぐ $\pm N_a$ の位置の帯2本として現れ, $A$ の各行の非ゼロ要素は高々5個です.

Code

Nad, Nxd = 12, 4
Ta = spdiagm(-1 => ones(Nad - 1), 0 => fill(2.0, Nad), 1 => ones(Nad - 1))
Bd = spdiagm(-1 => ones(Nxd - 1), 0 => fill(2.0, Nxd), 1 => ones(Nxd - 1))
Ad = blockdiag((Ta for _ in 1:Nxd)...) + kron(Bd, sparse(I, Nad, Nad))
nzrows, nzcols, _ = findnz(Ad)
scatter(nzcols, nzrows, yflip=true, aspect_ratio=1, markersize=2.5, color=:black,
    markerstrokewidth=0, legend=false, xlabel="Column", ylabel="Row",
    xlims=(0, Nad * Nxd + 1), ylims=(0, Nad * Nxd + 1), size=(380, 380))
vline!(collect(Nad:Nad:(Nad*(Nxd-1))) .+ 0.5, color=:gray, ls=:dash, alpha=0.6)
hline!(collect(Nad:Nad:(Nad*(Nxd-1))) .+ 0.5, color=:gray, ls=:dash, alpha=0.6)

図 B.4: Sparsity pattern of $A$ ($N_a = 12$, $N_x = 4$ for illustration).

行列のサイズは $N_a N_x$ の2乗で増えますが, 非ゼロ要素は $5 N_a N_x$ 程度にしか増えません. 疎行列ソルバはゼロでない要素だけを保存・計算するので, 後で出てくる連立一次方程式 (陰的スキームと KF 方程式) は巨大なグリッドでも一瞬で解けます. これが有限差分法の速さの核心です.

B.4.4 陰的スキーム

式 B.9 は $A$ が $v$ に依存する非線形方程式なので, 反復で解きます. 擬似的な時間ステップ $\Delta$ を入れた陰的スキーム (implicit scheme) を使います. \[ \frac{v^{n+1} - v^n}{\Delta} + \rho v^{n+1} = u\left(c^n\right) + A\left(v^n\right) v^{n+1} \quad \Longleftrightarrow \quad \left[\left(\frac{1}{\Delta} + \rho\right) I - A^n\right] v^{n+1} = u^n + \frac{v^n}{\Delta}. \tag{B.10}\] 各反復は疎な連立一次方程式を1回解くだけです. 陰的スキームは無条件に安定なので, $\Delta$ をいくらでも大きく取れます ($\Delta = 1000$ など). その結果, 反復は数十回で収束します. 素朴な価値関数反復が何百回もかかるのと比べると大きな違いです (離散時間との定量的な比較は実装のあとで行います).

B.4.5 KF 方程式はタダで解ける

分布はどう求めるのでしょうか. ここが有限差分法の最も美しい部分です. KF 作用素は HJB に現れる生成作用素の随伴 (adjoint) なので, 離散化すると式 B.4 は \[ A^{\top} g = 0 \tag{B.11}\] になります (Achdou et al. 2022, sec. 5.4). HJB を解いた時点で $A$ は手元にあるので, KF 方程式は転置行列の連立一次方程式を1回解くだけです. 分布のシミュレーションも, 遷移行列の別途構築も要りません. 式 B.11 は特異 (解が定数倍で不定) なので, 1つの成分を固定して解いてから, 密度の積分が1になるように正規化します. 借入制約に張り付く家計の質量は, 下端のグリッド点に立つ大きな値として現れます.

B.4.6 均衡

利子率を与えると企業の一階条件から $k\left(r\right) = \left(\alpha / \left(r + \delta\right)\right)^{1 / \left(1 - \alpha\right)}$ と $w\left(r\right)$ が決まるので, 資本の超過供給関数 \[ S\left(r\right) = \iint a\, g\, da\, dx - k\left(r\right) L \] の二分法で均衡利子率を求めます. $r$ を与えて HJB → KF → $S\left(r\right)$ と計算し, 正なら $r$ を下げ, 負なら上げるだけです.

まとめると, この手法が速く頑健な理由は4つです (Achdou et al. 2022, sec. 5.1).

借入制約は境界条件に封じ込められ, 内点の一階条件は常に等号で成り立つ (occasionally binding な場合分けが消える)
一階条件が静学的で, $c = \left(\partial_a v\right)^{-1/\gamma}$ と手で解ける. 求根も期待値の計算も不要
連続時間の過程はグリッド上を1点ずつしか動かないので, 行列が疎になる
KF 方程式が HJB の転置としてタダで手に入る

B.5 Julia 実装

カリブレーションは動的計画法の章の Aiyagari モデルを OU 過程に写したものです.

const γ, ρ = 2.0, -log(0.96)      # CRRA, discount rate (= -log β)
const α, δ = 0.36, 0.08           # production
const θ = -log(0.9)               # OU mean reversion:  e^{-θ} = 0.9
const σx = 0.2 / sqrt(1 - 0.9^2)  # stationary sd of x (matches the AR(1))
const σ = σx * sqrt(2θ)           # OU diffusion coefficient

const amin, amax, Na = 0.0, 50.0, 2000
const a = range(amin, amax, length=Na)
const da = step(a)
const Nx = 31
const x = range(-3σx, 3σx, length=Nx)
const dx = step(x)

u(c) = c^(1 - γ) / (1 - γ)
u′(c) = c^(-γ)
u′inv(v) = max(v, 1e-10)^(-1 / γ)  # guard against negative slopes early on

まず生産性方向の生成作用素 $B_x$ を式 B.8 のとおり組み立てます. その定常分布を使って $\mathbb{E}\left[z\right] = 1$ に基準化し, 離散化が OU 過程の定常標準偏差を再現しているかも確認しておきます.

function ou_generator()
    B = spzeros(Nx, Nx)
    for j in 1:Nx
        lo = σ^2 / (2dx^2) + max(θ * x[j], 0) / dx    # intensity of moving down
        hi = σ^2 / (2dx^2) + max(-θ * x[j], 0) / dx   # intensity of moving up
        if j > 1
            ;
            B[j, j-1] = lo;
            B[j, j] -= lo;
        end # boundary: drop the outward jump
        if j < Nx
            ;
            B[j, j+1] = hi;
            B[j, j] -= hi;
        end
    end
    return B
end
const Bx = ou_generator()

function stat_x(B)                              # stationary distribution over x
    BT = copy(sparse(transpose(B)))
    b = zeros(Nx);
    BT[1, :] .= 0;
    BT[1, 1] = 1.0;
    b[1] = 0.1
    ψ = BT \ b
    ψ ./ sum(ψ)
end
const ψ = stat_x(Bx)
const z = exp.(x) ./ dot(ψ, exp.(x))                  # normalize E[z] = 1
const Bswitch = kron(Bx, sparse(I, Na, Na))     # x-part on the full (a, x) grid
(sd_discretized=round(sqrt(dot(ψ, x .^ 2) - dot(ψ, x)^2), digits=4), sd_theory=round(σx, digits=4))

(sd_discretized = 0.4669, sd_theory = 0.4588)

HJB 方程式のソルバは, upwind による政策の計算, 遷移行列 $A$ の構築, 陰的更新式 B.10 の3ステップの反復です.

function solve_hjb(r, w; Δ=1000.0, maxit=400, tol=1e-8)
    inc = w .* z' .+ r .* a                  # inc[i, j] = w z_j + r a_i
    v = u.(inc) ./ ρ
    dVf = similar(v);
    dVb = similar(v)
    A = spzeros(Na * Nx, Na * Nx);
    c = similar(v)
    for it in 1:maxit
        dVf[1:(end-1), :] = diff(v, dims=1) ./ da
        dVf[end, :] = u′.(inc[end, :])       # unused: saving is negative up top
        dVb[2:end, :] = diff(v, dims=1) ./ da
        dVb[1, :] = u′.(inc[1, :])         # state-constraint boundary condition
        cf = u′inv.(dVf);
        sf = inc .- cf     # saving under the forward difference
        cb = u′inv.(dVb);
        sb = inc .- cb     # saving under the backward difference
        If = sf .> 0
        Ib = (sb .< 0) .& .!If
        # where neither drift applies the household is at rest, so s = 0
        c = cf .* If .+ cb .* Ib .+ inc .* (.!If .& .!Ib)

        X = -min.(sb, 0) ./ da         # intensity of moving one grid point down
        Z = max.(sf, 0) ./ da            # intensity of moving one grid point up
        A = blockdiag(
            (spdiagm(-1 => X[2:end, j], 0 => -X[:, j] .- Z[:, j], 1 => Z[1:(end-1), j])
             for j in 1:Nx)...) + Bswitch

        B = (1 / Δ + ρ) * sparse(I, Na * Nx, Na * Nx) - A
        vnew = reshape(B \ (vec(u.(c)) .+ vec(v) ./ Δ), Na, Nx)
        maximum(abs.(vnew .- v)) < tol && return (v=vnew, c=c, s=inc .- c, A=A, iter=it)
        v = vnew
    end
    error("HJB did not converge")
end

KF 方程式は転置行列の連立一次方程式です. 特異性は1成分を固定して外し, 最後に密度として正規化します.

function stationary_dist(A)
    AT = copy(sparse(transpose(A)))
    b = zeros(Na * Nx)
    AT[1, :] .= 0.0;
    AT[1, 1] = 1.0;
    b[1] = 0.1  # fix one entry of the singular system
    g = AT \ b
    g ./= sum(g) * da * dx                        # normalize to a density
    return reshape(g, Na, Nx)
end

最後に, 資本の超過供給に対する二分法で均衡利子率を求めます.

prices(r) = (k=(α / (r + δ))^(1 / (1 - α)); (k=k, w=(1 - α) * k^α))

function excess_supply(r)
    (; k, w) = prices(r)
    sol = solve_hjb(r, w)
    g = stationary_dist(sol.A)
    return sum(g .* a) * da * dx - k, sol, g
end

function equilibrium(rlo, rhi; n=20)
    for _ in 1:n
        rm = (rlo + rhi) / 2
        first(excess_supply(rm)) > 0 ? (rhi = rm) : (rlo = rm)
    end
    return (rlo + rhi) / 2
end

r★ = equilibrium(0.005, ρ - 1e-4)
ex, sol, g = excess_supply(r★)
(r=round(r★, digits=4), K=round(prices(r★).k, digits=2), hjb_iter=sol.iter)

(r = 0.0243, K = 6.93, hjb_iter = 43)

均衡利子率は割引率 $\rho = -\log 0.96 \approx 0.041$ を明確に下回ります. 予備的貯蓄による過剰な資本蓄積という Aiyagari モデルの中心的な結果が, 連続時間でもそのまま現れています. 離散時間の章で得た $r^* \approx 0.021$ とはカリブレーションの写像や離散化の違いから小さなずれがありますが, 定性的な結論は同じです. そして, グリッドが $N_a N_x = 62{,}000$ 点あるにもかかわらず, HJB 方程式の反復が数十回, 均衡全体も二分法 20 回で数分以内に収まっていることに注目してください. 速さの源泉は図 B.4 で見た疎性です. 実際の $A$ は $62{,}000 \times 62{,}000$ ですが, 非ゼロ要素は全体の 0.01% にも遠く及びません.

(nonzeros=nnz(sol.A), share_percent=round(100 * nnz(sol.A) / (Na * Nx)^2, sigdigits=2))

(nonzeros = 243957, share_percent = 0.0063)

B.5.1 政策関数

Code

jmid = (Nx + 1) ÷ 2
idx = a .<= 10.0
plot(a[idx], sol.s[idx, 1], label=L"x = -3\sigma_x", lw=2)
plot!(a[idx], sol.s[idx, jmid], label=L"x = 0", lw=2, ls=:dash)
plot!(a[idx], sol.s[idx, Nx], label=L"x = 3\sigma_x", lw=2, ls=:dashdot)
hline!([0], color=:gray, ls=:dot, label=false)
plot!(xlabel=L"a", ylabel="Saving", legend=:topright)

図 B.5: Saving policy functions. Low-productivity households at the constraint have $s = 0$.

生産性の低い家計は資産を取り崩しながら借入制約に向かって滑り落ち, 制約点でちょうど $s = 0$ になります. 内点で制約が binding にならない, という理論の構造がそのまま見えています. 生産性の高い家計は広い範囲で貯蓄を続けますが, 十分に資産を持つと ($\rho > r^*$ なので) やはり取り崩しに転じます.

B.5.2 定常分布

Code

ga = vec(sum(g, dims=2)) .* dx                # marginal density of assets
plot(a[2:end][a[2:end] .<= 30], ga[2:end][a[2:end] .<= 30],
    label=false, lw=2, xlabel=L"a", ylabel="Density")

図 B.6: Stationary asset distribution. The point mass at the constraint is omitted.

println("mass at the borrowing constraint: ", round(sum(g[1, :]) * da * dx, sigdigits=3))

mass at the borrowing constraint: 0.0322

離散時間版と同じく, 借入制約の近くに密度が集中し, 制約点そのものには質量が立ちます (図では質量を除いた密度部分を描いています). KF 方程式に一切の境界条件を与えていないのに, 質量まで含めて正しい分布が「転置を解くだけ」で出てくるのが, この手法の強みです.

B.5.3 離散時間モデルとの比較

連続時間版が本当に離散時間の Aiyagari モデルを再現しているのか, 同じカリブレーションの離散時間モデルを解いて確かめます. 実装は動的計画法の章と同じ, 単調性を利用したグリッドサーチの VFI に Howard の改善を加えたものです. 所得過程の離散化には Rouwenhorst 法を使います. Tauchen 法は少ない状態数で分散を過大評価するため (2.4節), 比較がその分だけ歪んでしまうからです. Rouwenhorst 法は AR(1) の2次モーメントを厳密に再現するので, 離散時間の所得過程は連続時間の OU 過程と (2次モーメントの意味で) 完全に揃います.

Code

const β = exp(-ρ)                        # discrete-time discount factor

function dt_setup(; na, nz, uniform=false)
    mc = rouwenhorst(nz, 0.9, 0.2)
    Π = mc.p
    zdt = exp.(collect(mc.state_values))
    zdt ./= dot(vcat(stationary_distributions(mc)...), zdt)   # E[z] = 1
    ag = uniform ? collect(range(0, amax, length=na)) :
         collect(range(0, 1, length=na)) .^ 2 .* amax
    (; Π, zdt, ag, na, nz)
end

function dt_vfi(r, w, m; howard=50, tol=1e-6, maxit=1000)
    (; Π, zdt, ag, na, nz) = m
    V = [u(w * zdt[iz] + r * ag[ia] + 1e-10) / (1 - β) for ia in 1:na, iz in 1:nz]
    pol = ones(Int, na, nz)
    iters = 0
    for it in 1:maxit
        iters = it
        EV = β .* (V * Π')
        Vnew = similar(V)
        for iz in 1:nz
            lo = 1
            for ia in 1:na
                res = w * zdt[iz] + (1 + r) * ag[ia]
                best = -Inf
                besti = lo
                for ja in lo:na                       # monotone grid search
                    c = res - ag[ja]
                    c <= 0 && break
                    val = u(c) + EV[ja, iz]
                    if val > best
                        best = val
                        besti = ja
                    else
                        break                       # concavity: passed the peak
                    end
                end
                Vnew[ia, iz] = best
                pol[ia, iz] = besti
                lo = besti
            end
        end
        for _ in 1:howard                             # Howard improvement
            EVh = β .* (Vnew * Π')
            for iz in 1:nz, ia in 1:na
                c = w * zdt[iz] + (1 + r) * ag[ia] - ag[pol[ia, iz]]
                Vnew[ia, iz] = u(c) + EVh[pol[ia, iz], iz]
            end
        end
        if maximum(abs.(Vnew .- V)) < tol
            V = Vnew
            break
        end
        V = Vnew
    end
    return V, pol, iters
end

function dt_dist(pol, m; tol=1e-10, maxit=100_000)
    (; Π, na, nz) = m
    n = na * nz
    Is = Int[];
    Js = Int[];
    Vs = Float64[]
    for iz in 1:nz, ia in 1:na
        row = ia + (iz - 1) * na
        for jz in 1:nz
            push!(Is, row)
            push!(Js, pol[ia, iz] + (jz - 1) * na)
            push!(Vs, Π[iz, jz])
        end
    end
    Tt = sparse(Js, Is, Vs, n, n)                 # transposed transition matrix
    λ = fill(1 / n, n)
    iters = 0
    for it in 1:maxit
        iters = it
        λn = Tt * λ
        if maximum(abs.(λn .- λ)) < tol
            λ = λn
            break
        end
        λ = λn
    end
    return reshape(λ, na, nz), iters
end

dt_assets(λ, m) = sum(λ .* m.ag)

まず, 資本の需要曲線と供給曲線を両方のモデルで重ねてみます (図 B.7).

Code

mdt = dt_setup(na=200, nz=7)

function dt_supply(r)
    w = prices(r).w
    _, pol, _ = dt_vfi(r, w, mdt)
    dt_assets(first(dt_dist(pol, mdt)), mdt)
end
ct_supply(r) = first(excess_supply(r)) + prices(r).k

rvals = range(0.005, 0.033, length=8)
A_dt = dt_supply.(rvals)
A_ct = ct_supply.(rvals)

rfine = range(0.005, 0.033, length=100)
plot(map(r -> prices(r).k, rfine), rfine, label=L"Demand $k(r)$",
    color=:gray, ls=:dot, lw=2)
plot!(A_ct, rvals, label="Supply (continuous time)", lw=2, marker=:circle)
plot!(A_dt, rvals, label="Supply (discrete time)", lw=2, ls=:dash, marker=:square)
plot!(xlabel="Assets / capital", ylabel=L"r", legend=:bottomright)

図 B.7: Capital demand and supply. Continuous vs. discrete time, same calibration.

供給曲線はほぼ重なり, 需要曲線との交点 (均衡利子率) も互いに近い位置にあります.

function dt_equilibrium(rlo, rhi; n=20)
    for _ in 1:n
        rm = (rlo + rhi) / 2
        dt_supply(rm) > prices(rm).k ? (rhi = rm) : (rlo = rm)
    end
    return (rlo + rhi) / 2
end
(r_discrete=round(dt_equilibrium(0.005, 0.035), digits=4), r_continuous=round(r★, digits=4))

(r_discrete = 0.025, r_continuous = 0.0243)

残る小さな差は, 意思決定のタイミング (離散時間では1年分まとめて意思決定する) と離散化そのものの違いによるものです. 連続時間版は離散時間の Aiyagari モデルの答えを再現しつつ, 滑らかな政策関数と分布を返している, と言ってよいでしょう.

次に計算時間です. 同じ均衡価格のもとで, 家計問題と定常分布の計算をそれぞれ計測します (収束判定はどちらも $\max \left|\Delta v\right| < 10^{-6}$).

mdt31 = dt_setup(na=Na, nz=Nx, uniform=true)   # same grid as continuous time
w★ = prices(r★).w
dt_vfi(r★, w★, mdt31)                                 # warmup

t_vfi7 = @elapsed ((_, pol7, it7) = dt_vfi(r★, w★, mdt))
t_dist7 = @elapsed ((λ7, itd7) = dt_dist(pol7, mdt))
t_vfi31 = @elapsed ((_, pol31, it31) = dt_vfi(r★, w★, mdt31))
t_dist31 = @elapsed ((λ31, itd31) = dt_dist(pol31, mdt31))
t_hjb = @elapsed (solb = solve_hjb(r★, w★; tol=1e-6))
t_kf = @elapsed stationary_dist(solb.A)

cell(t, n, unit) = string(round(t, digits=2), " 秒 / ", n, " ", unit)

表 B.1: 家計問題と定常分布の計算時間

benchrows = [
    "離散時間 VFI + Howard" raw"$200 \times 7$" cell(t_vfi7, it7, "反復") cell(t_dist7, itd7, "反復");
    "離散時間 VFI + Howard" raw"$2{,}000 \times 31$" cell(t_vfi31, it31, "反復") cell(t_dist31, itd31, "反復");
    "連続時間 有限差分" raw"$2{,}000 \times 31$" cell(t_hjb, solb.iter, "反復") cell(t_kf, 1, "回の求解")
]
pretty_table(benchrows;
    column_labels=["手法", raw"状態数 $N_a \times N_z$", "家計問題", "定常分布"],
    backend=:markdown, allow_markdown_in_cells=true)

手法	状態数 $N_a \times N_z$	家計問題	定常分布
離散時間 VFI + Howard	$200 \times 7$	0.01 秒 / 13 反復	0.01 秒 / 616 反復
離散時間 VFI + Howard	$2{,}000 \times 31$	0.32 秒 / 17 反復	1.03 秒 / 591 反復
連続時間有限差分	$2{,}000 \times 31$	4.37 秒 / 31 反復	0.29 秒 / 1 回の求解

結果は一見すると意外かもしれません. 家計問題単体では, チューニングされた離散時間の VFI も十分に速いのです. これは選択変数が1次元で, 政策の単調性と価値関数の凹性を使うとグリッドサーチがほぼ線形時間で済むという, この問題に固有の事情によります. 連続時間の優位がはっきり出るのは定常分布です. 離散時間では遷移行列の反復 (べき乗法) を数百回まわす必要があるのに対し, 連続時間では HJB の副産物 $A$ の転置を1回解くだけで済みます.

より本質的な違いは, 表の数字の外にあります.

離散時間の速さは「1次元の選択を単調なグリッドサーチで解ける」ことに依存する. 選択変数が増える (ポートフォリオ選択, 耐久財など) と最大化のコストが急速に増すのに対し, 連続時間の一階条件は静学的なまま変わらない
連続時間の政策関数はグリッドへの丸めがなく滑らかで, 借入制約もちょうど1点で処理される. 推定 (4.2節) のように解を何千回も呼ぶ用途では, この滑らかさが目的関数の滑らかさに直結する

B.6 まとめと拡張

この章のアルゴリズムをまとめると次のとおりです.

$r$ を所与として企業の一階条件から $\left(k, w\right)$ を求め, upwind + 陰的スキームで HJB 方程式式 B.9 を解く (反復数十回)
遷移行列の転置から KF 方程式式 B.11 を1発で解き, 定常分布を得る
資本市場の超過供給を見て $r$ を二分法で更新する

同じ骨格が, ほとんど変更なしに広いクラスのモデルへ拡張できます (Achdou et al. 2022, sec. 7).

より単純な所得過程: 所得を2状態のポアソン過程に落とすと2階微分の項が消えて HJB は常微分方程式になり, Achdou et al. (2022) の本文で解析される Huggett モデルになる. 借入制約近傍の貯蓄関数や資産分布の質量について, 鋭い理論的特徴づけが得られる
移行動学: 時間依存の HJB を終端条件から後ろ向きに, KF を初期分布から前向きに解き, 価格の経路を収束するまで更新する. MIT ショックの分析がそのまま可能
非凸性・停止問題: 住宅のような不可分性や, デフォルトのような停止問題も, 粘性解の枠組みのもとで同じ離散化で扱える

離散時間の動的計画法 (第3章) と対をなす計算技術として, この有限差分法は現代の異質エージェント・マクロの標準的な道具になっています. 実装のバリエーションは Achdou et al. (2022) の公開コードが網羅しており, この章のコードはその最小構成にあたります.

$L^2$ 収束は分布収束を含意する. $\mathcal{N}(0, v_n)$ の特性関数 $e^{-v_n u^2 / 2}$ は $v_n \to v$ より $\mathcal{N}(0, v)$ の特性関数 $e^{-v u^2 / 2}$ に各点収束するので, レヴィの連続性定理 (Lévy’s continuity theorem) により極限は $\mathcal{N}(0, v)$ に従う.↩︎
この「後ろ向きの HJB + 前向きの KF」という連立系は, 平均場ゲーム (Mean Field Game) の定常系そのものです. Achdou et al. (2022) の枠組みは, Lasry と Lions の MFG 理論を異質エージェント・マクロに持ち込んだものといえます.↩︎
なお Achdou et al. (2022) の本文はスキームの詳細を数値付録に委ねているため, この規則と後述の陰的スキームは, 公開されている数値付録と実装コードに従っています.↩︎
upwind のドリフト近似は分散をわずかに水増しします. Tauchen 法で見た「升目への丸めが余分な分散を生む」現象と同種のもので, $\Delta x \to 0$ で消えます. 下の実装では離散化された定常分布の標準偏差を理論値と比べて確認します.↩︎

--- engine: julia --- # 連続時間モデル {#sec-apdx-ct} ```{julia} #| label: setup #| output: false using LinearAlgebra using SparseArrays using Plots using LaTeXStrings using QuantEcon using PrettyTables ``` ```{julia} #| label: graphics #| include: false mypalette = ["#d5968c", "#c2676d", "#5c363a", "#995041", "#45939c", "#0f6a81"] default(size=(500, 308), titlefontsize=10, fmt=:svg, grid=false, foreground_color_legend=nothing, palette=mypalette[[6, 2, 4, 3, 1, 5]]) ``` 本文の章はすべて離散時間でモデルを組んできましたが, 現代のマクロ経済学では連続時間の定式化も広く使われています. この章の前半では, 連続時間の確率過程 (ブラウン運動, 伊藤積分, OU 過程) と最大値原理 (HJB 方程式) という理論の道具を整えます. 後半では, その道具を使って Aiyagari モデルの連続時間版を**有限差分法** (finite difference method) で解きます [@achdou2022]. ## 確率過程連続時間の確率過程は, 時間 $t$ に対してランダムな値をとる関数 $x_t$ の族 $(x_t)_{t \geq 0}$ として定義されます. 連続時間の確率過程として基礎になるのは, 連続的なノイズを表すブラウン運動 (Brownian motion) です. ### ブラウン運動 ::: {#def-brownian-motion} ## ブラウン運動 $b_0 = 0$ を満たす確率過程 $(b_t)_{t \geq 0}$ は次の性質を満たすとき, ブラウン運動と呼ばれる. - 独立増分: 任意の $0 \leq t_1 < t_2 < \cdots < t_n$ に対して, $b_{t_2} - b_{t_1}, b_{t_3} - b_{t_2}, \ldots, b_{t_n} - b_{t_{n-1}}$ は独立である. - 正規性: 任意の $t$ に対して, $b_t \sim \mathcal{N}(\mu t, \sigma^2 t)$ - 連続パス: $b_t$ は $t$ の関数として連続である. 特に, $\mu = 0, \sigma = 1$ の場合を標準ブラウン運動と呼ぶ. ::: なお標準ブラウン運動 $B_t$ は明らかに $B_{t+1} - B_t \sim \mathcal{N}(0, 1)$ を満たすので, 離散時間のホワイトノイズの連続時間版とみなすことができます. 後で伊藤積分や伊藤の補題で繰り返し使う $(dB_t)^2 = dt$ という計算規則は, ブラウン運動の二次変分 (quadratic variation) が時間 $t$ に等しいという性質に由来します. ::: {#prp-bm-qv} ## ブラウン運動の二次変分標準ブラウン運動 $B_t$ に対して, 区間 $[0, t]$ の分割 $0 = s_0 < s_1 < \cdots < s_n = t$ を考え, その幅 (メッシュ) を $\|P\| := \max_i (s_{i+1} - s_i)$ とする. メッシュを $0$ に近づけると, 増分の二乗和は $t$ に二乗平均収束する. $$ \sum_{i=0}^{n-1} \left(B_{s_{i+1}} - B_{s_i}\right)^2 \;\xrightarrow{\ L^2\ }\; t \qquad (\|P\| \to 0). $$ これを微分の記法で略記したものが $(dB_t)^2 = dt$ である. ::: ::: {.proof} 二乗和を $S := \sum_{i=0}^{n-1} \Delta_i^2$ (ただし $\Delta_i := B_{s_{i+1}} - B_{s_i}$, 区間幅 $\delta_i := s_{i+1} - s_i$) と書く. $S$ が $t$ に二乗平均収束 ($L^2$ 収束) するとは, 定義により $$ \mathbb{E}\left[(S - t)^2\right] \longrightarrow 0 \qquad (\|P\| \to 0) $$ が成り立つことである. これを示す. ブラウン運動の正規性と独立増分性から $\Delta_i \sim \mathcal{N}(0, \delta_i)$ で互いに独立である. Step1: $\mathbb{E}[\Delta_i^2] = \text{Var}(\Delta_i) = \delta_i$ なので, 区間幅の総和が $t$ であることから $$ \mathbb{E}[S] = \sum_i \delta_i = t. $$ Step2: $\Delta_i^2$ は独立なので分散は和になる. 正規分布 $\mathcal{N}(0, \delta_i)$ は $\mathbb{E}[\Delta_i^4] = 3\delta_i^2$ を満たすから, $\text{Var}(\Delta_i^2) = 3\delta_i^2 - \delta_i^2 = 2\delta_i^2$ であり, $$ \text{Var}(S) = \sum_i 2\delta_i^2 \le 2 \left(\max_i \delta_i\right) \sum_i \delta_i = 2 t\, \|P\|. $$ Step3: $\mathbb{E}[S] = t$ なので二乗平均誤差は分散に一致し, $$ \mathbb{E}\left[(S - t)^2\right] = \text{Var}(S) \le 2 t\, \|P\| \xrightarrow{\ \|P\| \to 0\ } 0. $$ よって $S \xrightarrow{\ L^2\ } t$ が成り立つ. ::: ### 伊藤積分ここからブラウン運動を含む確率微分方程式を扱います. このとき, 伊藤積分を用います. ::: {#def-ito-integral} ## 伊藤積分過去の情報のみに依存する関数 $f_s$ と標準ブラウン運動 $B_t$ に対して, 伊藤積分 $\int_0^t f_s\, dB_s$ を, 区間分割 $0 = s_0 < s_1 < \cdots < s_n = t$ に沿った和の二乗平均極限として定義する. $$ \int_0^t f_s\, dB_s := \lim_{\max_i |s_{i+1} - s_i| \to 0} \sum_{i=0}^{n-1} f_{s_i}\,\left(B_{s_{i+1}} - B_{s_i}\right). $$ このとき次の2つの性質が成り立つ. - 平均ゼロ: $\displaystyle \mathbb{E}\left[\int_0^t f_s\, dB_s\right] = 0$ - 伊藤の等長性 (Itô isometry): $\displaystyle \mathbb{E}\left[\left(\int_0^t f_s\, dB_s\right)^2\right] = \int_0^t \mathbb{E}\left[f_s^2\right]\, ds$ ::: 各増分 $B_{s_{i+1}} - B_{s_i}$ を区間の左端 $s_i$ の値 $f_{s_i}$ で評価する点が重要です. これにより, $f_{s_i}$ は区間 $[s_i, s_{i+1}]$ の将来の増分 $B_{s_{i+1}} - B_{s_i}$ に依存しない (独立である) という性質が使えます. この独立性を踏まえると, 平均ゼロは直感的です. $$ \mathbb{E}\left[f_{s_i} \Delta_i\right] = \mathbb{E}[f_{s_i}]\, \mathbb{E}[\Delta_i] = 0. $$ 等長性に関しては, 和の二乗を展開すると, 異なる区間の積 ($i \neq j$) は独立かつ平均ゼロより消え, 対角項だけが残ります. 各対角項は $\mathbb{E}[f_{s_i}^2\, \Delta_i^2] = \mathbb{E}[f_{s_i}^2]\, \mathbb{E}[\Delta_i^2]$ となるので, @prp-bm-qv の二次変分を使うと, $$ \sum_i \mathbb{E}[f_{s_i}^2\, \Delta_i^2] = \sum_i \mathbb{E}[f_{s_i}^2]\, (s_{i+1} - s_i) \xrightarrow{\ n \to \infty\ } \int_0^t \mathbb{E}[f_s^2]\, ds. $$ 微分の記法でいえば, $\mathbb{E}[dB_t] = 0$ が平均ゼロを, $(dB_t)^2 = dt$ が等長性を生んでいる, と見ることができます. ::: {#prp-ito-normal} ## 伊藤積分の正規性被積分関数 $f_s$ が決定論的な ($B_s$ に依存しない) とき, 伊藤積分は正規分布に従う. $$ \int_0^t f_s\, dB_s \;\sim\; \mathcal{N}\left(0,\ \int_0^t f_s^2\, ds\right). $$ ::: ::: {.proof} 伊藤積分は近似和 $\sum_i f_{s_i}\, \Delta_i$ ($\Delta_i := B_{s_{i+1}} - B_{s_i}$) の $L^2$ 極限である. $f_s$ が決定論的なら $f_{s_i}$ は定数なので, 各項は $f_{s_i}\, \Delta_i \sim \mathcal{N}\left(0,\, f_{s_i}^2 (s_{i+1} - s_i)\right)$ で正規である. さらにブラウン運動の独立増分性から $\Delta_i$ は互いに独立なので, 独立な正規確率変数の和である近似和もまた正規分布に従う. $$ \sum_i f_{s_i}\, \Delta_i \;\sim\; \mathcal{N}\left(0,\ \sum_i f_{s_i}^2 (s_{i+1} - s_i)\right). $$ この近似和は $\|P\| \to 0$ で伊藤積分に $L^2$ 収束し, その分散は $\sum_i f_{s_i}^2 (s_{i+1} - s_i) \to \int_0^t f_s^2\, ds$ に収束する. 正規分布の族は分布収束に関して閉じている^[$L^2$ 収束は分布収束を含意する. $\mathcal{N}(0, v_n)$ の特性関数 $e^{-v_n u^2 / 2}$ は $v_n \to v$ より $\mathcal{N}(0, v)$ の特性関数 $e^{-v u^2 / 2}$ に各点収束するので, レヴィの連続性定理 (Lévy's continuity theorem) により極限は $\mathcal{N}(0, v)$ に従う.] ので, 平均ゼロ性と合わせて, 極限である伊藤積分も正規分布 $\mathcal{N}\left(0,\, \int_0^t f_s^2\, ds\right)$ に従う. ::: ### Orstein-Uhlenbeck 過程 AR(1) 過程の連続時間版として, Orstein-Uhlenbeck (OU) 過程を考えます. OU過程は以下の確率微分方程式で定義されます. ::: {#def-ou-process} ## OU過程 OU過程 $x_t$ は以下の確率微分方程式で定義される. $$ dx_t = -\gamma (x_t - \mu) dt + \sigma dB_t. $$ ここで $B_t$ は標準ブラウン運動. ::: ::: {#prp-ou-solution} ## OU過程の解 OU過程 $x_t$ は次のように陽に解ける. $$ x_t = \mu + (x_0 - \mu) e^{-\gamma t} + \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s. $${#eq-ou-solution} ::: ::: {.proof} 確率微分方程式を積分因子 (integrating factor) $e^{\gamma t}$ を用いて解く. 積分因子は決定論的なので, 通常の積の微分法則が使えて $$ d\left(e^{\gamma t}(x_t - \mu)\right) = \gamma e^{\gamma t}(x_t - \mu)\, dt + e^{\gamma t}\, dx_t. $$ ここに $dx_t = -\gamma (x_t - \mu)\, dt + \sigma\, dB_t$ を代入するとドリフト項が打ち消し合い, $$ d\left(e^{\gamma t}(x_t - \mu)\right) = e^{\gamma t}\, \sigma\, dB_t. $$ 両辺を $0$ から $t$ まで積分する. 左辺は微分を積分した増分 $\left[e^{\gamma s}(x_s - \mu)\right]_0^t$ なので, $$ e^{\gamma t}(x_t - \mu) - (x_0 - \mu) = \sigma \int_0^t e^{\gamma s}\, dB_s. $$ 両辺を $e^{\gamma t}$ で割ると, 解が陽に求まる. $$ x_t = \mu + (x_0 - \mu)\, e^{-\gamma t} + \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s. $$ ::: ::: {#prp-ou-moment} ## OU過程の平均と分散 $$ \mathbb{E}[x_t] = \mu + (x_0 - \mu)e^{-\gamma t}, \quad \text{Var}(x_t) = \frac{\sigma^2}{2\gamma} \left(1 - e^{-2\gamma t}\right). $$ また, $t \to \infty$ のとき, $x_t$ は平均 $\mu$, 分散 $\frac{\sigma^2}{2\gamma}$ の正規分布に収束する. ::: ::: {.proof} Step1: @eq-ou-solution の期待値をとる. 右辺第3項は伊藤積分なので平均ゼロであり, $$ \mathbb{E}[x_t] = \mu + (x_0 - \mu)\, e^{-\gamma t}. $$ Step2: @eq-ou-solution より $x_t - \mathbb{E}[x_t] = \sigma \int_0^t e^{-\gamma (t - s)}\, dB_s$ である. 被積分関数 $e^{-\gamma(t-s)}$ は決定論的なので, 伊藤の等長性を用いると $$ \text{Var}(x_t) = \sigma^2\, \mathbb{E}\left[\left(\int_0^t e^{-\gamma (t - s)}\, dB_s\right)^{2}\right] = \sigma^2 \int_0^t e^{-2\gamma (t - s)}\, ds. $$ $\tau = t - s$ と置換すると $\displaystyle \int_0^t e^{-2\gamma (t - s)}\, ds = \int_0^t e^{-2\gamma \tau}\, d\tau = \frac{1 - e^{-2\gamma t}}{2\gamma}$ なので, $$ \text{Var}(x_t) = \frac{\sigma^2}{2\gamma}\left(1 - e^{-2\gamma t}\right). $$ Step 3: 決定論的な被積分関数に対する伊藤積分は正規分布に従うので, @eq-ou-solution より $x_t$ も正規分布である. $t \to \infty$ で平均は $\mu$, 分散は $\dfrac{\sigma^2}{2\gamma}$ に収束するから, $x_t$ は定常分布 $\mathcal{N}\left(\mu,\, \dfrac{\sigma^2}{2\gamma}\right)$ に収束する. ::: ## 最大値原理 ### HJB方程式 (Deterministic) 復習として, ノイズのない決定論的な最適制御問題の Hamilton-Jacobi-Bellman (HJB) 方程式を導きます. ::: {#prp-hjb-det} ## HJB 方程式 (決定論的) 状態 $x_t$ が制御 $u_t$ のもとで $$ \dot{x}_t = g(x_t, u_t) $$ に従い, 意思決定者は割引率 $\rho$ で $$ V(x_0) = \max_{\{u_t\}} \int_0^{\infty} e^{-\rho t} f(x_t, u_t)\, dt $$ を最大化するとする. このとき, 最適価値関数 $V(x)$ は次のHJB方程式を満たす. $$ \rho V(x) = \max_{u}\left\{ f(x, u) + V'(x)\, g(x, u) \right\}. $${#eq-hjb-det} ::: ::: {.proof} 微小区間 $[t, t + dt]$ に最適性原理 (principle of optimality) を適用する. いま $dt$ だけ制御 $u$ を選び, 以降を最適に続けると $$ V(x) = \max_{u}\left\{ f(x, u)\, dt + e^{-\rho\, dt}\, V(x_{t+dt}) \right\}. $$ 決定論的なので $x_{t+dt} = x + g(x, u)\, dt + o(dt)$ であり, テイラー展開より $$ V(x_{t+dt}) = V(x) + V'(x)\, g(x, u)\, dt + o(dt). $$ これに $e^{-\rho\, dt} \approx 1 - \rho\, dt$ を代入し, $dt$ の1次まで残すと $$ V(x) = \max_{u}\left\{ V(x) + \left(f(x, u) + V'(x)\, g(x, u) - \rho V(x)\right) dt \right\} + o(dt). $$ 両辺から $V(x)$ を引いて $dt$ で割り, $dt \to 0$ とすると @eq-hjb-det を得る. ::: ここに確率的なノイズを加えた場合, HJB方程式を導くには伊藤の補題が必要になります. ### 伊藤の補題 ::: {#lem-ito} ## 伊藤の補題 (Itô's lemma) $x_t$ が確率微分方程式 $dx_t = g_t\, dt + \sigma_t\, dB_t$ ($g_t$ はドリフト, $\sigma_t$ は拡散係数で, ともに時刻 $t$ までの情報に依存してよい) に従うとき, 2階連続微分可能な関数 $\phi$ に対して $$ \begin{aligned} d\phi(x_t) &= \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, (dx_t)^2 \\ &= \left[\phi'(x_t)\, g_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\right] dt + \phi'(x_t)\, \sigma_t\, dB_t. \end{aligned} $$ ::: ::: {.proof} $\phi(x_t)$ の微小時間 $[t, t + dt]$ での変化をテイラー展開し, 2次の項まで残すと $$ \phi(x_{t+dt}) - \phi(x_t) = \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, (dx_t)^2 + o\left((dx_t)^2\right). $$ ここに $dx_t = g_t\, dt + \sigma_t\, dB_t$ を代入すると $$ (dx_t)^2 = g_t^2 (dt)^2 + 2 g_t \sigma_t\, dt\, dB_t + \sigma_t^2 (dB_t)^2. $$ @prp-bm-qv の二次変分 $(dB_t)^2 = dt$ を用い, より高次の $(dt)^2$ と $dt\, dB_t$ を $o(dt)$ として落とすと $(dx_t)^2 = \sigma_t^2\, dt$ となる. これを代入すると $$ \begin{aligned} d\phi(x_t) &= \phi'(x_t)\, dx_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\, dt \\ &= \left[\phi'(x_t)\, g_t + \frac{1}{2}\phi''(x_t)\, \sigma_t^2\right] dt + \phi'(x_t)\, \sigma_t\, dB_t. \end{aligned} $$ ::: ### HJB方程式 (Stochastic) 状態が確率微分方程式 $dx_t = g(x_t, u_t)\, dt + \sigma(x_t, u_t)\, dB_t$ に従い, 割引率 $\rho$ で $$ V(x_0) = \max_{\{u_t\}} \mathbb{E}_0 \int_0^{\infty} e^{-\rho t} f(x_t, u_t)\, dt $$ を最大化する確率的な問題を考えます. $V(x)$ を最適価値関数と呼びます. 伊藤の補題を使うと, 決定論的な @eq-hjb-det が次のように拡張されます. ::: {#prp-hjb-stochastic} ## HJB 方程式 (確率的) 状態が $dx_t = g(x_t, u_t)\, dt + \sigma(x_t, u_t)\, dB_t$ に従うとき, 最適価値関数 $V$ は次の Hamilton-Jacobi-Bellman (HJB) 方程式を満たす. $$ \rho V(x) = \max_{u}\left\{ f(x, u) + V'(x)\, g(x, u) + \frac{1}{2}\,\sigma(x, u)^2\, V''(x) \right\}. $$ ::: ::: {.proof} 微小区間 $[t, t + dt]$ に最適性原理 (principle of optimality) を適用する. いま $dt$ だけ制御 $u$ を選び以降を最適に続けると, 状態が確率的なので将来の価値は期待値で評価される. $$ V(x) = \max_{u}\left\{ f(x, u)\, dt + e^{-\rho\, dt}\, \mathbb{E}\left[V(x_{t+dt}) \mid x_t = x\right] \right\}. $$ 伊藤の補題より $dV = \left[V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x)\right] dt + V'(x)\, \sigma(x, u)\, dB_t$ であり, 伊藤積分の項は平均ゼロなので, $$ \mathbb{E}\left[V(x_{t+dt}) \mid x_t = x\right] = V(x) + \left[V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x)\right] dt. $$ これと $e^{-\rho\, dt} \approx 1 - \rho\, dt$ を代入し, $dt$ の1次まで残すと $$ V(x) = \max_{u}\left\{ V(x) + \left(f(x, u) + V'(x)\, g(x, u) + \tfrac{1}{2}\sigma(x, u)^2\, V''(x) - \rho V(x)\right) dt \right\} + o(dt). $$ 両辺から $V(x)$ を引いて $dt$ で割り, $dt \to 0$ とすると命題の HJB 方程式を得る. ::: ## 連続時間の Aiyagari モデルここからは, 上で整えた道具を使って, 連続時間の異質エージェントモデルを実際に解きます. 題材は Aiyagari モデル (@sec-dp-hetero) の連続時間版です. 生産性の対数を OU 過程でモデル化するので, HJB 方程式には伊藤の補題が生む2階微分の項がそのまま現れます. 数値解法は @achdou2022 の有限差分法です. 離散時間版と比べながら, 連続時間の数値計算がなぜ速く, なぜ頑健なのかを見ていきます. ### 設定家計は割引率 $\rho$ で期待生涯効用 $\mathbb{E}_0 \int_0^\infty e^{-\rho t} u\left(c_t\right) dt$ を最大化します. 資産 $a$ は資本で, 予算制約と借入制約は $$ \dot{a}_t = w z_t + r a_t - c_t, \qquad a_t \geq 0 $$ です. 労働生産性 $z_t = e^{x_t}$ の対数は OU 過程 $$ dx_t = -\theta x_t\, dt + \sigma\, dB_t $$ に従います (@prp-ou-solution の記法で $\gamma = \theta$, $\mu = 0$ とし, 平均回帰の速さを CRRA 係数と区別して $\theta$ と書きます). OU 過程は AR(1) の連続時間版なので, 離散時間の章との対応は明快です. 1年間隔で観測すると自己相関 $e^{-\theta}$, 定常分散 $\sigma^2 / \left(2\theta\right)$ の AR(1) になるので (@prp-ou-moment), 動的計画法の章のカリブレーション $\left(\rho_z, \sigma_\varepsilon\right) = \left(0.9, 0.2\right)$ をそのまま $$ \theta = -\log 0.9, \qquad \frac{\sigma^2}{2\theta} = \frac{\sigma_\varepsilon^2}{1 - \rho_z^2} $$ で写し取れます. 割引率も同様に $\rho = -\log \beta = -\log 0.96$ とします. 企業は Cobb-Douglas 生産関数 $K^\alpha L^{1-\alpha}$ を持ち, $k = K / L$ として $r = \alpha k^{\alpha - 1} - \delta$, $w = \left(1 - \alpha\right) k^{\alpha}$ が成り立ちます. 労働供給は非弾力的で, $\mathbb{E}\left[z\right] = 1$ に基準化すれば $L = 1$ です. ::: {#def-hact-equilibrium} ## 定常均衡の方程式系定常均衡は, 価値関数 $v\left(a, x\right)$, 密度 $g\left(a, x\right)$, 価格 $\left(r, w\right)$ の組であって次を満たすものである. **HJB 方程式**: $$ \rho v = \max_{c}\ u\left(c\right) + \partial_a v \left(w e^x + ra - c\right) - \theta x\, \partial_x v + \frac{\sigma^2}{2} \partial_{xx} v. $${#eq-hjb} **Kolmogorov Forward (KF) 方程式** (定常): $$ 0 = -\partial_a \left[s\left(a, x\right) g\right] + \partial_x \left[\theta x\, g\right] + \frac{\sigma^2}{2} \partial_{xx} g, $${#eq-kf} ここで $s\left(a, x\right) = w e^x + ra - c\left(a, x\right)$ は貯蓄関数, $c\left(a, x\right) = \left(u'\right)^{-1}\left(\partial_a v\right)$. **資本市場の清算**: $$ \iint a\, g\left(a, x\right)\, da\, dx = K, \qquad r = \alpha \left(K / L\right)^{\alpha - 1} - \delta. $${#eq-clearing} ::: HJB 方程式が価値関数を「後ろ向き」に決め, その政策関数のもとで KF 方程式が分布を「前向き」に決め, 分布が価格 $r$ を決めて HJB に戻る, というループ構造です.^[この「後ろ向きの HJB + 前向きの KF」という連立系は, 平均場ゲーム (Mean Field Game) の定常系そのものです. @achdou2022 の枠組みは, Lasry と Lions の MFG 理論を異質エージェント・マクロに持ち込んだものといえます.] アルゴリズムの全体像を先に図で示しておきます (@fig-hact-loop). 以下の数値解法の節は, この図の各ブロックを上から順に実装していきます. ![Equilibrium loop of the finite difference method](/static/cetz/hact_loop.svg){#fig-hact-loop width=55%} **借入制約の入り方**が離散時間と決定的に違います. 連続時間では制約は内点で決して binding にならないため, 一階条件 $u'\left(c\left(a, x\right)\right) = \partial_a v\left(a, x\right)$ が $a > 0$ の**すべての点で等号で**成り立ちます. 制約は HJB 方程式そのものではなく, 境界条件 $$ \partial_a v\left(0, x\right) \geq u'\left(w e^x\right) $${#eq-state-constraint} として入ります. これは一階条件を通じて $s\left(0, x\right) \geq 0$ を含意するので, 制約が破られることはありません. KF 方程式には境界条件すら不要です. 最適政策が制約を尊重している限り, 分布は自動的に $\left[0, \infty\right)$ に収まるからです. $x$ 方向はグリッドを定常標準偏差の $\pm 3$ 倍で打ち切り, 端では反射壁 (reflecting barrier) を仮定します. OU 過程のドリフトは端では常に内側を向くので, この打ち切りの影響はわずかです. ## 有限差分法 ### Upwind スキーム資産グリッド $a_1 = 0 < a_2 < \cdots < a_{N_a}$ (等間隔, 幅 $\Delta a$) を取り, $x$ を固定して資産方向の微分を前方差分と後方差分で近似します. $$ \partial_a v_{i, F} = \frac{v_{i+1} - v_i}{\Delta a}, \qquad \partial_a v_{i, B} = \frac{v_i - v_{i-1}}{\Delta a}. $$ どちらを使うべきでしょうか. 答えは「**状態が動く向きに合わせる**」です. 貯蓄 (ドリフト) が正なら前方差分, 負なら後方差分を使う. これが upwind スキームで, $s_i^{+} = \max\left\{s_i, 0\right\}$, $s_i^{-} = \min\left\{s_i, 0\right\}$ として HJB 方程式の資産方向を $$ \frac{v_{i+1} - v_i}{\Delta a}\, s_i^{+} + \frac{v_i - v_{i-1}}{\Delta a}\, s_i^{-} $${#eq-upwind} と離散化します (@fig-upwind). ![Upwind scheme: the difference follows the drift](/static/cetz/upwind_scheme.svg){#fig-upwind width=75%} 有限差分スキームが真の解 (粘性解) に収束するための Barles--Souganidis の十分条件は単調性・整合性・安定性の3つで, このうち非自明なのが単調性です. Upwind は, まさに単調性を保証するための差分の選び方になっています [@achdou2022, sec. 5.3]. 実装上は, まず前方差分で仮の貯蓄 $s_{i,F} = w e^{x} + r a_i - \left(u'\right)^{-1}\left(\partial_a v_{i,F}\right)$ を, 後方差分で $s_{i,B}$ を計算し, - $s_{i,F} > 0$ なら前方差分を採用する - $s_{i,B} < 0$ なら後方差分を採用する - どちらでもない場合 ($s_{i,F} \leq 0 \leq s_{i,B}$) はその点を**定常点**とみなして $s_i = 0$, $c_i = w e^{x} + r a_i$ とするという規則にします.^[なお @achdou2022 の本文はスキームの詳細を数値付録に委ねているため, この規則と後述の陰的スキームは, 公開されている数値付録と実装コードに従っています.] 境界条件の扱いは驚くほど簡単です. 下端 $a_1 = 0$ では**後方差分だけ**に状態制約の境界条件 @eq-state-constraint を等号で課します: $\partial_a v_{1, B} = u'\left(w e^{x}\right)$. 前方差分は普通に計算し, どちらを使うかは upwind に選ばせます. 貯蓄が正なら境界条件は使われず, 負になりそうなときだけ発動して $s = 0$ で止める, という仕掛けです. 上端では貯蓄が負になるので前方差分は使われず, 境界条件は実質不要です. ### 生産性方向の離散化 OU 過程の項 $-\theta x\, \partial_x v + \frac{\sigma^2}{2} \partial_{xx} v$ も同じ発想で離散化します. ドリフト $-\theta x$ は upwind (符号に応じて片側差分), 拡散項は中心差分です. グリッド $x_1 < \cdots < x_{N_x}$ (幅 $\Delta x$) 上の点 $j$ について, 隣へ移る「強度」をまとめると $$ \lambda_j^{\downarrow} = \frac{\sigma^2}{2 \Delta x^2} + \frac{\max\left\{\theta x_j, 0\right\}}{\Delta x}, \qquad \lambda_j^{\uparrow} = \frac{\sigma^2}{2 \Delta x^2} + \frac{\max\left\{-\theta x_j, 0\right\}}{\Delta x} $${#eq-ou-intensity} となります ($x_j > 0$ なら平均回帰のドリフトは下向きに働く, という形です). 端のグリッドでは外向きの強度を落とします. これが反射壁にあたります. できあがる三重対角行列 $B_x$ は各行の和がゼロの**強度行列** (連続時間マルコフ連鎖の生成作用素) です. つまりこの離散化は, AR(1) を有限状態マルコフ連鎖で近似した Tauchen 法・Rouwenhorst 法 (@sec-discretization) の連続時間版です.^[upwind のドリフト近似は分散をわずかに水増しします. Tauchen 法で見た「升目への丸めが余分な分散を生む」現象と同種のもので, $\Delta x \to 0$ で消えます. 下の実装では離散化された定常分布の標準偏差を理論値と比べて確認します.] ### 遷移行列 $A$ 以上を合わせると, 離散化された HJB 方程式は $$ \rho v = u + A\left(v\right) v $${#eq-hjb-discrete} という形になります. $A$ は $N_a N_x \times N_a N_x$ の疎行列で, 資産方向の upwind 部分 (三重対角) と生産性方向の $B_x$ (帯2本) の和です. 重要な性質が2つあります. - **各行の和はゼロ**: $A$ は全体として連続時間マルコフ連鎖の強度行列である - 有限差分法は, 富と生産性の連続な動きを「隣のグリッド点へのポアソンジャンプ」で近似していると解釈できる. 離散時間で遷移確率行列を作ったことの, 連続時間版にあたる @fig-fd-stencil がこのジャンプの解釈を図解しています. 状態 $\left(a_i, x_j\right)$ から確率的に移れる先は, **上下左右の隣接グリッドの4点だけ**です. 資産方向は upwind が選んだ片側 (実際に正の強度を持つのは矢印2本のうち1本), 生産性方向は @eq-ou-intensity の2本です. ![Poisson-jump interpretation of the generator](/static/cetz/fd_stencil.svg){#fig-fd-stencil width=60%} この「隣にしか跳べない」構造が, そのまま $A$ の疎パターンになります. @fig-hact-sparsity は小さいグリッド ($N_a = 12$, $N_x = 4$) で $A$ の非ゼロ要素の位置を描いたものです. 資産方向のジャンプは各ブロック内の**三重対角** (対角とその上下1本ずつ), 生産性方向のジャンプはブロックを跨ぐ $\pm N_a$ の位置の**帯2本**として現れ, $A$ の各行の非ゼロ要素は高々5個です. ```{julia} #| label: fig-hact-sparsity #| fig-cap: "**Sparsity pattern of $A$** ($N_a = 12$, $N_x = 4$ for illustration)." #| code-fold: true Nad, Nxd = 12, 4 Ta = spdiagm(-1 => ones(Nad - 1), 0 => fill(2.0, Nad), 1 => ones(Nad - 1)) Bd = spdiagm(-1 => ones(Nxd - 1), 0 => fill(2.0, Nxd), 1 => ones(Nxd - 1)) Ad = blockdiag((Ta for _ in 1:Nxd)...) + kron(Bd, sparse(I, Nad, Nad)) nzrows, nzcols, _ = findnz(Ad) scatter(nzcols, nzrows, yflip=true, aspect_ratio=1, markersize=2.5, color=:black, markerstrokewidth=0, legend=false, xlabel="Column", ylabel="Row", xlims=(0, Nad * Nxd + 1), ylims=(0, Nad * Nxd + 1), size=(380, 380)) vline!(collect(Nad:Nad:(Nad*(Nxd-1))) .+ 0.5, color=:gray, ls=:dash, alpha=0.6) hline!(collect(Nad:Nad:(Nad*(Nxd-1))) .+ 0.5, color=:gray, ls=:dash, alpha=0.6) ``` 行列のサイズは $N_a N_x$ の2乗で増えますが, 非ゼロ要素は $5 N_a N_x$ 程度にしか増えません. 疎行列ソルバはゼロでない要素だけを保存・計算するので, 後で出てくる連立一次方程式 (陰的スキームと KF 方程式) は巨大なグリッドでも一瞬で解けます. これが有限差分法の速さの核心です. ### 陰的スキーム @eq-hjb-discrete は $A$ が $v$ に依存する非線形方程式なので, 反復で解きます. 擬似的な時間ステップ $\Delta$ を入れた**陰的スキーム** (implicit scheme) を使います. $$ \frac{v^{n+1} - v^n}{\Delta} + \rho v^{n+1} = u\left(c^n\right) + A\left(v^n\right) v^{n+1} \quad \Longleftrightarrow \quad \left[\left(\frac{1}{\Delta} + \rho\right) I - A^n\right] v^{n+1} = u^n + \frac{v^n}{\Delta}. $${#eq-implicit} 各反復は疎な連立一次方程式を1回解くだけです. 陰的スキームは無条件に安定なので, $\Delta$ をいくらでも大きく取れます ($\Delta = 1000$ など). その結果, 反復は**数十回**で収束します. 素朴な価値関数反復が何百回もかかるのと比べると大きな違いです (離散時間との定量的な比較は実装のあとで行います). ### KF 方程式はタダで解ける分布はどう求めるのでしょうか. ここが有限差分法の最も美しい部分です. KF 作用素は HJB に現れる生成作用素の**随伴** (adjoint) なので, 離散化すると @eq-kf は $$ A^{\top} g = 0 $${#eq-kf-discrete} になります [@achdou2022, sec. 5.4]. HJB を解いた時点で $A$ は手元にあるので, KF 方程式は**転置行列の連立一次方程式を1回解くだけ**です. 分布のシミュレーションも, 遷移行列の別途構築も要りません. @eq-kf-discrete は特異 (解が定数倍で不定) なので, 1つの成分を固定して解いてから, 密度の積分が1になるように正規化します. 借入制約に張り付く家計の質量は, 下端のグリッド点に立つ大きな値として現れます. ### 均衡利子率を与えると企業の一階条件から $k\left(r\right) = \left(\alpha / \left(r + \delta\right)\right)^{1 / \left(1 - \alpha\right)}$ と $w\left(r\right)$ が決まるので, 資本の超過供給関数 $$ S\left(r\right) = \iint a\, g\, da\, dx - k\left(r\right) L $$ の二分法で均衡利子率を求めます. $r$ を与えて HJB → KF → $S\left(r\right)$ と計算し, 正なら $r$ を下げ, 負なら上げるだけです. まとめると, この手法が速く頑健な理由は4つです [@achdou2022, sec. 5.1]. 1. 借入制約は境界条件に封じ込められ, 内点の一階条件は常に等号で成り立つ (occasionally binding な場合分けが消える) 1. 一階条件が静学的で, $c = \left(\partial_a v\right)^{-1/\gamma}$ と手で解ける. 求根も期待値の計算も不要 1. 連続時間の過程はグリッド上を1点ずつしか動かないので, 行列が疎になる 1. KF 方程式が HJB の転置としてタダで手に入る ## Julia 実装カリブレーションは動的計画法の章の Aiyagari モデルを OU 過程に写したものです. ```{julia} #| label: hact-params #| output: false const γ, ρ = 2.0, -log(0.96) # CRRA, discount rate (= -log β) const α, δ = 0.36, 0.08 # production const θ = -log(0.9) # OU mean reversion: e^{-θ} = 0.9 const σx = 0.2 / sqrt(1 - 0.9^2) # stationary sd of x (matches the AR(1)) const σ = σx * sqrt(2θ) # OU diffusion coefficient const amin, amax, Na = 0.0, 50.0, 2000 const a = range(amin, amax, length=Na) const da = step(a) const Nx = 31 const x = range(-3σx, 3σx, length=Nx) const dx = step(x) u(c) = c^(1 - γ) / (1 - γ) u′(c) = c^(-γ) u′inv(v) = max(v, 1e-10)^(-1 / γ) # guard against negative slopes early on ``` まず生産性方向の生成作用素 $B_x$ を @eq-ou-intensity のとおり組み立てます. その定常分布を使って $\mathbb{E}\left[z\right] = 1$ に基準化し, 離散化が OU 過程の定常標準偏差を再現しているかも確認しておきます. ```{julia} #| label: hact-ou function ou_generator() B = spzeros(Nx, Nx) for j in 1:Nx lo = σ^2 / (2dx^2) + max(θ * x[j], 0) / dx # intensity of moving down hi = σ^2 / (2dx^2) + max(-θ * x[j], 0) / dx # intensity of moving up if j > 1 ; B[j, j-1] = lo; B[j, j] -= lo; end # boundary: drop the outward jump if j < Nx ; B[j, j+1] = hi; B[j, j] -= hi; end end return B end const Bx = ou_generator() function stat_x(B) # stationary distribution over x BT = copy(sparse(transpose(B))) b = zeros(Nx); BT[1, :] .= 0; BT[1, 1] = 1.0; b[1] = 0.1 ψ = BT \ b ψ ./ sum(ψ) end const ψ = stat_x(Bx) const z = exp.(x) ./ dot(ψ, exp.(x)) # normalize E[z] = 1 const Bswitch = kron(Bx, sparse(I, Na, Na)) # x-part on the full (a, x) grid (sd_discretized=round(sqrt(dot(ψ, x .^ 2) - dot(ψ, x)^2), digits=4), sd_theory=round(σx, digits=4)) ``` HJB 方程式のソルバは, upwind による政策の計算, 遷移行列 $A$ の構築, 陰的更新 @eq-implicit の3ステップの反復です. ```{julia} #| label: hact-hjb #| output: false function solve_hjb(r, w; Δ=1000.0, maxit=400, tol=1e-8) inc = w .* z' .+ r .* a # inc[i, j] = w z_j + r a_i v = u.(inc) ./ ρ dVf = similar(v); dVb = similar(v) A = spzeros(Na * Nx, Na * Nx); c = similar(v) for it in 1:maxit dVf[1:(end-1), :] = diff(v, dims=1) ./ da dVf[end, :] = u′.(inc[end, :]) # unused: saving is negative up top dVb[2:end, :] = diff(v, dims=1) ./ da dVb[1, :] = u′.(inc[1, :]) # state-constraint boundary condition cf = u′inv.(dVf); sf = inc .- cf # saving under the forward difference cb = u′inv.(dVb); sb = inc .- cb # saving under the backward difference If = sf .> 0 Ib = (sb .< 0) .& .!If # where neither drift applies the household is at rest, so s = 0 c = cf .* If .+ cb .* Ib .+ inc .* (.!If .& .!Ib) X = -min.(sb, 0) ./ da # intensity of moving one grid point down Z = max.(sf, 0) ./ da # intensity of moving one grid point up A = blockdiag( (spdiagm(-1 => X[2:end, j], 0 => -X[:, j] .- Z[:, j], 1 => Z[1:(end-1), j]) for j in 1:Nx)...) + Bswitch B = (1 / Δ + ρ) * sparse(I, Na * Nx, Na * Nx) - A vnew = reshape(B \ (vec(u.(c)) .+ vec(v) ./ Δ), Na, Nx) maximum(abs.(vnew .- v)) < tol && return (v=vnew, c=c, s=inc .- c, A=A, iter=it) v = vnew end error("HJB did not converge") end ``` KF 方程式は転置行列の連立一次方程式です. 特異性は1成分を固定して外し, 最後に密度として正規化します. ```{julia} #| label: hact-kf #| output: false function stationary_dist(A) AT = copy(sparse(transpose(A))) b = zeros(Na * Nx) AT[1, :] .= 0.0; AT[1, 1] = 1.0; b[1] = 0.1 # fix one entry of the singular system g = AT \ b g ./= sum(g) * da * dx # normalize to a density return reshape(g, Na, Nx) end ``` 最後に, 資本の超過供給に対する二分法で均衡利子率を求めます. ```{julia} #| label: hact-equilibrium prices(r) = (k=(α / (r + δ))^(1 / (1 - α)); (k=k, w=(1 - α) * k^α)) function excess_supply(r) (; k, w) = prices(r) sol = solve_hjb(r, w) g = stationary_dist(sol.A) return sum(g .* a) * da * dx - k, sol, g end function equilibrium(rlo, rhi; n=20) for _ in 1:n rm = (rlo + rhi) / 2 first(excess_supply(rm)) > 0 ? (rhi = rm) : (rlo = rm) end return (rlo + rhi) / 2 end r★ = equilibrium(0.005, ρ - 1e-4) ex, sol, g = excess_supply(r★) (r=round(r★, digits=4), K=round(prices(r★).k, digits=2), hjb_iter=sol.iter) ``` 均衡利子率は割引率 $\rho = -\log 0.96 \approx 0.041$ を明確に下回ります. 予備的貯蓄による過剰な資本蓄積という Aiyagari モデルの中心的な結果が, 連続時間でもそのまま現れています. 離散時間の章で得た $r^* \approx 0.021$ とはカリブレーションの写像や離散化の違いから小さなずれがありますが, 定性的な結論は同じです. そして, グリッドが $N_a N_x = 62{,}000$ 点あるにもかかわらず, HJB 方程式の反復が数十回, 均衡全体も二分法 20 回で数分以内に収まっていることに注目してください. 速さの源泉は @fig-hact-sparsity で見た疎性です. 実際の $A$ は $62{,}000 \times 62{,}000$ ですが, 非ゼロ要素は全体の 0.01% にも遠く及びません. ```{julia} #| label: hact-nnz (nonzeros=nnz(sol.A), share_percent=round(100 * nnz(sol.A) / (Na * Nx)^2, sigdigits=2)) ``` ### 政策関数 ```{julia} #| label: fig-hact-policy #| fig-cap: "**Saving policy functions.** Low-productivity households at the constraint have $s = 0$." #| code-fold: true jmid = (Nx + 1) ÷ 2 idx = a .<= 10.0 plot(a[idx], sol.s[idx, 1], label=L"x = -3\sigma_x", lw=2) plot!(a[idx], sol.s[idx, jmid], label=L"x = 0", lw=2, ls=:dash) plot!(a[idx], sol.s[idx, Nx], label=L"x = 3\sigma_x", lw=2, ls=:dashdot) hline!([0], color=:gray, ls=:dot, label=false) plot!(xlabel=L"a", ylabel="Saving", legend=:topright) ``` 生産性の低い家計は資産を取り崩しながら借入制約に向かって滑り落ち, 制約点でちょうど $s = 0$ になります. 内点で制約が binding にならない, という理論の構造がそのまま見えています. 生産性の高い家計は広い範囲で貯蓄を続けますが, 十分に資産を持つと ($\rho > r^*$ なので) やはり取り崩しに転じます. ### 定常分布 ```{julia} #| label: fig-hact-dist #| fig-cap: "**Stationary asset distribution.** The point mass at the constraint is omitted." #| code-fold: true ga = vec(sum(g, dims=2)) .* dx # marginal density of assets plot(a[2:end][a[2:end] .<= 30], ga[2:end][a[2:end] .<= 30], label=false, lw=2, xlabel=L"a", ylabel="Density") ``` ```{julia} #| label: hact-mass println("mass at the borrowing constraint: ", round(sum(g[1, :]) * da * dx, sigdigits=3)) ``` 離散時間版と同じく, 借入制約の近くに密度が集中し, 制約点そのものには質量が立ちます (図では質量を除いた密度部分を描いています). KF 方程式に一切の境界条件を与えていないのに, 質量まで含めて正しい分布が「転置を解くだけ」で出てくるのが, この手法の強みです. ### 離散時間モデルとの比較連続時間版が本当に離散時間の Aiyagari モデルを再現しているのか, 同じカリブレーションの離散時間モデルを解いて確かめます. 実装は動的計画法の章と同じ, 単調性を利用したグリッドサーチの VFI に Howard の改善を加えたものです. 所得過程の離散化には Rouwenhorst 法を使います. Tauchen 法は少ない状態数で分散を過大評価するため (@sec-discretization), 比較がその分だけ歪んでしまうからです. Rouwenhorst 法は AR(1) の2次モーメントを厳密に再現するので, 離散時間の所得過程は連続時間の OU 過程と (2次モーメントの意味で) 完全に揃います. ```{julia} #| label: dt-solver #| code-fold: true #| output: false const β = exp(-ρ) # discrete-time discount factor function dt_setup(; na, nz, uniform=false) mc = rouwenhorst(nz, 0.9, 0.2) Π = mc.p zdt = exp.(collect(mc.state_values)) zdt ./= dot(vcat(stationary_distributions(mc)...), zdt) # E[z] = 1 ag = uniform ? collect(range(0, amax, length=na)) : collect(range(0, 1, length=na)) .^ 2 .* amax (; Π, zdt, ag, na, nz) end function dt_vfi(r, w, m; howard=50, tol=1e-6, maxit=1000) (; Π, zdt, ag, na, nz) = m V = [u(w * zdt[iz] + r * ag[ia] + 1e-10) / (1 - β) for ia in 1:na, iz in 1:nz] pol = ones(Int, na, nz) iters = 0 for it in 1:maxit iters = it EV = β .* (V * Π') Vnew = similar(V) for iz in 1:nz lo = 1 for ia in 1:na res = w * zdt[iz] + (1 + r) * ag[ia] best = -Inf besti = lo for ja in lo:na # monotone grid search c = res - ag[ja] c <= 0 && break val = u(c) + EV[ja, iz] if val > best best = val besti = ja else break # concavity: passed the peak end end Vnew[ia, iz] = best pol[ia, iz] = besti lo = besti end end for _ in 1:howard # Howard improvement EVh = β .* (Vnew * Π') for iz in 1:nz, ia in 1:na c = w * zdt[iz] + (1 + r) * ag[ia] - ag[pol[ia, iz]] Vnew[ia, iz] = u(c) + EVh[pol[ia, iz], iz] end end if maximum(abs.(Vnew .- V)) < tol V = Vnew break end V = Vnew end return V, pol, iters end function dt_dist(pol, m; tol=1e-10, maxit=100_000) (; Π, na, nz) = m n = na * nz Is = Int[]; Js = Int[]; Vs = Float64[] for iz in 1:nz, ia in 1:na row = ia + (iz - 1) * na for jz in 1:nz push!(Is, row) push!(Js, pol[ia, iz] + (jz - 1) * na) push!(Vs, Π[iz, jz]) end end Tt = sparse(Js, Is, Vs, n, n) # transposed transition matrix λ = fill(1 / n, n) iters = 0 for it in 1:maxit iters = it λn = Tt * λ if maximum(abs.(λn .- λ)) < tol λ = λn break end λ = λn end return reshape(λ, na, nz), iters end dt_assets(λ, m) = sum(λ .* m.ag) ``` まず, 資本の需要曲線と供給曲線を両方のモデルで重ねてみます (@fig-hact-vs-dt). ```{julia} #| label: fig-hact-vs-dt #| fig-cap: "**Capital demand and supply.** Continuous vs. discrete time, same calibration." #| code-fold: true mdt = dt_setup(na=200, nz=7) function dt_supply(r) w = prices(r).w _, pol, _ = dt_vfi(r, w, mdt) dt_assets(first(dt_dist(pol, mdt)), mdt) end ct_supply(r) = first(excess_supply(r)) + prices(r).k rvals = range(0.005, 0.033, length=8) A_dt = dt_supply.(rvals) A_ct = ct_supply.(rvals) rfine = range(0.005, 0.033, length=100) plot(map(r -> prices(r).k, rfine), rfine, label=L"Demand $k(r)$", color=:gray, ls=:dot, lw=2) plot!(A_ct, rvals, label="Supply (continuous time)", lw=2, marker=:circle) plot!(A_dt, rvals, label="Supply (discrete time)", lw=2, ls=:dash, marker=:square) plot!(xlabel="Assets / capital", ylabel=L"r", legend=:bottomright) ``` 供給曲線はほぼ重なり, 需要曲線との交点 (均衡利子率) も互いに近い位置にあります. ```{julia} #| label: dt-ct-rstar function dt_equilibrium(rlo, rhi; n=20) for _ in 1:n rm = (rlo + rhi) / 2 dt_supply(rm) > prices(rm).k ? (rhi = rm) : (rlo = rm) end return (rlo + rhi) / 2 end (r_discrete=round(dt_equilibrium(0.005, 0.035), digits=4), r_continuous=round(r★, digits=4)) ``` 残る小さな差は, 意思決定のタイミング (離散時間では1年分まとめて意思決定する) と離散化そのものの違いによるものです. 連続時間版は離散時間の Aiyagari モデルの答えを再現しつつ, 滑らかな政策関数と分布を返している, と言ってよいでしょう. 次に計算時間です. 同じ均衡価格のもとで, 家計問題と定常分布の計算をそれぞれ計測します (収束判定はどちらも $\max \left|\Delta v\right| < 10^{-6}$). ```{julia} #| label: hact-bench #| output: false mdt31 = dt_setup(na=Na, nz=Nx, uniform=true) # same grid as continuous time w★ = prices(r★).w dt_vfi(r★, w★, mdt31) # warmup t_vfi7 = @elapsed ((_, pol7, it7) = dt_vfi(r★, w★, mdt)) t_dist7 = @elapsed ((λ7, itd7) = dt_dist(pol7, mdt)) t_vfi31 = @elapsed ((_, pol31, it31) = dt_vfi(r★, w★, mdt31)) t_dist31 = @elapsed ((λ31, itd31) = dt_dist(pol31, mdt31)) t_hjb = @elapsed (solb = solve_hjb(r★, w★; tol=1e-6)) t_kf = @elapsed stationary_dist(solb.A) cell(t, n, unit) = string(round(t, digits=2), " 秒 / ", n, " ", unit) ``` ```{julia} #| label: tbl-hact-bench #| output: asis #| tbl-cap: "家計問題と定常分布の計算時間" benchrows = [ "離散時間 VFI + Howard" raw"$200 \times 7$" cell(t_vfi7, it7, "反復") cell(t_dist7, itd7, "反復"); "離散時間 VFI + Howard" raw"$2{,}000 \times 31$" cell(t_vfi31, it31, "反復") cell(t_dist31, itd31, "反復"); "連続時間有限差分" raw"$2{,}000 \times 31$" cell(t_hjb, solb.iter, "反復") cell(t_kf, 1, "回の求解") ] pretty_table(benchrows; column_labels=["手法", raw"状態数 $N_a \times N_z$", "家計問題", "定常分布"], backend=:markdown, allow_markdown_in_cells=true) ``` 結果は一見すると意外かもしれません. **家計問題単体では, チューニングされた離散時間の VFI も十分に速い**のです. これは選択変数が1次元で, 政策の単調性と価値関数の凹性を使うとグリッドサーチがほぼ線形時間で済むという, この問題に固有の事情によります. 連続時間の優位がはっきり出るのは**定常分布**です. 離散時間では遷移行列の反復 (べき乗法) を数百回まわす必要があるのに対し, 連続時間では HJB の副産物 $A$ の転置を1回解くだけで済みます. より本質的な違いは, 表の数字の外にあります. - 離散時間の速さは「1次元の選択を単調なグリッドサーチで解ける」ことに依存する. 選択変数が増える (ポートフォリオ選択, 耐久財など) と最大化のコストが急速に増すのに対し, 連続時間の一階条件は静学的なまま変わらない - 連続時間の政策関数はグリッドへの丸めがなく滑らかで, 借入制約もちょうど1点で処理される. 推定 (@sec-md) のように解を何千回も呼ぶ用途では, この滑らかさが目的関数の滑らかさに直結する ## まとめと拡張この章のアルゴリズムをまとめると次のとおりです. 1. $r$ を所与として企業の一階条件から $\left(k, w\right)$ を求め, upwind + 陰的スキームで HJB 方程式 @eq-hjb-discrete を解く (反復数十回) 1. 遷移行列の転置から KF 方程式 @eq-kf-discrete を1発で解き, 定常分布を得る 1. 資本市場の超過供給を見て $r$ を二分法で更新する同じ骨格が, ほとんど変更なしに広いクラスのモデルへ拡張できます [@achdou2022, sec. 7]. - **より単純な所得過程**: 所得を2状態のポアソン過程に落とすと2階微分の項が消えて HJB は常微分方程式になり, @achdou2022 の本文で解析される Huggett モデルになる. 借入制約近傍の貯蓄関数や資産分布の質量について, 鋭い理論的特徴づけが得られる - **移行動学**: 時間依存の HJB を終端条件から**後ろ向き**に, KF を初期分布から**前向き**に解き, 価格の経路を収束するまで更新する. MIT ショックの分析がそのまま可能 - **非凸性・停止問題**: 住宅のような不可分性や, デフォルトのような停止問題も, 粘性解の枠組みのもとで同じ離散化で扱える離散時間の動的計画法 (@sec-dynamic-programming) と対をなす計算技術として, この有限差分法は現代の異質エージェント・マクロの標準的な道具になっています. 実装のバリエーションは @achdou2022 の公開コードが網羅しており, この章のコードはその最小構成にあたります.

B.1 確率過程

B.1.1 ブラウン運動

B.1.2 伊藤積分

B.1.3 Orstein-Uhlenbeck 過程

B.2 最大値原理

B.2.1 HJB方程式 (Deterministic)

B.2.2 伊藤の補題

B.2.3 HJB方程式 (Stochastic)

B.3 連続時間の Aiyagari モデル

B.3.1 設定

B.4 有限差分法

B.4.1 Upwind スキーム

B.4.2 生産性方向の離散化

B.4.3 遷移行列 \(A\)

B.4.4 陰的スキーム

B.4.5 KF 方程式はタダで解ける

B.4.6 均衡

B.5 Julia 実装

B.5.1 政策関数

B.5.2 定常分布

B.5.3 離散時間モデルとの比較

B.6 まとめと拡張