W04 실습 · 라그랑주 승수는 가격이다

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager

# 한국어 글꼴을 후보 목록에서 있는 것으로 고른다
_cands = ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic', 'Noto Sans CJK KR', 'Malgun Gothic']
_have = {f.name for f in font_manager.fontManager.ttflist}
for _f in _cands:
    if _f in _have:
        plt.rcParams['font.family'] = _f
        break
plt.rcParams['axes.unicode_minus'] = False

# 세 효용함수를 손으로 쌓는다 — 로그의 경계는 아주 작은 수로 밀어 둔다
def ln(x):
    return np.log(np.maximum(np.asarray(x, dtype=float), 1e-300))


def u_log(x1, x2, a=0.5):
    return a * ln(x1) + (1.0 - a) * ln(x2)


def u_qlin(x1, x2):
    return ln(x1) + np.asarray(x2, dtype=float)


def u_pow(x1, x2, a=0.5):
    return np.asarray(x1, dtype=float) ** a * np.asarray(x2, dtype=float) ** (1.0 - a)


A, P1, P2, W0 = 0.5, 1.0, 2.0, 100.0
print('글꼴:', plt.rcParams['font.family'][0])
print(f'u_log(50, 25) = {u_log(50.0, 25.0):.6f}     u_pow(50, 25) = {u_pow(50.0, 25.0):.6f}')

# 예산선 위로 1변수화해 격자를 좁혀 가며 최대점을 찾는다 (양 끝 모서리 포함)
def solve(u, p, w, q=None, n=801, rounds=6):
    hi0 = w / p[0] if q is None else min(q, w / p[0])
    lo, hi = 0.0, hi0
    for _ in range(rounds + 1):
        xs = np.linspace(lo, hi, n)
        vals = u(xs, np.maximum((w - p[0] * xs) / p[1], 0.0))
        k = int(np.argmax(vals))
        h = (hi - lo) / (n - 1)
        lo, hi = max(0.0, xs[k] - h), min(hi0, xs[k] + h)
    return np.array([xs[k], max((w - p[0] * xs[k]) / p[1], 0.0)]), float(vals[k])


_x, _v = solve(u_log, (P1, P2), W0)
print(f'검산 solve : x*(w=100) = ({_x[0]:.6f}, {_x[1]:.6f}), v = {_v:.6f}   닫힌형 (50, 25), 3.565449')

# 한계효용은 차분으로 쌓고, 승수는 1계 조건에서 읽는다
def grad(u, x, h=1e-6):
    g = np.zeros(2)
    for i in (0, 1):
        s = h * max(1.0, abs(x[i]))
        a, b = np.array(x, dtype=float), np.array(x, dtype=float)
        b[i] += s
        a[i] = max(0.0, a[i] - s)
        g[i] = (float(u(*b)) - float(u(*a))) / (b[i] - a[i])
    return g


def mults(u, p, x):
    g = grad(u, x)
    lam = float(np.max(g / np.asarray(p, dtype=float)))
    return lam, lam * np.asarray(p, dtype=float) - g


def dvdw(u, p, w, h=1e-3):
    return (solve(u, p, w + h)[1] - solve(u, p, w - h)[1]) / (2.0 * h)


def pad(s, n):
    return ' ' * max(0, n - sum(2 if '\uac00' <= c <= '\ud7a3' else 1 for c in s)) + s


_g = grad(u_log, np.array([50.0, 25.0]))
print(f'검산 grad  : ∇u_log(50, 25) = ({_g[0]:.8f}, {_g[1]:.8f})   닫힌형 (0.01, 0.02)')
print(f'검산 mults : λ* = {mults(u_log, (P1, P2), np.array([50.0, 25.0]))[0]:.8f}   닫힌형 0.01')

# 기준점 — 격자 해와 닫힌형 (eq-w04-13)을 맞춘다
x_log, v0 = solve(u_log, (P1, P2), W0)
lam0, mu0 = mults(u_log, (P1, P2), x_log)
g0 = grad(u_log, x_log)
v_cf = np.log(W0) + A * np.log(A / P1) + (1 - A) * np.log((1 - A) / P2)
print(f'x* (격자)      = ({x_log[0]:.8f}, {x_log[1]:.8f})      닫힌형 (50, 25)')
print(f'v(p, 100)      = {v0:.8f}                     닫힌형 {v_cf:.8f}')
print(f'λ* = u_1/p_1   = {lam0:.8f}                     1/w = {1.0 / W0:.8f}')
print(f'∂v/∂w (중심차분) = {dvdw(u_log, (P1, P2), W0):.8f}')
print(f'μ* = λ*p − ∇u  = ({mu0[0]:+.2e}, {mu0[1]:+.2e})        내부해이므로 0')
print(f'1원당 한계효용  u_1/p_1 = {g0[0] / P1:.8f},  u_2/p_2 = {g0[1] / P2:.8f}   (eq-w04-17)')

