動学モデルの基礎

Author

Affiliation

柳本和春

神戸大学

Published

October 20, 2025

Modified

October 27, 2025

using Plots
using LaTeXStrings
using Distributions
import Random
using Roots
using Optim
using LinearAlgebra
using Interpolations
using SpecialFunctions
default(size=(500, 309), titlefontsize=10, fmt=:svg)

グリッドと補完

未知なる関数 \(f(x)\) の一般形が分からないが, 具体的な \(x\) に関しては \(f(x)\) の値が計算できる状況というのがよくあります. ちょっと極端な例を考えてみましょう. 次の効用最大化問題を考えます.

\[ \max_{n} \log c - \Phi(n) \text{ s.t.　}　c = w n. \]

ここで, 労働供給の不効用 \(\Phi(n) = \phi \log (n!)\) とします. この時の最適な \(n^*\) を求めてみましょう.

Φ(n::Int; ϕ=0.1) = ϕ * sum(log(i) for i in 1:n)
u(n::Int; w=1.) = log(n) - Φ(n)

Grid Search

最も単純な方法は, \(n = 1, \dots, 10\) ぐらいまで試してみて, 最も効用 \(u(c, n)\) が大きくなる \(n\) を選ぶ方法です. これをグリッドサーチと呼びます.

Code

scatter(1:10, u, label=false, xlabel=L"n", ylabel=L"u(c, n)")

Figure 1: Grid Search for Optimal n. \(ϕ = 0.1\).

argmax(u, 1:10)

補完 (Interpolation)

これを実数の世界に拡張してみましょう. もし, \(n\) が実数の時の \(\Phi(n)\) の値が計算できれば, 通常の一変数関数の最大化問題として解くことができます. そのために, \(\Phi(n)\) を補完 (interpolation) してみましょう. Juliaでは Interpolations.jl パッケージを用います.

n_grid = 1:10
Φ_grid = Φ.(n_grid)
Φ_interp = linear_interpolation(n_grid, Φ_grid)
Φ_spline = Interpolations.scale(
    interpolate(Φ_grid, BSpline(Cubic(Line(OnGrid())))),
    n_grid
)

Code

scatter(1:10, Φ, label="Grid Points", xlabel=L"n", ylabel=L"\Phi(n)")
plot!(1:0.1:10, n -> Φ_interp(n), label="Linear Interpolation", legend=:topleft)
plot!(1:0.1:10, n -> Φ_spline(n), label="Cubic Spline", linestyle=:dash)

Figure 2: Linear and Spline Interpolation of Cost Function \(\Phi(n)\).

ここでは線形補間 (linear interpolation) と3次スプライン補間 (cubic spline interpolation) の両方を定義しました. 線形補間は2点間を直線で結ぶだけですが, スプライン補間はより滑らかな曲線で補間します. 詳細は北尾, 砂川, and 山田 (2024) の付録Bが参考になります.

このように実数上に定義された関数を用いれば, 通常の最適化問題のように解くことができます.¹

u_interp(n) = log(n) - Φ_interp(n)
u_spline(n) = log(n) - Φ_spline(n)

sol_interp = optimize(n -> -u_interp(n), 1.0, 10.0)
sol_spline = optimize(n -> -u_spline(n), 1.0, 10.0)

println("Linear Interpolation: n = $(round(sol_interp.minimizer, digits=4))")
println("Cubic Spline: n = $(round(sol_spline.minimizer, digits=4))")

Linear Interpolation: n = 5.5811
Cubic Spline: n = 5.5516

今回の例では, 階乗はガンマ関数を用いて実数に拡張できるので, より厳密な解を求めることもできます.

Code

Φ(n::Real; ϕ=0.1) = ϕ * loggamma(n + 1)
u_exact(n) = log(n) - Φ(n)
sol_exact = optimize(n -> -u_exact(n), 1.0, 10.0)
println("Exact Solution: n = $(round(sol_exact.minimizer, digits=4))")

Exact Solution: n = 5.5512

ベルマン方程式の数値解法

以下の無限期間代表的個人のモデルを考えます.

\[ \max_{\{c_t, k_{t+1}\}_{t=0}^{\infty}} \sum_{t=0}^{\infty} \beta^t u(c_t) \quad \text{s.t.} \quad c_t + k_{t+1} = f(k_t) + (1-\delta)k_t. \]

