\[ \newcommand{\R}{\mathbb{R}} \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\G}{\mathcal{G}} \newcommand{\abs}[1]{\lvert#1\rvert} \newcommand{\Abs}[1]{\left|#1\right|} \newcommand{\norm}[1]{\|#1\|} \newcommand{\Norm}[1]{\left\|#1\right\|} \newcommand{\SquareBrac}[1]{\left[#1\right]} \newcommand{\Brace}[1]{\left\{#1\right\}} \newcommand{\Paren}[1]{\left(#1\right)} \newcommand{\brac}[1]{\langle#1\rangle} \newcommand{\Brac}[1]{\left\langle#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\ov}[1]{\overline{#1}} \newcommand{\un}[1]{\underline{#1}} \newcommand{\wt}[1]{\widetilde{#1}} \newcommand{\wh}[1]{\widehat{#1}} \newcommand{\pp}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\ppp}[3]{\frac{\partial #1}{\partial #2\partial #3}} \newcommand{\dd}[2]{\frac{d #1}{d #2}} \newcommand{\floor}[1]{\lfloor#1\rfloor} \newcommand{\Floor}[1]{\left\lfloor#1\right\rfloor} \newcommand{\ceil}[1]{\lceil#1\rceil} \newcommand{\Ceil}[1]{\left\lceil#1\rceil\right} \newcommand{\iso}{\overset{\sim}{\to}} \newcommand{\eqDescribe}[2]{\underbrace{#1}_{\text{#2}}} \newcommand{\bm}[1]{\boldsymbol{#1}} \newcommand{\argmax}{\operatorname*{argmax}} \newcommand{\argmin}{\operatorname*{argmin}} \]
今年の6月に『定量的マクロ経済学と数値計算』という本が出版された (北尾早霧, 砂川武貴, and 山田知明 2024)．
近年マクロ経済学のモデルが複雑になるにつれて，コンピュータを使って近似的にモデルの性質を理解するアプローチが広まっている．
しかし，日本のマクロ経済学教育の中で，それを可能にする数値計算 (numerical computation) という分析手法を学ぶ機会は乏しい．
この状況を憂いた著者陣が，『経済セミナー』という雑誌に投稿した連載を書籍化したのが本書である．
本書が扱う内容は，現代的なマクロ経済学のベースラインとなりつつあるトピックであり，教育が手薄になるのは好ましくない．
本書が扱う内容は現在のマクロ経済学のメインストリームの１つといっても過言ではない。＜中略＞本書を通じて現代のマクロ経済分析に必要な新たな道具を身に付けるサポートをすることが、我々の目的である。(p.8)
（気持ちネットワークを扱うモデルをメインにやっているとはいえ，）マクロ経済学徒の一人としては現代的な分析手法を学ぶ必要があると感じたため，勉強ノートを作ろうと思い立った次第である．
このウェブサイトのフレームワークとなっている Quarto は，このようなハンズオンで学ぶスタイルの勉強ノート作成に適している． Quarto もまだまだ勉強中の身ゆえ，このブログは Quarto の勉強も兼ねている．
まず，シンプルな２期間モデルを用いて数値計算手法の基本的な使い方のイメージをつかむ．
現在のマクロ経済学でベンチマークとなっている新古典派成長モデルは，この２期間モデルを多期間に拡張していき，その極限をとって無限期間にしたものである．
したがって，まずはこのシンプルな２期間モデルで数値計算の基礎をおさえておくことが重要である．
ある経済主体の人生全体での消費・貯蓄行動をモデル化しよう．
経済主体は若年期に働いて所得 \(w\) を獲得し，その所得を若年期の消費 \(c_1\) に充てるか，それとも老後のための貯蓄 \(a\) に残すかを決める問題に直面している．
若年期の予算制約は，
\[ c_1 + a = w \tag{1}\]
であり，老年期の予算制約は若年期の貯蓄に金利 \(r\) が付き，また遺産動機はないとすると，
\[ c_2 = a(1 + r) \tag{2}\]
で与えられる．\(\beta > 0\) を割引因子 (discount factor) とすると，経済主体の生涯効用は次式で与えられる：
\[ U(c_1, c_2) = u(c_1) + \beta u(c_2). \tag{3}\]
右辺第２項だが，主体は若年期に意思決定を行うため，将来の消費 \(c_2\) から得られる効用 \(u(c_2)\) は割り引かれる．
ここで，効用関数について \(u'(c) > 0\) および \(u''(c) < 0\) を仮定し，主体は消費の平準化 (consumption smoothing) を望むものとする1．
以上の設定の下で解くべき意思決定問題は，次のように定式化できる．
\[ \begin{aligned} \max_{c_1, c_2, a} & \quad U(c_1, c_2) = u(c_1) + \beta u(c_2) \\ \text{s.t.} & \quad c_1 + a = w, \\ & \quad c_2 = (1+r)a. \\ \end{aligned} \tag{4}\]
ラグランジアンは
\[ \mathcal{L} = u(c_1) + \beta u(c_2) + \lambda_1 (w - c_1 - a) + \lambda_2[(1+r)a - c_2] \]
より，一階条件は
\[ \begin{aligned} &0 = \frac{\partial\mathcal{L}}{\partial c_1} = u'(c_1) - \lambda_1,\\ &0 = \frac{\partial\mathcal{L}}{\partial c_2} = \beta u'(c_2) - \lambda_2,\\ &0 = \frac{\partial\mathcal{L}}{\partial a} = -\lambda_1 + \lambda_2(1+r). \end{aligned} \]
整理すると 式 5 のオイラー条件（オイラー方程式）を得る．
\[ u'(c_1) = \beta (1+r) u'(c_2). \tag{5}\]
オイラー条件（式 5）の直感的な意味は，若年期と老年期で極端に消費水準が変動することを嫌うリスク回避的な主体は，若年期の消費 \(c_1\) から得られる限界効用 \(u'(c_1)\) と，老年期の消費 \(c_2\) から得られる限界効用の割引現在価値 \(\beta (1+r) u'(c_2)\) が一致するように消費計画を決めるということである．
ところで，何をもってモデルを「解いた」と言えるのか．マクロモデルはさまざまな経済変数が登場するため，この点を整理しておかないとすぐに迷子になる．
このモデルの内生変数は \(c_1\) と \(a\) であり，外生変数は \(w\) であるから2，モデルを解くにはある所得 \(w\) のもとでの主体の貯蓄関数 \(g(w)\) と消費関数 \(h(w)\) を導出すれば良い：
\[ \begin{aligned} a &= g(w), \\ c_1 &= h(w). \end{aligned} \]
この例のように，
\[ \text{内生変数} = f(\text{外生変数}) \]
と内生変数が外生変数の関数として書けるとき，「モデルが解けた」という．
２期間モデルは人間にとってはシンプルなモデルだが，コンピュータにとってはそうではない．
というのも，効用関数 \(u(c)\) という表現は抽象的で，コンピュータには理解できないからである．
そこで，コンピュータが理解できる形で効用関数を「特定化」する必要がある．
ここでは，マクロ経済学で頻繁に使われる相対的リスク回避度一定 (constant relative risk aversion: CRRA) 型効用関数を仮定する．
\[ u(c) = \frac{c^{1-\gamma}}{1-\gamma}. \tag{6}\]
リスクが存在するモデルであれば，\(\gamma\) は相対的リスク回避度 (coefficient of relative risk aversion) であり，同時に異時点間の代替の弾力性 (intertemporal elasticity of substitution) の逆数である．
式 6 をコンピュータが認識できるように Python で書くと，次のようになる．
import numpy as np
def CRRA(c, gamma):
if gamma != 1.0:
util = c**(1.0 - gamma) / (1.0 - gamma)
else:
util = np.log(c)
return util実際にプログラムを動かす際は，関数形の特定に加えて具体的にパラメータの値を決める必要があるが，このように関数形を特定してパラメータを定める一連の作業をカリブレーション (calibration) と呼ぶ．
異なる gamma の下で CRRA 型効用関数（式 6）を可視化したものが 図 1 である．
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
# c の範囲を指定して曲線をプロット
c = np.linspace(start=0.1, stop=1.0, num=100)
gammas = [1, 1.5, 2]
for gamma in gammas:
u = CRRA(c, gamma)
ax.plot(c, u, label=f'gamma = {gamma}')
ax.set_title("Utility Function for Different gamma values")
ax.set_xlabel(r"$c$")
ax.set_ylabel(r"$u(c)$")
ax.legend()
ax.grid()
pass
ベンチマーク・モデルでは人生を 2 期間に分けているので，モデル上の 1 期間を 30 年と想定する． そのため，割引因子 \(\beta\) と金利 \(r\) は年率ではなく 30 年間の値を使う．
今回は実際に 1 期間を 30 年でカリブレートしている Song, Storesletten, and Zilibotti (2012) の値を拝借する．割引因子は年率で \(\beta = 0.985\) として，1 期間は 30 年なので 30 乗する（\(\beta = 0.985^{30}\)）．金利は年率で \(2.5\%\) と設定すると \(1+r = 1.025^{30}\) となる．相対的リスク回避度は，よく使われる値の \(\gamma = 2\) としておく．
beta = 0.985**30 # 割引因子
gamma = 2.0 # 相対的危険回避度
r = 1.025**30 - 1.0 # 金利数値計算に入る前に，２期間モデルの解析解の性質を簡単に確認しておく．今回のモデルはシンプルなので，貯蓄関数 \(a = g(w)\) を手計算で求めることができる．
\[ % \begin{equation} a = \frac{w}{1+(1+r)\{\beta(1+r)\}^{-1/\gamma}} % \end{equation} \tag{7}\]
これがモデルから導出される真の貯蓄関数であり，内生変数 \(a\) が外生変数 \(w, r, \beta, \gamma\) で表されているので，これでモデルは解けている．
貯蓄関数（式 7）は若年期の所得 \(w\) の連続な線形の増加関数になっており，そのグラフを可視化したものが 図 2 である．
import japanize_matplotlib
# 傾き
slope = 1 / (1 + (1+r) * (beta*(1+r)) **(-1/gamma))
# 貯蓄関数
def saving(w, beta, r, gamma):
return w / (1 + (1+r) * (beta*(1+r)) **(-1/gamma))
fig, ax = plt.subplots()
w = np.linspace(start=0, stop=1.0, num=100)
a = saving(w, beta, r, gamma)
ax.plot(w, a, c="#FF7A72")
ax.set(title="貯蓄関数", xlabel="若年期の所得: "+r"$w$", ylabel="若年期の貯蓄: "+r"$a=g(w)$", xlim=(0,1), ylim=(0,0.4))
ax.grid()
pass
貯蓄関数（式 7）は連続関数であるが，コンピュータは連続という概念をそのままの形では理解できない．
そこで，数値計算においては基本的に，連続な変数を有限の \(N\) 個の点に離散化 (discretize) して考える必要がある．
若年期の所得 \(w\) がとりうる値は，\(w_i \in \{w_1,\ldots,w_N\}\) の範囲にあるとする．
コンピュータは無限を扱えないので，\(w \in [0, \infty)\) や \(w \in (-\infty, \infty)\) のような範囲を扱うことはできない． 数値計算において，定義域は必ず有界 (finite)である必要がある．
所得を \(N\) 個に離散化するということは，若年期の所得に応じて \(N\) 種類の経済主体が存在している状況を作り出すということである．
この離散的な点の集まりをグリッド (grid) あるいはノード (node) と呼ぶ． また，それぞれの点はグリッドポイント (grid point) あるいは評価点 (evaluation point) と呼ばれている．
ここでは単純に，若年期の所得 \(w\) は \([0.1, 1]\) 区間の間に \(0.1\) 刻みで \(10\) 個の点
\[ w_i \in \{0.1, 0.2,\ldots, 1.0\} \]
として存在していると考える．
nw = 10 # 所得グリッドの数
w_min = 0.1 # 所得の最大値
w_max = 1.0 # 所得の最小値
grid_w = np.linspace(start=0.1, stop=1.0, num=10)
print(grid_w)[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]この意思決定問題では \(w_i = 0\) を含めると若年期も老年期も消費ができなくなってしまうため，最小値は正値にしておく必要がある．
最適化問題を解く際にすでに決まっている若年期の所得 \(w\) を状態変数 (state variable) と呼び，意思決定において選択する若年期の貯蓄 \(a\) を制御変数 (control variable) と呼ぶ． また，貯蓄関数のような意思決定関数は政策関数 (policy function) と呼ばれる．
まずは，所得 \(w\) と資産 \(a\) がともに離散的な場合に２期間モデルを解く方法を見る．
すなわち，状態変数 \(w\) だけでなく，制御変数 \(a\) についても離散化を行う．
ここでは
\[ a_j \in \{0.025, 0.05, 0.075, 0.1, \ldots, 1.0\} \]
と \(0.025\) 刻みで \(40\) 個のグリッドをとる3．
na = 40 # 貯蓄グリッドの数
a_min = 0.025 # 貯蓄の最小値
a_max = 1.0 # 貯蓄の最大値
grid_a = np.linspace(start=a_min, stop=a_max, num=na)
print(grid_a)[0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3
0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6
0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9
0.925 0.95 0.975 1. ]このような離散化を行うと，解くべき効用最大化問題（式 4）は，所与の \(w_i\) について，
\[ \max_{a_j} \quad \frac{[w_i-a_j]^{1-\gamma}}{1-\gamma} + \beta\frac{[(1+r)a_j]^{1-\gamma}}{1-\gamma} \tag{8}\]
と書き直すことができる．
組み合わせの数は，所得 \(w_i\) と資産 \(a_j\) がそれぞれ \(10\) 通りと \(40\) 通りあるので， \(400\) 通りとなる．
そして各所得 \(w_i\) のもとで \(40\) 種類の資産 \(a_j\) の中から生涯効用（式 8）を最大化する資産 \(a^*\) を求めれば良いので，これくらいであれば「総当たり」で求めることができる．
obj = np.zeros((na, nw))
# 状態変数と制御変数のすべての組み合わせ (w,a) について生涯効用を計算
for i in range(nw):
for j in range(na):
c = grid_w[i] - grid_a[j]
if c > 0:
obj[j, i] = CRRA(c, gamma) + beta * CRRA((1+r)*grid_a[j], gamma)
else:
obj[j, i] = -10000.0 na \(\times\) nw ゼロ行列
試しに，w = 0.5, 0.8, 1.0 のときの生涯効用を可視化してみよう（図 3）．
fig, ax = plt.subplots()
ax.plot(grid_a, obj[:,4], c="b", label="w=0.5")
ax.plot(grid_a, obj[:,7], c="r", label="w=0.8")
ax.plot(grid_a, obj[:,9], c="g", label="w=1.0")
ax.set(xlabel="若年期の貯蓄(老年期の資産): "+r"$a$", ylabel="生涯効用: "+r"$U(c_1,c_2)$", xlim=(0,1), ylim=(-10,0))
ax.legend()
ax.grid(ls="--")
pass
図 3 では若年期の所得 \(w_i\) の下での生涯効用 \(U(a_j; w_i)\) が凹関数として可視化されている．
直感的にそれぞれの山の頂上に対応する貯蓄点が，求めたい最適貯蓄だとわかるだろう．
そこで，若年期の所得 \(w_i\) ごとに最適な貯蓄点 \(a_j\) を計算し，求めた貯蓄関数を可視化する（図 4）．
pol = np.zeros(nw)
for i in range(nw):
maximizer = np.argmax(obj[:, i])
pol[i] = grid_a[maximizer]
fig, ax = plt.subplots()
ax.plot(grid_w, saving(grid_w, beta, r, gamma), c="#FF7A72",label="解析解")
ax.plot(grid_w, pol, c="#78C2AD", marker="o", label="グリッドサーチ")
ax.set(title="貯蓄関数", xlabel="若年期の所得: "+r"$w$", ylabel="若年期の貯蓄: "+r"$a=g(w)$")
ax.grid(ls="--")
ax.legend()
pass
図 4 を見ると，グリッドサーチで求めた最適貯蓄（緑色の点）が，解析解（赤色の直線）からあまり離れていないことがわかる．
グリッドサーチのグラフは点同士を直線でつないだ折れ線グラフのように見えるが，これは plot() 関数の仕様によるものである．実際には，状態変数 \(w\) のグリッドポイント上でしか最適貯蓄 \(a\) を求めていないことに注意されたい．
グリッドサーチで精度を上げるには，離散化するグリッドの数を増やせば良さそうだが，考えているモデルによっては次のセクションで扱う「次元の呪い」という問題にぶつかる．
状態変数と制御変数をどちらも離散化して，とりうる全ての組み合わせを計算するというグリッドサーチのアプローチは，単純ではあるが複雑な非線形モデルにも使えるのでなかなか侮れない．
たとえば，政策関数の形状や性質がよくわかっていないモデルを解く場合，とりあえずグリッドを切ってモデルの性質を大雑把に掴むというアプローチは有効である．
しかし，このアプローチは次元の呪い (curse of dimensionality) の影響をダイレクトに受けるので，計算に時間がかかる傾向にある．
状態変数の数が \(n\) 個のとき，最適貯蓄を計算する必要のあるグリッドポイントの数は
\[ (\text{状態変数のグリッド数})^n \]
で与えられる．したがって，たとえば教育水準の違いや性別の違いを考慮できるように若年期における状態変数を増やしていくと，計算時間は指数的に増えてしまう．
状態変数を増やせば増やすほど，計算時間は指数的に増える．そのため，現実的な計算時間で複雑なモデルを高い精度で解くための工夫が常に求められる．
数値計算では，コンピュータが理解できるように関数を特定化してパラメータを設定するカリブレーションという作業が必要なのであった．
ここまでの実装では，逐一関数を定義してパラメータを指定していた．
しかし，その方法ではパラメータが変わった時にモデルの振る舞いがどう変わるかを分析する比較静学のような場面で苦労しそうなのは容易に想像できる．
そこで，オブジェクト指向プログラミングのクラス (class) という機能でモデルのテンプレートを作っておくと便利である．
```{python}
class Models:
def __init__(self, beta, gamma, r, nw, w_min, w_max, grid_w, na, a_min, a_max, grid_a):
self.beta = beta
self.gamma = gamma
self.r = r
self.nw = nw
self.w_min = w_min
self.w_max = w_max
self.grid_w = grid_w
self.na = na
self.a_min = a_min
self.a_max = a_max
self.grid_a = grid_a
```カリブレーションを行う際は Models を定義した上で，次の Calibration() 関数を実行する．
```{python}
def Calibration():
beta = 0.985 ** 30
gamma = 2.0
r = 1.025**30 - 1.0
nw = 10
w_min = 0.1
w_max = 1.0
na = 40
a_min = 0.025
a_max = 1.0
grid_w = np.linspace(w_min, w_max, nw)
grid_a = np.linspace(a_min, a_max, na)
return Models(beta, gamma, r, nw, w_min, w_max, grid_w, na, a_min, a_max, grid_a)
```このノートでもこれ以降は，クラスを使ってさまざまなパラメータを管理する．
グリッドサーチはモデルを拡張した場合に計算量が指数的に増加するのに加え，グリッド数を節約すると精度が悪化するという問題を抱えていた．
実際，真の貯蓄関数は 図 2 のように線形である一方，グリッドサーチで求めた貯蓄関数は 図 4 にあるように直線ではない．
そこで，状態変数 \(w\) は § 3.2 と同様に離散化するが，連続な制御変数 \(a\) を許容するより洗練されたアプローチを採用しよう．
グリッドサーチで考えた問題（式 8）で離散化していた \(a_j\) を連続値にするので，今考えている最適化問題は次式で表される．
\[ \max_{a\in\R} \quad \frac{[w_i-a]^{1-\gamma}}{1-\gamma} + \beta\frac{[(1+r)a]^{1-\gamma}}{1-\gamma} \tag{9}\]
式 9 のような最適化問題を数値計算で解くには，数値計算ソフトに標準的に実装されている最適化ライブラリを用いれば良い．
大抵は，目的関数とその関数のパラメータを入力として受け取り，目的関数の最大値や最小値を探索する仕様になっている．
最適化問題（式 9）であれば，カリブレーションしたパラメータ \(\{\beta, \gamma, r\}\) を所与として，各状態変数 \(w_i \in \Brace{w_1, \ldots, w_N}\) の下で生涯効用の最大値とそれを与える最大元 \(a\) を探索してくれる．
まず，カリブレーションを § 3.2.2 で導入した方法で行う．
class Models:
def __init__(self, beta, gamma, r, nw, w_min, w_max, grid_w):
self.beta = beta
self.gamma = gamma
self.r = r
self.nw = nw
self.w_min = w_min
self.w_max = w_max
self.grid_w = grid_w
def Calibration():
beta = 0.985 ** 30
gamma = 2.0
r = 1.025**30 - 1.0
nw = 10
w_min = 0.1
w_max = 1.0
grid_w = np.linspace(w_min, w_max, nw)
return Models(beta, gamma, r, nw, w_min, w_max, grid_w)Calibration() 関数を実行することで，Calibration() 内で定義されたパラメータを情報として持つ Models インスタンスが生成されるので，params に代入しておく．
params = Calibration()ここでは，scipy の optimize モジュールにある fminbound() 関数を使って最適化問題（式 9）を解く4．
from scipy import optimizeでは，optimize.fminbound() に渡す目的関数を定義しよう．最小化問題に変換するために \(-1\) をかけている．
\[ \text{obj}(a, w_i;\beta, \gamma, r) = -\Brace{\frac{[w_i-a]^{1-\gamma}}{1-\gamma} + \beta\frac{[(1+r)a]^{1-\gamma}}{1-\gamma}} \]
def obj(a, w_i, params):
c = w_i - a
if c > 0:
life_util = CRRA(c, params.gamma) + params.beta * CRRA((1+params.r)*a, params.gamma)
else:
life_util = -100000.0
return -1.0 * life_util上の obj() は a, w_i, params の３つの引数を持つが，実際に optimize.fminbound() に渡すときは，最適化を行う a だけの関数にしておく必要がある．
所与の状態変数 \(w_i\) に対して制御変数 \(a\) についての最適化問題を解けば良いので，grid_w についての for ループ内で optimize.fminbound() を使えば良い．
opt_a = np.zeros(params.nw)
for i, w_i in enumerate(params.grid_w):
obj_specified = lambda a: obj(a, w_i, params)
opt_a[i] = optimize.fminbound(obj_specified, w_i*0.01, w_i*2.0)
print(f"最適貯蓄の配列: \n{opt_a}")optimize.fminbound() に渡すために，obj() を a だけの関数にする
最適貯蓄の配列:
[0.03550017 0.07100034 0.10650385 0.14200403 0.17750421 0.21300505
0.24850589 0.28400673 0.31950757 0.35500841]opt_a には grid_w の各 \(w_i\) に対する最適貯蓄 \(a\) が格納してある．これを可視化したものが 図 5 である．
fig, ax = plt.subplots()
ax.grid(ls="--")
ax.scatter(params.grid_w, saving(params.grid_w, params.beta, params.r, params.gamma), c="#FF7A72",label="解析解")
ax.plot(params.grid_w, opt_a, c="#78C2AD",label="アルゴリズム")
ax.set(title="貯蓄関数", xlabel="若年期の所得: "+r"$w$", ylabel="若年期の貯蓄: "+r"$a=g(w)$")
ax.legend()
pass
同じグリッドポイントを使っているが，グリッドサーチで求めた 図 4 と違って，最適化アルゴリズムで導出した 図 5 は解析的な解（図 2）が示す綺麗な直線になっており，計算精度が大幅に改善された様子がわかる5．
これは制御変数 \(a\) がとりうる値を連続値にしたことに起因する．
モデルの解が満たすべきオイラー条件（式 5）をうまく使って数値計算を行うアプローチもある．
ここでは，オイラー条件を求根問題に落とし込んで解く方法（§ 3.4.1）と，政策関数自体をパラメトリックに近似する方法（§ 3.4.2）を扱う．
オイラー条件（式 5）に予算制約を代入すると次式を得る．
\[ u'(w-a) = \beta(1+r)u'((1+r)a) \tag{10}\]
状態変数 \(w\) を \(w_i\) と離散化して，それぞれの変数の役割を見てみると
\[ u'(\eqDescribe{w_i}{given} - \eqDescribe{a}{control}) = \eqDescribe{\beta(1+r)}{parameter}u'(\eqDescribe{(1+r)}{parameter}\ \eqDescribe{a}{control}) \]
なので，未知変数は \(a\) だけである．そこで，式 10 を変形して次式で残差関数 (residual function) を定義する：
\[ R(a;w_i) \equiv \beta(1+r)\frac{u'((1+r)a)}{u'(w_i-a)} - 1. \tag{11}\]
式 11 を使うと，オイラー条件を満たす制御変数 \(a\) を見つける問題を，
\[ R(a;w_i) = 0 \]
となる \(a\) を探すという求根（ゼロ点）問題 (root-finding problem) に変換することができる．
一般に，オイラー条件を変換して得た残差関数（式 11）は，複雑な形をした非線形方程式である可能性があるが，非線形方程式のゼロ点を探すアルゴリズムの研究は長い歴史を持つため，様々なアプローチが考案されている．
求根問題は scipy.optimize.fsolve() 関数で解くことができる． ここではカリブレーションは § 3.3 のものを再び使用する．
まずはソルバーに渡す残差関数（式 11）を実装しよう．残差関数内にある限界効用 (marginal utility) は
\[ u'(c) = c^{-\gamma} \]
で与えられるので，Python で表すと次のようになる．
def marginal_util(c, gamma):
return c ** (- gamma)これを用いると残差関数は次のように実装できる．
def resid(a, w_i, params):
c = w_i - a
if c > 0:
mu_y = marginal_util(c, params.gamma) # 若年期の限界効用
else:
mu_y = 10000.0
mu_o = marginal_util((1+params.r)*a, params.gamma) # 老年期の限界効用
return params.beta * (1+params.r) * (mu_o/mu_y) - 1.0§ 3.3 と同様に，ソルバーを使って各状態変数 \(w_i\) の下で最適貯蓄を求めていく．
opt_a = np.zeros(params.nw)
for i, w_i in enumerate(params.grid_w):
resid_specified = lambda a: resid(a, w_i, params)
opt_a[i] = optimize.fsolve(resid_specified, x0=0.01)[0]
print(f"最適貯蓄の配列: \n{opt_a}")optimize.fsolve() に渡すために，resid() を a だけの関数にする
optimize.fsolve() の返り値の最初の要素に格納されている
最適貯蓄の配列:
[0.03550089 0.07100178 0.10650266 0.14200355 0.17750444 0.21300533
0.24850621 0.2840071 0.31950799 0.35500888]結果を図示した 図 6 を見ると，アルゴリズムで求めた解が解析解の直線上に乗っている様子がわかる．
fig, ax = plt.subplots()
ax.grid(ls="--")
ax.plot(params.grid_w, saving(params.grid_w, params.beta, params.r, params.gamma), c="#FF7A72",label="解析解")
ax.scatter(params.grid_w, opt_a, c="#78C2AD", label="アルゴリズム")
ax.set(title="貯蓄関数", xlabel="若年期の所得: "+r"$w$", ylabel="若年期の貯蓄: "+r"$a=g(w)$")
ax.legend()
pass
最適化（§ 3.3）と求根アルゴリズム（§ 3.4.1）を使った手法はどちらも，本来連続値をとる現在の所得水準 \(w\) を有限個に離散化して，そのグリッドポイント上で最適貯蓄を計算するという点は共通している．
これに対して，射影法 (projection method) では求めたい政策関数そのものをパラメトリックに近似するというアプローチをとる6．
今考えている政策関数は若年期の所得 \(w\) を変数にとる貯蓄関数 \(a=g(w)\) である．貯蓄関数をパラメトリックに近似するとは、貯蓄関数を基底関数 (basis function) \(\{\Psi_m\}_{m=0}^M\) の線形結合 \(\hat{g}(w;\bm{\theta})\) で近似するということである7．
\[ a = g(w) \approx \hat{g}(w;\bm{\theta}) = \sum_{m=0}^M \theta_m\Psi_m(w) \tag{12}\]
ここで \(\bm{\theta} \equiv \Paren{\theta_m}_{m=0}^M\) は未知の \(M+1\) 次元係数ベクトルであり，推定する対象である．
この段階では，未知の貯蓄関数 \(g(w)\) はきっとパラメトリックな関数 \(\hat{g}(w;\bm{\theta})\) で表されるだろうと思っているに過ぎず，何も解決していない．
しかし，真の貯蓄関数 \(g(w)\) はオイラー条件（式 10）を満たすので，それを近似した関数 \(\hat{g}(w;\bm{\theta})\) もできるだけオイラー条件を満たしておいて欲しいと思うのは自然な要請である．
つまり射影法とは，オイラー条件を満たす政策関数を探す問題を，近似関数 \(\hat{g}(w;\bm{\theta})\) ができるだけオイラー条件を満たすように係数ベクトル \(\bm{\theta}\) を決める問題に置き換える手法である．
以下では基底関数を \(\Psi_m(w) = w^m\) と指定し，貯蓄関数を次の多項式 (polynomial) で近似して議論を進める．
\[ a \approx \hat{g}(w;\bm{\theta}) = \sum_{m=0}^M \theta_m w^m \tag{13}\]
ところで，上述の「できるだけオイラー条件を満たす」ことをどうやって評価すれば良いだろうか．ここで，§ 3.4.1 で定義した残差関数（式 11）を思い出そう．
仮に近似関数 \(\hat{g}(w;\bm{\theta})\) が真の貯蓄関数 \(g(w)\) と完全に一致する場合は残差はゼロになり，うまく近似しているのであれば残差はゼロに近いはずである．つまり，
\[ R(\bm{\theta}; w) \equiv \beta (1+r)\frac{u'((1+r)\hat{g}(w;\bm{\theta}))}{u'(w-\hat{g}(w;\bm{\theta}))} - 1 \approx 0 \tag{14}\]
が「あらゆる \(w\) で」成り立つ．
「」をつけたのは，実際に数値計算を行う際コンピュータは連続的な \(w\) を扱えないため，任意にとった評価点 \(\{w_i\}_{i=1}^N\) 上での残差がゼロに限りなく近くなるようなベクトル \(\bm{\theta}\) を見つける必要があるからである．
要素に評価点 \(w_i\) 上での残差 \(R(\bm{\theta}; w_i)\) を持つベクトルを \(\bm{R}(\bm{w}; \bm{\theta}) \equiv \Paren{R(\bm{\theta}; w_i)}_{i=1}^N\) とする． ただし，\(\bm{w} = \Paren{w_i}_{i=1}^N\) である．
\(\rho(\cdot, \cdot)\) を距離関数 (metric function) とすると，考えている問題は
\[ \bm{\theta}^* = \argmin_\bm{\theta} \rho(\bm{R}(\bm{w}; \bm{\theta}), \bm{0}) \tag{15}\]
と定式化できる．このように，評価点上でのみ距離を測る方法を選点法 (collocation method) と呼ぶ．
カリブレーションは § 3.3 のものを再び使用する．
ここでは多項式 式 13 の次数を \(M=1\) として，１次関数 \(\hat{g}(w;\bm{\theta}) = \theta_0 + \theta_1 w\) で政策関数を近似しよう．
各評価点 \(w_i\) での近似値 \(\hat{g}(w_i; \bm{\theta})\) を要素にもつ列ベクトルを \(\widehat{\bm{g}}\) とすると
\[ \widehat{\bm{g}} = \Paren{\hat{g}(w_i;\bm{\theta})}_{i=1}^N = \begin{pmatrix} \theta_0 + \theta_1 w_1 \\ \vdots \\ \theta_0 + \theta_1 w_i \\ \vdots \\ \theta_0 + \theta_1 w_N \end{pmatrix} = \eqDescribe{\begin{pmatrix} 1 & w_1 \\ \vdots & \vdots \\ 1 & w_i \\ \vdots & \vdots \\ 1 & w_N \end{pmatrix}}{denoted by X below} \begin{pmatrix} \theta_0 \\ \theta_1 \end{pmatrix} \]
と書けるが，これを実装すると次のようになる．
def approx_g(theta: np.ndarray, w: np.ndarray) -> np.ndarray:
dim = len(theta)
nw = len(w)
X = np.zeros((nw, dim))
for j in range(dim):
X[:, j] = w ** j
return X @ theta次に，残差関数 resid()（Code 1）を利用して，式 14 に従って残差ベクトル \(\bm{R}(\bm{w}; \bm{\theta})\) を実装する．
def resid_vec(theta, params):
g_hats = approx_g(theta, params.grid_w)
R = np.zeros(nw)
for g, w_i, i in zip(g_hats, params.grid_w, np.arange(nw)):
R[i] = resid(g, w_i, params)
return Rここで，距離関数 \(\rho\) にユークリッド距離を採用すると，考えている問題（式 15）は残差二乗和を最小にする \(\bm{\theta}\) を見つける非線形最小二乗法に帰着する．
非線形最小二乗法は scipy.optimize.least_squares() で実装することができる．fun 引数に残差ベクトルを計算する関数を渡せば良い．
def projection(params, initial_guess=[0.1, 0.35]):
resid_specified = lambda theta: resid_vec(theta, params)
result = optimize.least_squares(fun=resid_specified, x0=initial_guess, method="lm")
estimate_g = approx_g(result.x, params.grid_w)
return result.x, result.success, estimate_goptimize.least_squares() に渡すために，resid_vec() を theta だけの関数にする
では，実装した projection() 関数を使って射影法を実際に行おう．
result = projection(params)
print(f"convergence: {result[1]}")
print(f"The estimated parameter: {result[0]}")convergence: True
The estimated parameter: [-1.90171907e-10 3.55008878e-01]解析解（式 7）によれば，切片はゼロで，傾きは \(1/1+(1+r)\{\beta(1+r)\}^{-1/\gamma}\) から計算すると 0.355 であるから，非常に精度良く近似することができている．
推定されたパラメータ \(\widehat{\bm{\theta}}\) の下での近似関数のグリッド \(\widehat{\bm{g}}(\bm{w}; \widehat{\bm{\theta}})\) と真の貯蓄関数を可視化したものが 図 7 である．
fig, ax = plt.subplots()
ax.grid(ls="--")
ax.plot(params.grid_w, saving(params.grid_w, params.beta, params.r, params.gamma), c="#FF7A72",label="解析解")
ax.scatter(params.grid_w, result[2], c="#78C2AD", label="射影法")
ax.set(title="貯蓄関数", xlabel="若年期の所得: "+r"$w$", ylabel="若年期の貯蓄: "+r"$a=g(w)$")
ax.legend()
pass
効用関数が凹関数であるとき，経済主体はリスク回避的である．↩︎
経済学では，方程式の解として決まる変数を内生変数 (endogenous variable)，方程式を解く前に値がすでに決まっている変数を外生変数 (exogenous variable) と呼ぶ(神谷和也 and 浦井憲 1996)．外生変数は単にパラメータとか状態変数とも言われる．↩︎
\(w_i\) と同じ離散化をすると，\(w = 0.1\) のときに 1 期目の消費 \(c_1 = w - a\) がゼロか負値しかとれず，最適貯蓄あるいは最適消費が存在しなくなってしまうため．↩︎
fminbound() のように数値計算によって方程式を解くアルゴリズムのことをソルバー (solver) と呼ぶ．↩︎
図 5 では，アルゴリズムで求めた解のグラフが直線に見えるが，これは plot() の仕様によるもので，実際にはグリッドポイント上でしか最適貯蓄を計算していないことに注意．↩︎
筆者は射影法を初めて学んだとき，未知の関数の形状をパラメータで規定し，ある判断基準に基づいてパラメータを推定するという点にどこか計量経済学のような雰囲気を感じた．↩︎
チェビシェフ多項式 (Chebyshev polynomials) が基底関数として頻繁に用いられている．↩︎