# 같은 것을 소득 격자 위에서 세 길로 겨룬다 — 1계 조건, 중심차분, 닫힌형 1/w
print()
print(pad('w', 7) + pad('x1*', 10) + pad('x2*', 10) + pad('u_1/p_1', 12)
      + pad('u_2/p_2', 12) + pad('∂v/∂w', 12) + pad('1/w', 12) + pad('최대 어긋남', 14))
for w in [25.0, 50.0, 100.0, 200.0, 400.0]:
    xw, _ = solve(u_log, (P1, P2), w)
    gw = grad(u_log, xw)
    dw = dvdw(u_log, (P1, P2), w)
    gap = max(abs(gw[0] / P1 - 1.0 / w), abs(gw[1] / P2 - 1.0 / w), abs(dw - 1.0 / w))
    print(f'{w:7.0f}{xw[0]:10.4f}{xw[1]:10.4f}{gw[0] / P1:12.6f}{gw[1] / P2:12.6f}'
          f'{dw:12.6f}{1.0 / w:12.6f}{gap:14.2e}')

# 예측 1 — 유한차분 몫이 λ*의 어느 쪽에 놓이는가, 오차는 Δ에 어떻게 비례하는가
print(pad('Δ', 7) + pad('차분 몫', 15) + pad('오차', 14) + pad('2차항 −Δ/(2w²)', 18))
rows1 = []
for d in [0.1, 1.0, 10.0, 50.0, -50.0]:
    qd = (solve(u_log, (P1, P2), W0 + d)[1] - v0) / d
    rows1.append((d, qd, qd - lam0))
    print(f'{d:7.1f}{qd:15.6f}{qd - lam0:14.6f}{-d / (2.0 * W0 ** 2):18.6f}')

q1 = rows1[1][1]
r10 = [rows1[i + 1][2] / rows1[i][2] for i in (0, 1)]
print()
print(f'[v(101) − v(100)]/1 = {q1:.6f}  {"<" if q1 < lam0 else ">"}  λ* = {lam0:.6f}'
      f'   (할선은 접선 {"아래" if q1 < lam0 else "위"})')
print(f'Δ 10배 시 오차 배율 : 0.1→1 에서 {r10[0]:.3f}배,  1→10 에서 {r10[1]:.3f}배')
print(f'Δ = −50 의 몫 {rows1[4][1]:.6f} > λ* — 같은 크기의 소득 감소가 더 크게 움직인다')

# 가치함수의 접선·할선과 오차의 기울기를 그린다
wg = np.linspace(60.0, 160.0, 61)
vg = np.array([solve(u_log, (P1, P2), ww)[1] for ww in wg])
ds = np.logspace(-2.0, 1.7, 16)
errs = np.array([abs((solve(u_log, (P1, P2), W0 + d)[1] - v0) / d - lam0) for d in ds])
w1, v1 = 150.0, solve(u_log, (P1, P2), 150.0)[1]

fig, ax = plt.subplots(1, 2, figsize=(9.6, 3.6))
ax[0].plot(wg, vg, lw=2.0, color='#3b6ea5', label='$v(p,w)$')
ax[0].plot(wg, v0 + lam0 * (wg - W0), '--', lw=1.2, color='#c0392b',
           label='접선 기울기 $\\lambda^*=%.4f$' % lam0)
ax[0].plot([W0, w1], [v0, v1], lw=1.4, color='#2e8b57',
           label='할선 기울기 %.4f' % ((v1 - v0) / (w1 - W0)))
ax[0].plot([W0], [v0], 'o', color='k', ms=4)
ax[0].set_xlabel('소득 $w$')
ax[0].set_ylabel('가치함수 $v$')
ax[0].legend(fontsize=8, loc='lower right')
ax[0].set_title('(a) 할선은 접선 아래에 놓인다')
ax[1].loglog(ds, errs, 'o-', lw=1.6, ms=3.5, color='#3b6ea5', label='|차분 몫 $-\\lambda^*$|')
ax[1].loglog(ds, ds / (2.0 * W0 ** 2), ':', lw=1.4, color='k',
             label='2차항 $\\Delta/(2w^2)$ — 기울기 1')