これをベルマン方程式 \(V(k)\) に書き換えた時の数値解法を考えます.

\[ V(k) = \max_{k'} \left\{u(f(k) + (1-\delta)k - k') + \beta V(k')\right\} \]

大きく分けて, 2つの方法があります.

Value Function Iteration
- ベルマン方程式の不動点を直接求める方法
- 収束性が保証されているが, 収束速度は比較的遅い
Policy Function Iteration
- 政策関数の不動点を求める方法
- 収束性は保証されていないが, 収束速度は比較的速い

数値計算に用いるために, 以下の関数形とパラメータを仮定します.

\(u(c) = \log c\)
\(f(k) = k^\alpha\)
\(\beta = 0.96\)
\(\delta = 0.1\)
\(\alpha = 0.36\)

@kwdef struct Model{TF<:AbstractFloat,TI<:Integer}
    # Parameters
    α::TF = 0.4
    β::TF = 0.96
    δ::TF = 0.1

    # Value Function
    n_k::TI = 101
    k_min::TF = 0.05
    k_max::TF = 0.5
    k_grid::Vector{TF} = collect(range(k_min, k_max, length=n_k))
    V::Vector{TF} = zeros(n_k)
    h::Vector{TF} = copy(k_grid)
end

u(c; m) = log(c)
u′(c; m) = 1 / c
f(k; m) = k^m.α
f′(k; m) = m.α * k^(m.α - 1)

Value Function Iteration

\(k\) の定義域を \([k_{\min}, k_{\max}]\) を \(n\) 個のグリッドポイントに分割する.
\(V(k)\) の初期値 \(V^0\) を適当に設定する. つまり, \(V^0 = (V_0^0, V_1^0, \dots, V_n^0)\) に適当な値を設定する.
\(V^0\) をベルマン方程式の右辺に代入し, 各 \(k_i\) に対して \(V_i^1\) を計算する.
2-3を繰り返す. 収束条件は以下のように設定する.
- \(|V^N - V^{N-1}| < \varepsilon\)
- \(|V^N - V^{N-1}|\) はベクトルのノルム
- \(\varepsilon\) は十分小さな値を設定する. 例えば, \(\epsilon = 10^{-6}\) とする

なおこのような収束の閾値 \(\varepsilon\) を tolerance と呼ぶことがあります.

function vfi!(m::Model; tol=1e-6, max_iter=1000, verbose=true)
    (; k_grid, β, δ) = m

    iter, dist = 0, Inf
    V_new = similar(m.V)
    while dist > tol && iter < max_iter
        for (i_k, k) in enumerate(k_grid)
            c_min = 1e-9 # consumption cannot be negative
            V_new[i_k] = maximum(
                log(max(f(k; m) + (1 - δ) * k - k′, c_min)) + β * m.V[i_k′]
                for (i_k′, k′) in enumerate(k_grid)
            )
        end

        dist = maximum(abs, V_new - m.V)
        m.V .= V_new
        iter += 1
    end

    if verbose
        if iter == max_iter
            println("Warning: Maximum iterations reached.")
        else
            println("Converged in $iter iterations.")
        end
    end

    return nothing
end

m = Model()
vfi!(m)
plot(m.k_grid, m.V, label=false, xlabel=L"Capital $k$", ylabel=L"Value $V(k)$")

Converged in 315 iterations.

Figure 3: Value Function Iteration

なお, 各イタレーションでは, \(k\) に対する最適な \(k'\) や \(c\) を求めています. それらを記録しておくことで, 政策関数 \(c = h(k)\) や \(k' = g(k)\) を求めることができます.

Policy Function Iteration

上記のベルマン方程式から, 以下のオイラー方程式が導出されます.

\[ u'(c) = \beta u'(c')\left(f'(k') + (1-\delta)\right). \]

ここで, 政策関数 \(c' = h(k')\) と予算制約 \(k' = f(k) + (1-\delta)k - c\) を導入すると, 上記のオイラー方程式は以下のように書き換えられます.

\[ u'(c) = \beta u'(h(f(k) + (1-\delta)k - c))\left(f'(f(k) + (1-\delta)k - c) + (1-\delta)\right). \]

この時, \(h(k)\) を所与とした時に \(c\) について解くと新しい政策関数 \(h(k)\) が得られます. この政策関数繰り返し更新することで, 政策関数の不動点を求めることができます.

Policy Function Iteration

\(k\) の定義域を \([k_{\min}, k_{\max}]\) を \(n\) 個のグリッドポイントに分割する
\(h(k)\) の初期値 \(h^0\) を適当に設定する. つまり, \(h^0 = (h_0^0, h_1^0, \dots, h_n^0)\) に適当な値を設定する.
\(h^0\) をオイラー方程式の右辺に代入し, 各 \(k_i\) に対して \(h_i^1\) を計算する.
2-3を繰り返す. 収束条件は以下のように設定する.
- \(|h^N - h^{N-1}| < \varepsilon\)
- \(|h^N - h^{N-1}|\) はベクトルのノルム
- \(\varepsilon\) は十分小さな値を設定する. 例えば, \(\epsilon = 10^{-6}\) とする

function euler_eq(c, k, h; m)
    (; β, δ, k_grid) = m

    k′ = max(f(k; m) + (1 - δ) * k - c, k_grid[begin])
    h_interp = linear_interpolation(k_grid, h)
    c′ = h_interp(k′)

    return u′(c; m) - β * u′(c′; m) * (f′(k′; m) + 1 - δ)
end

function pfi!(m::Model; tol=1e-6, max_iter=1000, verbose=true)
    (; k_min, k_max, k_grid, α, β, δ) = m

    h_new = similar(m.h)
    iter, dist = 0, Inf
    while dist > tol && iter < max_iter

        for (i_k, k) in enumerate(k_grid)
            c_min = max(1e-9, f(k; m) + (1 - δ) * k - k_grid[end])
            c_max = f(k; m) + (1 - δ) * k - k_grid[begin]
            ee_left = euler_eq(c_min, k, m.h; m)
            ee_right = euler_eq(c_max, k, m.h; m)
            if ee_left > 0 && ee_right > 0
                h_new[i_k] = c_max
            elseif ee_left < 0 && ee_right < 0
                h_new[i_k] = c_min
            else
                h_new[i_k] = find_zero(c -> euler_eq(c, k, m.h; m), (c_min, c_max))
            end
        end

        dist = maximum(abs, h_new .- m.h)
        m.h .= h_new
        iter += 1
    end

    if verbose
        if iter == max_iter
            println("Warning: Maximum iterations reached.")
        else
            println("Converged in $iter iterations.")
        end
    end

    return nothing
end

m = Model()
pfi!(m)
plot(m.k_grid, m.h, label=false, xlabel=L"Capital $k$", ylabel=L"Policy $h(k)$")

Converged in 5 iterations.

Figure 4: Policy Function Iteration.

ここから \(k' := g(k) = f(k) + (1-\delta)k - h(k)\) が求められます. さらに, \(h(k), g(k)\) から価値関数 \(V(k)\) を求める方法は主に2つあります.

\(V(k) \leftarrow u(h(k)) + \beta V(g(k))\) を用いて \(V(k)\) の不動点を求める方法 (VFI)
線形方程式
1. 離散化された \(k\) から \(k'\) への遷移行列 \(P\) を求める. グリッド \(k_i\) が \(k_j\) と \(k_{j+1}\) の間に遷移する場合, 以下のように定義します.
  - \(P_{i, j} = \frac{k_{j+1} - g(k_i)}{k_{j+1} - k_j}\)
  - \(P_{i, j+1} = \frac{g(k_i) - k_j}{k_{j+1} - k_j}\)
  - \(P_{i, k} = 0\), \(k \neq j, j+1\)
2. \(u, V\) を \(u(h(k)), V(k)\) からなる離散化されたベクトルとすると, \(V = u + \beta P V\) となる
3. \(V = (I - \beta P)^{-1} u\) により \(V\) を求めることができる

ここでは, より高速な線型方程式を用いる方法を実装します. より速度を求める場合, 遷移行列はほとんどの要素がゼロであるため, 疎行列 (Sparse Matrix) を用いることで高速化する可能性があります. Juliaの場合, SparseArrays.jl で実装されています.

function compute_vf!(m::Model)
    (; h, n_k, k_grid, h, β, δ) = m
    P = zeros(n_k, n_k)
    for (i_k, k) in enumerate(k_grid)
        g = f(k; m) + (1 - δ) * k - h[i_k]
        j = searchsortedlast(k_grid, g)
        if j == 1 || j == n_k
            P[i_k, j] = 1.0
        else
            P[i_k, j] = (k_grid[j+1] - g) / (k_grid[j+1] - k_grid[j])
            P[i_k, j+1] = (g - k_grid[j]) / (k_grid[j+1] - k_grid[j])
        end
    end

    m.V .= (I - β * P) \ (u.(h; m))

    return nothing
end

compute_vf!(m)
plot(m.k_grid, m.V, label=false, xlabel=L"Capital $k$", ylabel=L"Value $V(k)$")

Figure 5: Value Function from PFI and Linear Equation

また, VFIよりもPFIの方が収束速度が速いことも確認できますが, 実行速度も速いことも確認できます.

using BenchmarkTools
@benchmark vfi!(m, verbose=false) setup = (m = Model())

BenchmarkTools.Trial: 66 samples with 1 evaluation per sample.
 Range (min … max):  74.393 ms …  77.433 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     76.274 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   76.230 ms ± 561.599 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

                              ▅    ▂▂    █▂ ▂█ ▂▂    ▂          
  ▅▁▁▁▁▁▁▁▁▁▁▅▁▁▁▁▅▅▁▁▁█▁▁▅▅▁▁██▅▁▁████▅▅█████▁██▅▁▅▅█▅█▅▁▁█▅▅ ▁
  74.4 ms         Histogram: frequency by time         77.2 ms <

 Memory estimate: 286.38 KiB, allocs estimate: 632.

@benchmark (pfi!(m, verbose=false); compute_vf!(m)) setup = (m = Model())

BenchmarkTools.Trial: 917 samples with 1 evaluation per sample.
 Range (min … max):  4.818 ms …  16.852 ms  ┊ GC (min … max): 0.00% … 69.74%
 Time  (median):     5.380 ms               ┊ GC (median):    7.12%
 Time  (mean ± σ):   5.452 ms ± 690.859 μs  ┊ GC (mean ± σ):  6.57% ±  5.49%

              ▁▄█▅▃▄▅▄ ▂▂▄▂▁                                   
  ▄▅▄▆▅▅▆▆▆▆▄▆██████████████▇█▆▇▄▅▅▆▃▄▂▃▃▃▂▂▄▂▂▂▃▂▁▂▁▂▂▁▁▂▁▁▂ ▄
  4.82 ms         Histogram: frequency by time        6.59 ms <

 Memory estimate: 17.28 MiB, allocs estimate: 38533.

確率的動学モデル

この節では確率的な不確実性を導入する. 生産関数 \(f(k, z) = z k^\alpha\) とし, 生産性 \(z\) が確率的に変化するモデルを考えます. 生産性 \(z\) は以下のような AR(1) 過程に従うとします.

\[ \log z' = (1 - \rho) \mu + \rho \log z + \sigma \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, 1). \]

この時, ベルマン方程式は以下のように書き換えられます.

\[ V(k, z) = \max_{k'} \left\{u(f(k, z) + (1-\delta)k - k') + \beta \mathbb{E}[V(k', z')| z]\right\}. \]

この \(V(k, z)\) の数値計算するためには, \(z\) の状態空間を離散化する必要があります. ここでは, Tauchen Method (Tauchen 1986) を用いて \(z\) の状態空間を離散化します.

Tauchen Method

実数列 \(x_t\) が以下のAR(1)過程に従うとします.

\[ x_{t+1} = (1 - \rho) \mu + \rho x_t + \sigma \varepsilon_{t+1}, \quad \varepsilon_{t+1} \sim \mathcal{N}(0, 1). \]

これを \(n\) 個のグリッドポイント \(\{x_1, \dots, x_n\}\) におけるマルコフ過程として離散化することを考えます. すなわち, 以下のような遷移確率行列 \(\Lambda\) を考えます.

\[ \begin{pmatrix} x_{1, t+1} \\ \vdots \\ x_{n, t+1} \end{pmatrix} = \begin{pmatrix} \lambda_{1, 1} & \cdots & \lambda_{1, n} \\ \vdots & \ddots & \vdots \\ \lambda_{n, 1} & \cdots & \lambda_{n, n} \end{pmatrix} \begin{pmatrix} x_{1, t} \\ \vdots \\ x_{n, t} \end{pmatrix} \]

ここで, \(\lambda_{i, j} = \Pr(x_{t+1} = x_j \mid x_t = x_i)\) であり, \(\sum_{j=1}^{n} \lambda_{i, j} = 1\) が成り立ちます. この遷移確率行列を求める方法として, Tauchen Method が知られています.

Tauchen Method

\(x\) の定義域を \([x_{\min}, x_{\max}]\) を \(n\) 個のグリッドポイントに分割する. 通常は \(x\) が従う分布の \(3\sigma\) 程度の範囲を考える.
グリッドポイント \(x_i\) から \(x_j\) に遷移する確率を以下のように定義する.

\[ \begin{aligned} \lambda_{i, 1} &= \Phi\left(\frac{x_1 + \frac{d}{2} - (1 - \rho) \mu - \rho x_i}{\sigma}\right) \\ \lambda_{i, j} &= \Phi\left(\frac{x_j + \frac{d}{2} - (1 - \rho) \mu - \rho x_i}{\sigma}\right) - \Phi\left(\frac{x_{j} - \frac{d}{2} - (1 - \rho) \mu - \rho x_i}{\sigma}\right) & (2 \leq j \leq n-1) \\ \lambda_{i, n} &= 1 - \Phi\left(\frac{x_{n} - \frac{d}{2} - (1 - \rho) \mu - \rho x_i}{\sigma}\right) \end{aligned} \]

Figure 6: **Visualization of Tauchen Method**. \(\mu = 0\).

function tauchen_method(; n_std=3.0, n=5, ρ=0.9, μ=0.1, σ=1.0)

    x_min = μ - n_std * sqrt(σ^2 / (1 - ρ^2))
    x_max = μ + n_std * sqrt(σ^2 / (1 - ρ^2))
    x = range(x_min, x_max, length=n)

    d = (x_max - x_min) / (n - 1)
    Λ = zeros(n, n)
    for i in 1:n
        Λ[i, 1] = cdf(Normal(0.0, σ), x[1] + d / 2 - (1 - ρ) * μ - ρ * x[i])
        Λ[i, n] = 1 - cdf(Normal(0.0, σ), x[n] - d / 2 - (1 - ρ) * μ - ρ * x[i])
        for j in 2:n-1
            Λ[i, j] = cdf(Normal(0.0, σ), x[j] + d / 2 - (1 - ρ) * μ - ρ * x[i]) -
                      cdf(Normal(0.0, σ), x[j] - d / 2 - (1 - ρ) * μ - ρ * x[i])
        end
    end

    return x, Λ
end


x, Λ = tauchen_method()
@show x
Λ

x = -6.782472016116854:3.4412360080584268:6.982472016116853

5×5 Matrix{Float64}:
 0.849051     0.150945     3.84556e-6  1.22125e-15  0.0
 0.0194737    0.896192     0.0843336   7.26002e-7   1.11022e-16
 1.22258e-7   0.04266      0.91468     0.04266      1.22258e-7
 7.34696e-17  7.26002e-7   0.0843336   0.896192     0.0194737
 3.45903e-30  1.23783e-15  3.84556e-6  0.150945     0.849051

実用上は QuantEcon.jl の tauchen 関数を利用するのがいいでしょう.

using QuantEcon

mc = tauchen(5, 0.9, 1.0, 0.1) # tauchen(N, ρ, σ, μ)
@show mc.state_values
mc.p

mc.state_values = -5.8824720161168536:3.4412360080584268:7.8824720161168536

5×5 Matrix{Float64}:
 0.849051     0.150945     3.84556e-6  1.22125e-15  0.0
 0.0194737    0.896192     0.0843336   7.26002e-7   1.11022e-16
 1.22258e-7   0.04266      0.91468     0.04266      1.22258e-7
 7.34696e-17  7.26002e-7   0.0843336   0.896192     0.0194737
 3.45903e-30  1.23783e-15  3.84556e-6  0.150945     0.849051

ただし, ここでは AR(1) process が以下のように定義されていることに注意してください.

\[ y_{t+1} = \mu + \rho y_{t} + \varepsilon_{t+1}, \quad \varepsilon_{t+1} \sim \mathcal{N}(0, \sigma^2). \]

そのため, \(\mathbb{E}[y_{t}] = \frac{\mu}{1-\rho}\) となります. 上の例では, \(\mu = 0.1\), \(\rho = 0.9\), \(\sigma = 1.0\) としているため, \(\mathbb{E}[y_{t}] = 1.0\) となります.

ps = vcat(stationary_distributions(mc)...)
ps' * collect(mc.state_values)

1.0

AR(1) processの期待値を任意の値 \(\tilde{\mu}\) に設定するためには, \(\mu = (1-\rho)\tilde{\mu}\) とすればよいです.

mc = tauchen(5, 0.9, 1.0, 0.01)
@show mc.state_values
mc.p

mc.state_values = -6.782472016116854:3.4412360080584268:6.982472016116853

5×5 Matrix{Float64}:
 0.849051     0.150945     3.84556e-6  1.22125e-15  0.0
 0.0194737    0.896192     0.0843336   7.26002e-7   1.11022e-16
 1.22258e-7   0.04266      0.91468     0.04266      1.22258e-7
 7.34696e-17  7.26002e-7   0.0843336   0.896192     0.0194737
 3.45903e-30  1.23783e-15  3.84556e-6  0.150945     0.849051

Value Function Iteration

AR(1) process を \(\Lambda\) で離散化した後のベルマン方程式は以下のようになります.

\[ V(k, z) = \max_{k'} \left\{u(f(k, z) + (1-\delta)k - k') + \beta \sum_{z'} \Lambda(z, z') V(k', z')\right\}. \]

VFIを用いる場合, 実装は価値関数 \(V\) の次元を \(z\) のために1つ増やし, ベルマン方程式の通り実装するだけです.

@kwdef struct StochasticModel{TF<:AbstractFloat,TI<:Integer}
    # Parameters
    α::TF = 0.4
    β::TF = 0.96
    δ::TF = 0.1
    # AR(1) process
    ρ::TF = 0.6
    μ::TF = 0.0
    σ::TF = 0.4

    # Value Function
    n_k::TI = 101
    n_z::TI = 5
    k_min::TF = 0.05
    k_max::TF = 0.5
    k_grid::Vector{TF} = collect(range(k_min, k_max, length=n_k))
    mc::MarkovChain = tauchen(n_z, ρ, σ, μ)
    z_grid::Vector{TF} = collect(exp.(mc.state_values))
    Λ::Matrix{TF} = mc.p
    V::Matrix{TF} = zeros(n_k, n_z)
end

f(k, z; m) = z * k^m.α
function vfi!(m::StochasticModel; tol=1e-6, max_iter=1000)
    (; k_grid, z_grid, β, δ, Λ) = m

    iter, dist = 0, Inf
    V_new = similar(m.V)
    while dist > tol && iter < max_iter
        for (i_z, z) in enumerate(z_grid), (i_k, k) in enumerate(k_grid)
            c_min = 1e-9 # consumption cannot be negative
            V_new[i_k, i_z] = maximum(
                log(max(f(k, z; m) + (1 - δ) * k - k′, c_min)) +
                β * sum(Λ[i_z, i_z′] * m.V[i_k′, i_z′]
                        for (i_z′, z′) in enumerate(z_grid))
                for (i_k′, k′) in enumerate(k_grid)
            )
        end

        dist = maximum(abs, V_new .- m.V)
        m.V .= V_new
        iter += 1
    end

    if iter == max_iter
        println("Warning: Maximum iterations reached.")
    else
        println("Converged in $iter iterations.")
    end

    return nothing
end


m = StochasticModel()
vfi!(m)
p = plot(m.k_grid, m.V[:, 5], label=L"z_5", xlabel=L"Capital $k$", ylabel=L"Value $V(k, z)$")
for i in 4:-1:1
    plot!(m.k_grid, m.V[:, i], label=L"z_%$i")
end

p

Converged in 316 iterations.

Figure 7: Value Function of Stochastic Model.

References

Tauchen, George. 1986. “Finite State Markov-Chain Approximations to Univariate and Vector Autoregressions.” Economics Letters 20 (2): 177–81. https://doi.org/10.1016/0165-1765(86)90168-0.

北尾早霧, 砂川武貴, and 山田知明. 2024. 定量的マクロ経済学と数値計算. 日本評論社.

Footnotes

線形補間では一階微分までしかとれませんが, スプライン補間を用いれば二階微分が取れるので, derivative-based な最適化手法も使えます. ここでは1変数関数なので, bracketing method のような微分のいらない手法を用いています.↩︎