---
engine: julia
---
# Foundations
```{julia}
#| label: setup-proof
#| output: false
#| code-fold: true
using Plots
using Distributions
using Roots
using LaTeXStrings
import Random
Random.seed!(1234)
mypalette = ["#d5968c", "#c2676d", "#5c363a", "#995041", "#45939c", "#0f6a81"]
default(size=(500, 309), titlefontsize=10, fmt=:svg, grid=false,
foreground_color_legend=nothing, palette=mypalette[[6, 2, 4, 3, 1, 5]])
```
## AR(1) Process {#sec-ar1-process}
::: {.callout-note}
## AR(1) 過程
AR(1) 過程 $\{x_{t}\}_{t=0}^{\infty}$ は以下のように表される.
$$
x_{t+1} = \bar{x} + \rho x_t + \varepsilon_{t}, \quad \varepsilon_t \overset{\text{iid}}{\sim} (0, \sigma^2).
$$
:::
平均 $0$ で分散 $\sigma^2$ (有限) の i.i.d. な確率変数 ($\varepsilon_t$) は特にホワイトノイズ (white noise) と言います. 以下の内容は, ホワイトノイズであれば同様に成り立ちますが, 実務上は正規分布 $\varepsilon_t \sim \mathcal{N}(0, \sigma^2)$ を仮定することが多いです.
```{julia}
#| label: fig-ar1-process
#| fig-cap: AR(1) 過程のシミュレーション. $\bar{x} = 0, \rho = 0.9, \sigma = 0.01$.
#| code-fold: true
x̄ = 0.0
ρ = 0.9
σ = 0.01
T = 100
xs = fill(x̄, T)
for t in 2:T
xs[t] = x̄ + ρ * xs[t-1] + rand(Normal(0, σ))
end
plot(1:T, xs, label=false, size=(500, 309))
```
::: {#prp-mean-var-ar1}
## AR(1) 過程の平均と分散
$$
\mathbb{E}[x_t] = \frac{\bar{x}}{1-\rho}, \quad \text{Var}(x_t) = \frac{\sigma^2}{1-\rho^2}.
$$
:::
_証明_
$$
\begin{aligned}
\mathbb{E}[x_{t+1}] &= \mathbb{E}[\bar{x} + \rho x_t + \varepsilon_{t}]= \bar{x} + \rho \mathbb{E}[x_t] + 0\\
\text{Var}(x_{t+1}) &= \text{Var}(\bar{x} + \rho x_t + \varepsilon_{t}) = 0 + \rho^2 \text{Var}(x_t) + \sigma^2.
\end{aligned}
$$
$\mathbb{E}[x_{t+1}] = \mathbb{E}[x_t]$ かつ, $\text{Var}(x_{t+1}) = \text{Var}(x_t)$ より, 証明された.
::: {#prp-ac-ar1}
## AR(1) 過程の自己共分散 (Autocovariance)
$$
\text{Cov}(x_{t}, x_{t+k}) = \frac{\rho^{k}}{1-\rho^2}\sigma^2.
$$
:::
_証明_
$k = 1$ の時,
$$
\begin{aligned}
\text{Cov}(x_{t}, x_{t+1})
&= \text{Cov}(x_{t}, \bar{x} + \rho x_t + \varepsilon_{t}) \\
&= \text{Cov}\rho(x_{t}, x_t) + \text{Cov}(x_{t}, \varepsilon_{t})\\
&= \rho \text{Var}(x_t) + 0 \\
&= \frac{\rho}{1-\rho^2}\sigma^2.
\end{aligned}
$$
$k = l + 1$ の時, $\text{Cov}(x_{t}, x_{t+l})= \frac{\rho^{l}}{1-\rho^2}\sigma^2$ として,
$$
\begin{aligned}
\text{Cov}(x_{t}, x_{t+l+1})
&= \text{Cov}(x_{t}, \bar{x} + \rho x_{t+l} + \varepsilon_{t+l}) \\
&= \text{Cov}(x_{t}, \bar{x}) + \text{Cov}(x_{t}, \rho x_{t+l}) + \text{Cov}(x_{t}, \varepsilon_{t+l}) \\
&= 0 + \rho \text{Cov}(x_{t}, x_{t+l}) + 0 \\
&= \frac{\rho^{l+1}}{1-\rho^2}\sigma^2.
\end{aligned}
$$
## Labor Supply
### Static Model
まずは静学的な労働供給モデルを考えます. 家計は以下の効用最大化問題を解きます:
$$
\max_{c, h} u(c, h) \quad \text{ s.t. } c \leq w h + I.
$${#eq-labor-static}
ここで, 消費 $c \geq 0$, 労働供給 $h \geq 0$, 賃金率 $w > 0$, 非労働所得 $I \geq 0$ です. また, 賃金率 $w$ と非労働所得 $I$ に対する労働供給を $h(w, I)$ とします. この時, 労働供給の弾力性は以下で定義されます:
$$
\begin{aligned}
\eta_{h, w}^M &= \frac{\partial h(w, I)}{\partial w} \frac{w}{h(w, I)} \\
\eta_{h, I} &= \frac{\partial h(w, I)}{\partial I} \frac{I}{h(w, I)}
\end{aligned}
$$
特に, 所得を一定とした時の賃金弾力性 $\eta_{h, w}^M$ をマーシャリアン弾力性 (Marshallian elasticity) と呼びます. 次に, 労働供給の代替効果と所得効果を分離するために, ヒックス弾力性 (Hicksian elasticity) を導入します. 家計は以下の支出最小化問題を解きます:
$$
\min_{c, h} c - wh \quad \text{ s.t. } u(c, h) \geq \bar{u}.
$$
ここで, 効用水準 $\bar{u}$ に対する労働供給を $h^H(w, \bar{u})$ とします. この時, ヒックス弾力性は以下で定義されます:
$$
\eta_{h, w}^H = \frac{\partial h^H(w, \bar{u})}{\partial w} \frac{w}{h^H(w, \bar{u})}.
$$
::: {.callout-note}
## スルツキーの恒等式 (Slutsky Equation)
所得 $I$, 価格ベクトル $\mathbf{p}$ のもとで次の家計の効用最大化問題を考える.
$$
\max_{\mathbf{x}} u(\mathbf{x}) \quad \text{ s.t. } \mathbf{p}^{\top} \mathbf{x} \leq I.
$$
この時, 財 $i$ の最適消費を $x_i(\mathbf{p}, I)$ とする. また, 効用水準 $\bar{u}$ に対する支出最小化問題を考える.
$$
\min_{\mathbf{x}} \mathbf{p}^{\top} \mathbf{x} \quad \text{ s.t. } u(\mathbf{x}) \geq \bar{u}.
$$
この時, 財 $i$ の最適消費を $x_i^H(\mathbf{p}, \bar{u})$ とし, 支出関数を $I(\mathbf{p}, \bar{u})$ とする. 連鎖律から次の式を得る.
$$
\frac{\partial x_i^H(\mathbf{p}, \bar{u})}{\partial p_i} = \frac{\partial x_i(\mathbf{p}, I(\mathbf{p}, \bar{u}))}{\partial p_i} + \frac{\partial x_i(\mathbf{p}, I(\mathbf{p}, \bar{u}))}{\partial I(\mathbf{p}, \bar{u})} \frac{\partial I(\mathbf{p}, \bar{u})}{\partial p_i}
$$
また, シェファードの補題 (Shephard's lemma) $\frac{\partial I(\mathbf{p}, \bar{u})}{\partial p_i} = x_i^H(\mathbf{p}, \bar{u})$ より,
$$
\frac{\partial x_i}{\partial p_i} = \frac{\partial x_i^H(\mathbf{p}, \bar{u})}{\partial p_i} - \frac{\partial x_i(\mathbf{p}, I)}{\partial I} x_i^H(\mathbf{p}, \bar{u}).
$$
これをスルツキーの恒等式 (Slutsky equation) と呼ぶ.
:::
労働供給の問題は価格1の財と価格 $-w$ の労働の2財問題とみなせるため,
$$
\frac{\partial h(w, I)}{\partial w} = \frac{\partial h^H(w, \bar{u})}{\partial w} + \frac{\partial h(w, I)}{\partial I} h^H(w, \bar{u}).
$$
ここで, 価格が負であるためシェファードの補題が $\frac{\partial I(w, \bar{u})}{\partial w} = -h^H(w, \bar{u})$ となることに注意してください. 全体を $\frac{w}{h}$ で掛けると^[
均衡において $h(w, I) = h^H(w, \overline{u})$ であるため, 右辺の $h(w, I)$ を $h^H(w, \overline{u})$ に置き換えても問題ありません.
],
$$
\underbrace{\frac{\partial h(w, I)}{\partial w} \frac{w}{h}}_{\eta_{h, w}^M} = \underbrace{\frac{\partial h^H(w, \bar{u})}{\partial w} \frac{w}{h}}_{\eta_{h, w}^H} + \underbrace{\frac{\partial h(w, I)}{\partial I} \frac{I}{h}}_{\eta_{h, I}} \frac{wh}{I}.
$$
この式は, 労働供給の弾力性が次の代替効果と所得効果の和に分解されることを示しています.
- **代替効果**: 賃金の増加によって労働のリターンが増加し, 労働供給を増やす効果. 余暇の価格 (労働の機会費用) の増加に対する財への代替と考えてもよい
- **所得効果**: 所得の増加によって働く必要がなくなり, 労働供給を減らす効果. 余暇が正常財である場合, 所得の増加によって余暇の消費が増加し, 労働供給が減少すると考えてもよい
::: {#exm-labor-loglog}
## Cobb-Douglas Utility
$$
\max_{c, l} (1-\phi) \log c + \phi \log l \quad \text{ s.t. } c \leq w(1-l) + I.
$$
$$
\eta_{h, w}^M = \frac{\phi}{w}\frac{I}{h},
$$
- 代替効果: $\eta_{h, w}^H = \frac{\phi}{w}\frac{I}{h} + \phi$
- 所得効果: $\eta_{h, I}\frac{wh}{I}=-\phi$
$I = 0$ の時, 代替効果 $\phi$ と所得効果 $\phi$ が打ち消しあい, 労働供給の賃金弾力性は0となる.
:::
::: {.proof}
最適条件より,
$$
h(w, I) = 1 - \phi\left(1 + \frac{I}{w}\right).
$$
よって,
$$
\begin{aligned}
\eta_{h, w}^M &= \frac{\partial h(w, I)}{\partial w} \frac{w}{h(w, I)} = \frac{\phi I}{w^2} \frac{w}{h(w, I)} = \frac{\phi}{w}\frac{I}{h(w, I)}\\
\eta_{h, I} &= \frac{\partial h(w, I)}{\partial I} \frac{I}{h(w, I)} = -\frac{\phi}{w} \frac{I}{h(w, I)}.
\end{aligned}
$$
スルツキーの恒等式より,
$$
\eta_{h, w}^H = \eta_{h, w}^M - \eta_{h, I} \frac{wh(w, I)}{I} = \frac{\phi}{w}\frac{I}{h(w, I)} + \phi.
$$
**費用最小化問題**
スルツキーの恒等式を確認するために, ヒックス弾力性を定義から計算する.
$$
\min_{c, l} c - w(1-l) \quad \text{ s.t. } (1-\phi) \log c + \phi \log l \geq \bar{u}.
$$
ラグランジアン $\mathcal{L} = c - w(1-l) + \lambda \left(\bar{u} - (1-\phi) \log c - \phi \log l \right)$ の最適条件は次の通り.
$$
\begin{aligned}
\frac{\partial \mathcal{L}}{\partial c} &= 1 - \lambda \frac{1-\phi}{c} = 0, \\
\frac{\partial \mathcal{L}}{\partial l} &= w - \lambda \frac{\phi}{l} = 0, \\
\frac{\partial \mathcal{L}}{\partial \lambda} &= \bar{u} - (1-\phi) \log c - \phi \log l = 0.
\end{aligned}
$$
一階条件から, $c = \frac{1-\phi}{\phi}wl$ より, 効用制約に代入すると,
$$
\bar{u} = (1-\phi) \log \left(\frac{1-\phi}{\phi}wl\right) + \phi \log l = (1-\phi) \log \left(\frac{1-\phi}{\phi}w\right) + \log l.
$$
よって, 関数 $A(\bar{u}; \phi)$ を用いて, $l = A(\bar{u}; \phi) w^{-(1-\phi)}$ と表せるので,
$$
h^H(w, \bar{u}) = 1 - A(\bar{u}; \phi) w^{-(1-\phi)}.
$$
ヒックシアン弾力性は定義から,
$$
\begin{aligned}
\eta_{h, w}^H &= \frac{\partial h^H(w, \bar{u})}{\partial w} \frac{w}{h^H(w, \bar{u})}\\
&= (1-\phi) A(\bar{u}; \phi) w^{-(1-\phi)-1} \frac{w}{h^H(w, \bar{u})} \\
&= (1-\phi) \frac{1 - h^H(w, \bar{u})}{h^H(w, \bar{u})} & \because 1 - h^H(w, \bar{u}) = A(\bar{u}; \phi) w^{-(1-\phi)} \\
\end{aligned}
$$
ここで一階条件から $1 - h^H(w, \bar{u}) = \frac{\phi}{1-\phi}\frac{c}{w}$ であり, 均衡における予算制約 $c = w h^H(w, \bar{u}) + I$ を用いると,
$$
\begin{aligned}
\eta_{h, w}^H &= \frac{\phi}{w} \frac{c}{h^H(w, \bar{u})} \\
&= \frac{\phi}{w} \frac{w h^H(w, \bar{u}) + I}{h^H(w, \bar{u})} \\
&= \frac{\phi}{w} \frac{I}{h^H(w, \bar{u})} + \phi.
\end{aligned}
$$
:::
### Dynamic Model
次の動学的労働供給モデルを考えます:
$$
V(a) = \max_{c, h, a'} u(c, h) + \beta V(a') \quad \text{s.t. } c + a' = (1 + r)a + wh.
$$
ここで, $I := (1+r)a - a'$ と置くと, 静学的な労働供給モデル (-@eq-labor-static) と同じ形になります. 動学的な労働供給モデルにおける賃金弾力性には, 主に次の3種類があります.
1. Marshallian 弾力性: 所得を一定とした下での労働供給の弾力性
2. Hicksian 弾力性: 効用のレベルを一定とした下での労働供給の弾力性
3. Frisch 弾力性: 限界効用を一定にした下での労働供給の弾力性
**Marshallian & Hicksian Elasticities**
Marshallian 弾力性は, 静学的な労働供給と同様に定義することが可能です. 直感的には, 賃金が下落した際に貯金を切り崩すといった埋め合わせ (compensation) 行動は可能ですが, その点を考慮せずに, 所得 $I$ を一定とした下での労働供給の弾力性として定義されます. その意味でMarshallian 弾力性は uncompesated elasticity とも呼ばれます.
Hicksian 弾力性も静学的な労働供給と同様に定義することが可能です. 賃金が下落した際に貯金を切り崩すことで効用水準を一定に保つ状況を考えるという意味で, Hicksian 弾力性は compensated elasticity とも呼ばれます.
::: {#exm-labor-crra}
## CRRA Utility
$$
u(c, h) = \frac{c^{1-\sigma}-1}{1-\sigma} - \psi \frac{h^{1+\frac{1}{\gamma}}}{1+\frac{1}{\gamma}}.
$$
$$
\eta_{h, w}^M = \frac{1 - \sigma(1-s)}{\frac{1}{\gamma} + \sigma(1-s)}.
$$
- 代替効果: $\eta_{h, w}^H = \frac{1}{\frac{1}{\gamma} + \sigma(1-s)}$
- 所得効果: $\eta_{h, I}\frac{wh}{I} = -\frac{\sigma - \sigma s}{\frac{1}{\gamma} + \sigma(1-s)}$
ここで $s = \frac{I}{c}$ とおいた.
:::
::: {.proof}
Intra-temporal な最適条件より,
$$
c^{-\sigma} w = \psi h^{\frac{1}{\gamma}}.
$${#eq-labor-crra-intra}
対数をとって, $\log w$ で微分すると,
$$
-\sigma \frac{\partial \log c}{\partial \log w} + 1 = \frac{1}{\gamma} \frac{\partial \log h}{\partial \log w}.
$${#eq-labor-crra-intra-wage}
予算制約 $\log c = \log (w h + I)$ を $\log w$ で微分すると,
$$
\begin{aligned}
\frac{\partial \log c}{\partial \log w} &= \frac{1}{wh + I}\frac{\partial (wh+I)}{\partial w}\frac{\partial w}{\partial \log w} \\
&= \frac{wh}{wh + I}\left(1 + \frac{\partial \log h}{\partial \log w}\right)\\
&= (1 - s) \left(1 + \frac{\partial \log h}{\partial \log w}\right).
\end{aligned}
$$
これを @eq-labor-crra-intra-wage に代入すると,
$$
\eta_{h, w}^M = \frac{\partial \log h}{\partial \log w} = \frac{1 - \sigma(1-s)}{\frac{1}{\gamma} + \sigma(1-s)}.
$$
@eq-labor-crra-intra の対数を取って, $\log I$ で微分すると,
$$
-\sigma \frac{\partial \log c}{\partial \log I} = \frac{1}{\gamma} \frac{\partial \log h}{\partial \log I}.
$${#eq-labor-crra-intra-income}
予算制約 $\log c = \log (w h + I)$ を $\log I$ で微分すると,
$$
\begin{aligned}
\frac{\partial \log c}{\partial \log I} &= \frac{1}{wh + I}\frac{\partial (wh+I)}{\partial \log I} \\
&= \frac{1}{wh + I}\left(wh\frac{\partial \log h}{\partial \log I} + I\right)\\
&= \frac{1}{c}\left((c-I)\frac{\partial \log h}{\partial \log I}+ I\right) \\
&= (1 - s) \frac{\partial \log h}{\partial \log I}+ s.
\end{aligned}
$$
これを @eq-labor-crra-intra-income に代入すると,
$$
\eta_{h, I} = \frac{\partial \log h}{\partial \log I} = -\frac{\sigma s}{\frac{1}{\gamma} + \sigma(1-s)}.
$$
スルツキーの恒等式より,
$$
\begin{aligned}
\eta_{h, w}^H &= \eta_{h, w}^M - \eta_{h, I} \frac{wh}{I} \\
&= \frac{1 - \sigma(1-s)}{\frac{1}{\gamma} + \sigma(1-s)} + \frac{\sigma s}{\frac{1}{\gamma} + \sigma(1-s)} \left(\frac{1}{s} - 1\right) \\
&=\frac{1}{\frac{1}{\gamma} + \sigma(1-s)}.
\end{aligned}
$$
:::
**Frisch Elasticity**
Frisch 弾力性は, 消費の限界効用を一定にした下での労働供給の弾力性として定義されます. また, しばしば lifetime wealth に対する限界効用を一定にした下での労働供給の弾力性としても定義されます. これは恒常所得仮説 (permanent income hypothesis) に基づくと, 消費の限界効用が lifetime wealth のみに依存すると考えるためです.
直感的には, 一時的な賃金の上昇は lifetime wealth にほとんど影響を与えないため, Frisch 弾力性は一時的な賃金変動に対する労働供給の反応を捉えると考えられます. そのため, RBCモデルなどで確率過程に基づく賃金変動が考慮される場合, Frisch 弾力性が重要な役割を果たします.
一階条件から, ラグランジュ乗数は $\lambda = u_c(c, h)$ であり, Frish 弾力性はこの $\lambda$ を一定にした下での労働供給の弾力性として定義されます. すなわち, $u_h(c, h) = -\lambda w$ より,
$$
u_{cc} \frac{\partial c}{\partial w} + u_{ch} \frac{\partial h}{\partial w} = 0 \text{ and }
u_{hc} \frac{\partial c}{\partial w} + u_{hh} \frac{\partial h}{\partial w} = -\lambda.
$$
これを $\frac{\partial h}{\partial w}$ について解くと,
$$
\frac{\partial h}{\partial w} = \frac{u_{h}}{u_{hh} - \frac{u_{ch}^2}{u_{cc}}} \frac{1}{w}.
$$
Frisch 弾力性 $\eta_{h, w}^F := \frac{\partial h}{\partial w}\frac{w}{h}|_{\lambda}$ より,
$$
\eta_{h, w}^F = \frac{u_{h}}{h\left(u_{hh} - \frac{u_{ch}^2}{u_{cc}}\right)}.
$$
なお, この式から自明に @exm-labor-crra において $\eta_{h, w}^F = \gamma$ となります. @chetty2011 はメタ分析を行い, Frisch 弾力性の平均的な推定値を報告しました.
- **ミクロ分析** (擬似実験, quasi-experiment): 0.82
- **マクロモデル**: 2.84
マクロモデルの結果は, RBCモデルなどで雇用は実質賃金よりも volatile であることを反映しています. この「マクロモデルの Frisch 弾力性はミクロ分析よりもかなり大きい」というPuzzleに関しては今も議論が続いています (e.g., @erosa2016).
## Root Finding
定常状態の計算や, 一階条件から解を求める際, 数値的に非線形方程式の根を求める必要があります. ここではその基礎である一変数の非線形方程式の根の探索方法を紹介します.
例えば以下の最適化問題を考えてみましょう.
$$
\max_{c, l} \frac{c^{1-\gamma_c}}{1-\gamma_c} + \alpha_l \frac{l^{1-\gamma_l}}{1-\gamma_l} \quad \text{s.t.} \quad c = w (1-l)
$${#eq-max-utility}
この問題は以下の一階条件の解を求めることで解くことができます.
$$
w^{1-\gamma_c} (1-l)^{-\gamma_c} - \alpha_l l^{-\gamma_l} = 0.
$$
```{julia}
#| label: fn-foc-crra
#| output: false
foc(l; w=1.0, γ_c=1.5, α_l=1.2, γ_l=1.5) =
w^(1 - γ_c) * (1 - l)^(-γ_c) - α_l * l^(-γ_l)
```
```{julia}
#| code-fold: true
#| label: fig-foc-crra
#| fig-cap: FOC for CRRA utility
plot(0.1:0.01:0.9, foc, label=L"w^{1-\gamma_c} (1-l)^{-\gamma_c} - \alpha_l l^{\gamma_l}", lw=2)
hline!([0.0], ls=:dash, lw=2, label=false)
```
これは, 一般に解析解がありません ($\gamma_c = \gamma_l = 1$などを除く.) そのため, 数値計算を用いて解く必要があります. この様な一変数の非線形方程式 $f(x) = 0$ の解法として, 以下の2つの選択肢があります.
1. **Non-bracketing method**: 初期値 $x_0$ からスタートして, $f(x_n)$ が十分0に近づくまで反復的に計算. (ニュートン法など)
1. **Bracketing method**: 区間 $[a, b]$ を選び, 区間を狭めていくことで解を求める. (二分法など)
一般に, Non-bracketing method は収束が速い代わりにその収束は保証されません. 一方で, Bracketing method は収束がわずかに遅い代わりに, 解が存在する区間を保証することができます. 以下では, その理由とそれぞれの適した利用場面を解説します.
### Non-bracketing Methods
Non-bracketing method は基本的には一階微分を利用し, 初期値 $x_0$ からスタートして反復的に解を求めます.
ここでは, 最も古典的で簡単なニュートン法を紹介します.
::: {.callout-tip}
## ニュートン法
1. 初期値 $x_0$ を選ぶ. また停止条件として十分小さい $\epsilon > 0$ を選ぶ.
1. $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$ を計算する.
1. $|f(x_{n+1})| < \epsilon$ ならば停止し, $x_{n+1}$ を解として返す. そうでなければ, $n \leftarrow n+1$ として2に戻る.
:::
{#fig-newton-method width=80%}
```{julia}
#| label: find-zero-newton
dfoc(l; w=1.0, γ_c=1.5, α_l=1.2, γ_l=1.5) =
w^(1 - γ_c) * γ_c * (1 - l)^(-γ_c - 1) + α_l * γ_l * l^(-γ_l - 1)
find_zero((foc, dfoc), 0.5, Roots.Newton())
```
**Hybrid method**
Juliaの`Roots.jl`のデフォルトのNon-bracketing method (`Order0()`)は, 厳密にはBracketing methodも一部利用します. また一階微分も必要としないため, より汎用的に利用することができます. ただし, 収束が速く収束が保証されない性質はNon-bracketing methodと同じです.
```{julia}
#| label: find-zero-hybrid
l = find_zero(foc, 0.5)
```
**定義域の変換**
Non-bracketing method では, $x_{n+1}$ の値を更新する際に, 必ずしも $x_{n+1}$ が定義域内にあることが保証されません. 例えば, 今回の例では 余暇時間ですので, $l \in (0, 1)$ ですが, 初期値や関数形によっては, $l < 0$ や $l > 1$ となる可能性があります.
この場合, 端点解が存在するかどうかを確認する必要があります. 端点解が想定される場合は, 次節で説明するように Bracketing method を利用することが適しています. 今回の例では, 端点解が存在しないため ($l = 0, 1$ のどちらにおいても効用が負の無限大に発散するため), Non-bracketing method での解の収束が保証されています. 以下の様な定義域の変換を用いることで, 解の収束を保証することができます.
::: {.callout-note}
## 定義域の変換
$x$ が区間 $(a, b)$ で定義されるとき, 以下の変換を用いて定義域を $(-\infty, \infty)$ にすることができる.
$$
x = \frac{a + b \exp y}{1 + \exp y}
$$
明らかに, $y \rightarrow -\infty$ のとき $x \rightarrow a$, $y \rightarrow \infty$ のとき $x \rightarrow b$ となる.
:::
今回の例では, $l = \frac{1}{1 + \exp(y)}$ (つまり[シグモイド関数](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%82%B0%E3%83%A2%E3%82%A4%E3%83%89%E9%96%A2%E6%95%B0)) という変換を用いることができます.
```{julia}
#| label: find-zero-transform
y = find_zero(y -> foc(1 / (1 + exp(y))), 0.0)
l = 1 / (1 + exp(y))
```
### Bracketing Methods
Bracketing method は, 解が存在する区間 $[a, b]$ を選び, その区間を狭めていくことで解を求めます. ここでは最も基本的な二分法を用いることで解を求めることができます.
::: {.callout-tip}
## 二分法
1. 区間 $[a, b]$ を選ぶ. また停止条件として十分小さい $\epsilon > 0$ を選ぶ.
1. $f(a) \cdot f(b) < 0$ ならば, $c = \frac{a + b}{2}$ を計算する.
1. $|f(c)| < \epsilon$ ならば, $c$ を解として返す. そうでなければ, $f(a) \cdot f(c) < 0$ ならば, $b \leftarrow c$ として2に戻る. そうでなければ, $a \leftarrow c$ として2に戻る.
:::
{#fig-bisection-method width=80%}
例えば, 今回の例では $l \in (0, 1)$ ですので, $l = 0$ と $l = 1$ で効用が発散することから, 解が存在する区間は $(0, 1)$ です. なお, 端点でFOCが定義されてない場合も多いので, 端点より微小に内側の区間を選ぶことが多いです.
```{julia}
#| label: find-zero-bracketing
ϵ = 1e-9
l = find_zero(foc, (ϵ, 1 - ϵ))
```
**端点解の存在**
以下の様なカップルの意思決定問題を考えてみましょう.
$$
\max_{c_m, c_f, l_m, l_f} u(c_m, l_m) + u(c_f, l_f)
$${#eq-couple-problem}
subject to
$$
c_m + c_f = w_m (1 - l_m) + w_f (1 - l_f).
$$
一階条件から以下の関係が導けます:
$$
\frac{l_m}{l_f} = \left(\frac{w_f}{w_m}\right)^{\frac{1}{\gamma_l}}.
$$
この関係から一方の相対賃金が十分に大きい場合, 例えば $w_m \gg w_f$ の場合, $l_f > 1$ となる可能性があります ($l_m$ は余暇時間のため, 0よりある程度大きい値をとるということが想像できます.) 内点解が存在する場合の $l_f$ に関する一階条件は以下のようになります.
$$
w_m\left(1-\left(\frac{w_f}{w_m}\right)^{\frac{1}{\gamma_l}}l_f\right) + w_f(1-l_f) - 2\left(\frac{w_f}{\alpha_l}\right)^{\frac{1}{\gamma_c}}l_f^{\frac{\gamma_l}{\gamma_c}} = 0.
$$
```{julia}
#| code-fold: true
#| label: fig-foc-bargaining
#| fig-cap: "**FOC for bargaining problem.** $w_m = 1.0$, $\\alpha_l = 1.2$, $\\gamma_c = \\gamma_l = 1.5$."
function foc_bargaining(l_f; w_m=1.0, w_f=0.5, γ_c=1.5, α_l=1.2, γ_l=1.5)
return w_m * (1 - (w_f / w_m)^(1 / γ_l) * l_f) + w_f * (1 - l_f) -
2 * (w_f / α_l)^(1 / γ_c) * l_f^(γ_l / γ_c)
end
lf_grid = 0.01:0.01:0.99
plot(lf_grid, l_f -> foc_bargaining(l_f, w_f=1.0), label=L"w_f = 1.0", xlabel=L"l_f")
plot!(lf_grid, l_f -> foc_bargaining(l_f, w_f=0.5), label=L"w_f = 0.5", linestyle=:dash)
plot!(lf_grid, l_f -> foc_bargaining(l_f, w_f=0.1), label=L"w_f = 0.1", linestyle=:dot)
hline!([0.0], lw=2, label=false, color=:black)
```
また, $l_f = 1$ における一階条件は以下のようになります.
$$
w_m(1-l_m) - 2\left(\frac{w_m}{\alpha_l}\right)^{\frac{1}{\gamma_c}} l_m^{\frac{\gamma_l}{\gamma_c}} = 0.
$$
この様な場合, Bracketing method の考えが有効です.
1. 端点 $l_m = 0, 1$ (実用上は $\epsilon, 1-\epsilon$) において, 一階条件の値を計算.
1. 一階条件の値が異符号であれば, 二分法を用いて解を求める.
1. 一階条件の値が同符号であれば, 端点解 (この場合, $l_f = 1$) 上の解を求める.
```{julia}
#| label: find-zero-bracketing-corner
#| output: false
function solve_bargaining(; w_m=1.0, w_f=0.5, γ_c=1.5, α_l=1.2, γ_l=1.5)
ϵ = 1e-9
left = w_m * (1 - (w_f / w_m)^(1 / γ_l) * ϵ) + w_f * (1 - ϵ) -
2 * (w_f / α_l)^(1 / γ_c) * ϵ^(γ_l / γ_c)
right = w_m * (1 - (w_f / w_m)^(1 / γ_l) * (1 - ϵ)) + w_f * ϵ -
2 * (w_f / α_l)^(1 / γ_c) * (1 - ϵ)^(γ_l / γ_c)
if left * right < 0
l_f = find_zero(
l_f -> w_m * (1 - (w_f / w_m)^(1 / γ_l) * l_f) + w_f * (1 - l_f) -
2 * (w_f / α_l)^(1 / γ_c) * l_f^(γ_l / γ_c),
(ϵ, 1 - ϵ))
l_m = (w_f / w_m)^(1 / γ_l) * l_f
else
l_f = 1.0
x = find_zero(
x -> w_m * (1 - 1 / (1 + exp(x))) -
2 * (w_m / α_l)^(1 / γ_c) * (1 / (1 + exp(x)))^(γ_l / γ_c),
0.0)
l_m = 1 / (1 + exp(x))
end
c = w_m * (1 - l_m) + w_f * (1 - l_f)
c_m = 0.5 * c
c_f = 0.5 * c
U = (c_m^(1 - γ_c) / (1 - γ_c)) + α_l * (l_m^(1 - γ_l) / (1 - γ_l)) +
(c_f^(1 - γ_c) / (1 - γ_c)) + α_l * (l_f^(1 - γ_l) / (1 - γ_l))
return (c_m=c_m, c_f=c_f, l_m=l_m, l_f=l_f, U=U)
end
```
```{julia}
#| code-fold: true
#| label: fig-sol-bargaining
#| fig-cap: "**Solution of bargaining problem.** $w_m = 1.0, \\alpha_l = 1.2, \\gamma_c = \\gamma_l = 1.5$."
plot(0.1:0.01:0.9, w_f -> solve_bargaining(w_f=w_f).l_m, label=L"l_m", xlabel=L"w_f")
plot!(0.1:0.01:0.9, w_f -> solve_bargaining(w_f=w_f).l_f, label=L"l_f", linestyle=:dash)
```