ax[1].set_xlabel('소득 변화 $\\Delta$')
ax[1].set_ylabel('오차')
ax[1].tick_params(labelfontfamily='DejaVu Sans')   # 로그 눈금의 - 는 라틴 글꼴로 찍는다
ax[1].legend(fontsize=8, loc='upper left')
ax[1].set_title('(b) 오차는 $\\Delta$에 비례한다')
plt.tight_layout()
plt.show()

# 예측 2 — 준선형에서 두 승수를 소득별로 쌓는다
print(pad('w', 7) + pad('x1*', 9) + pad('x2*', 9) + pad('λ*', 11) + pad('∂v/∂w', 11)
      + pad('μ2*', 10) + pad('μ2*/λ*', 11) + pad('v', 11))
for w in [0.5, 1.0, 1.5, 2.0, 4.0, 8.0]:
    xw, vw = solve(u_qlin, (P1, P2), w)
    lam, mu = mults(u_qlin, (P1, P2), xw)
    print(f'{w:7.2f}{xw[0]:9.4f}{xw[1]:9.4f}{lam:11.6f}{dvdw(u_qlin, (P1, P2), w):11.6f}'
          f'{mu[1]:10.6f}{mu[1] / lam:11.4f}{vw:11.6f}')

# 촘촘한 소득 격자에서 μ2*가 꺼지는 문턱을 찾는다
wq = np.linspace(0.3, 8.0, 155)
lam_q = np.zeros(wq.size)
mu2_q = np.zeros(wq.size)
for i, w in enumerate(wq):
    lam_q[i], mu = mults(u_qlin, (P1, P2), solve(u_qlin, (P1, P2), w)[0])
    mu2_q[i] = mu[1]
on = mu2_q > 1e-6          # 격자·차분의 바닥 잡음(~1e-8)보다 큰 것만 센다
lam_hi = float(lam_q[wq >= 2.0].mean())
print()
print(f'μ2* > 0 인 구간 : w ≤ {wq[on].max():.2f},  μ2* = 0 이 시작되는 격자점 w = {wq[~on].min():.2f}'
      f'   (격자 간격 {wq[1] - wq[0]:.2f}, 본문 문턱 w̄ = 2)')
print(f'λ* : w < 2 에서 1/w (최대 어긋남 {np.max(np.abs(lam_q[wq < 2.0] - 1.0 / wq[wq < 2.0])):.2e}), '
      f'w ≥ 2 에서 {lam_hi:.6f} 고정')

# 예측 3 — 같은 선호의 지수형 표현을 소득별로 겨룬다
print(pad('w', 7) + pad('x1*(로그)', 13) + pad('x1*(지수)', 13) + pad('λ*(로그)', 13)
      + pad('λ*(지수)', 13) + pad('φ\'·λ*(로그)', 15) + pad('∂v/∂w(지수)', 15))
for w in [25.0, 50.0, 100.0, 200.0]:
    xl, vl = solve(u_log, (P1, P2), w)
    xp, vp = solve(u_pow, (P1, P2), w)
    ll = mults(u_log, (P1, P2), xl)[0]
    lp = mults(u_pow, (P1, P2), xp)[0]
    print(f'{w:7.0f}{xl[0]:13.4f}{xp[0]:13.4f}{ll:13.6f}{lp:13.6f}'
          f'{np.exp(vl) * ll:15.6f}{dvdw(u_pow, (P1, P2), w):15.6f}')

x_pow, v_pow = solve(u_pow, (P1, P2), W0)
lam_pow = mults(u_pow, (P1, P2), x_pow)[0]
print()
print(f'w = 100 : x*(지수) = ({x_pow[0]:.4f}, {x_pow[1]:.4f}), '
      f'x*(로그) = ({x_log[0]:.4f}, {x_log[1]:.4f}) — 같다')
print(f'          λ*(지수) = {lam_pow:.6f} = φ\'(v)·λ*(로그) = {np.exp(v0) * lam0:.6f}, '
      f'로그형의 {lam_pow / lam0:.2f}배')
print(f'          v(지수)/w = {v_pow / W0:.6f} — w에 무관한 상수다')

# 보론(4절 Q2) — 배급 상한 x1 ≤ q 의 승수 ν를 1계 조건과 차분 두 길로 쌓는다
print(pad('q', 6) + pad('x1*', 9) + pad('x2*', 9) + pad('λ*', 11) + pad('ν (1계)', 12)
      + pad('ν (차분)', 13) + pad('ν/λ*', 10) + pad('지불용의', 11) + pad('v', 11))
for q in [20.0, 40.0, 50.0, 60.0]:
    xq, vq = solve(u_log, (P1, P2), W0, q=q)
    gq = grad(u_log, xq)
    lam = gq[1] / P2                       # 배급이 걸리지 않은 재화 2의 1계 조건
    nu = max(gq[0] - lam * P1, 0.0)        # 재화 1의 부족분이 배급 제약의 가격
    h = 1e-3
    nud = (solve(u_log, (P1, P2), W0, q=q + h)[1]
           - solve(u_log, (P1, P2), W0, q=q - h)[1]) / (2.0 * h)
    print(f'{q:6.1f}{xq[0]:9.4f}{xq[1]:9.4f}{lam:11.6f}{nu:12.6f}{nud:13.6f}'
          f'{nu / lam:10.4f}{P1 + nu / lam:11.4f}{vq:11.6f}')
print()
print('q > 50 이면 배급이 묶이지 않고 ν = 0 이다 — 묶이지 않는 제약의 가격은 0 (eq-w04-3)')

# 두 승수의 소득 경로와 척도를 바꾼 승수를 그린다
wp = np.linspace(20.0, 200.0, 91)
lam_l = np.array([mults(u_log, (P1, P2), solve(u_log, (P1, P2), ww)[0])[0] for ww in wp])
lam_p = np.array([mults(u_pow, (P1, P2), solve(u_pow, (P1, P2), ww)[0])[0] for ww in wp])

fig, ax = plt.subplots(1, 2, figsize=(9.6, 3.6))
ax[0].plot(wq, lam_q, lw=2.0, color='#3b6ea5', label='$\\lambda^*(w)$ — 예산의 가격')
ax[0].plot(wq, mu2_q, lw=2.0, color='#c0392b', label='$\\mu_2^*(w)$ — 비음 제약의 가격')
ax[0].axvline(2.0, ls=':', lw=1.0, color='0.4')
ax[0].text(2.1, 2.6, '$\\bar w=2$', fontsize=9)
ax[0].set_xlabel('소득 $w$')
ax[0].set_ylabel('승수')
ax[0].legend(fontsize=8)
ax[0].set_title('(a) 준선형 — 문턱에서 $\\mu_2^*$가 꺼진다')
ax[1].semilogy(wp, lam_l, lw=2.0, color='#3b6ea5', label='로그형 $\\lambda^*=1/w$')
ax[1].semilogy(wp, lam_p, lw=2.0, color='#2e8b57',
               label='지수형 $\\lambda^*=%.4f$ (상수)' % lam_p.mean())
ax[1].set_xlabel('소득 $w$')
ax[1].set_ylabel('$\\lambda^*$')
ax[1].tick_params(labelfontfamily='DejaVu Sans')   # 로그 눈금의 - 는 라틴 글꼴로 찍는다
ax[1].legend(fontsize=8)
ax[1].set_title('(b) 같은 $x^*$, 다른 $\\lambda^*$ — 척도의 선택')
plt.tight_layout()
plt.show()

# 예측 항목 세 개에 대응하는 수치를 한자리에 모은다
print(f'[예측 1] [v(101)−v(100)]/1 = {q1:.6f} < λ* = {lam0:.6f}  (오차 {q1 - lam0:+.6f})')
print(f'         Δ 10배 → 오차 {r10[0]:.2f}배, {r10[1]:.2f}배 — 오차는 Δ에 비례한다')
print(f'[예측 2] λ* : w < 2 에서 1/w, w ≥ 2 에서 {lam_hi:.4f} 고정 '
      f'(w = 0.5, 1, 2, 8 에서 {", ".join(f"{mults(u_qlin, (P1, P2), solve(u_qlin, (P1, P2), w)[0])[0]:.3f}" for w in [0.5, 1.0, 2.0, 8.0])})')
print(f'         μ2* > 0 인 구간 : w ≤ {wq[on].max():.2f} (문턱 w̄ = 2), '
      f'w = 1 에서 μ2* = {mults(u_qlin, (P1, P2), solve(u_qlin, (P1, P2), 1.0)[0])[1][1]:.4f}')
print(f'[예측 3] x*(지수) = ({x_pow[0]:.4f}, {x_pow[1]:.4f}) = x*(로그) — 바뀌지 않는다')
print(f'         λ*(지수) = {lam_pow:.6f} 상수 (로그형 {lam0:.6f}의 {lam_pow / lam0:.2f}배 = φ\'(v))')