W06 실습 · 음함수 정리와 비교정학 - Understanding World on Equation

노트북은 예측 → 계산 → 대조 세 부분으로 고정한다. 새 개념은 도입하지 않는다. 본문에서 이미 유도한 것 — (eq-w06-7)·(eq-w06-11)·(eq-w06-15) — 을 수치로 확인할 뿐이다.

1. 예측 (코드를 쓰기 전에)¶

아래 세 항목의 답을 코드를 쓰기 전에 적는다. 계산하지 않는다.

1. $f(K)=K^{0.5}$ , $p=1$ , $r=0.5$ 에서 (eq-w06-7)은 $dK^{*}/dr=-4$ 다. 전진차분 $(K^{*}(r+h)-K^{*}(r))/h$ 의 오차는 $h$ 에 비례하고(1차), 중심차분의 오차는 $h^{2}$ 에 비례한다(2차). 구체 예측: $h=0.1$ 에서 전진 -3.06, 중심 -4.34; $h=0.01$ 에서 전진 -3.88, 중심 -4.003; $h=0.001$ 에서 전진 -3.99, 중심 -4.000. 계산은 $F(K,r)=0.5K^{-0.5}-r$ 의 근을 이분법으로 직접 찾는 함수로 한다 — 닫힌 해를 쓰지 않는다. 오차를 log-log로 그려 기울기 1과 2가 나오는지 본다(W03의 근사 차수).

예측: ____

2. $u=\ln x_1+\ln x_2$ , $p=(1,1)$ , $w=2$ 에서 $\partial x_1/\partial p_1=-1$ 이고 대체효과 -0.5, 소득효과 -0.5다. $\det\bar H=2$ . 계산은 (eq-w06-9)의 세 방정식을 뉴턴법으로 푼다 — 야코비안은 (eq-w06-8)의 $\bar H$ 를 직접 코드로 쌓는다. $p_1=1\pm0.01$ 에서 $x_1$ 을 구해 중심차분으로 $\partial x_1/\partial p_1$ 을 얻는다. 보상: $p_1=1.01$ 에서 $u(x(p',w'))=u_0$ 가 되는 $w'$ 를 이분법으로 찾아 $h_1$ 을 얻고 대체효과를, 나머지를 소득효과로 분해한다.

예측: ____

3. $\alpha\in\{0.5,0.8,0.95\}$ 에서 탄력성 -2, -5, -20. $\Delta r/r=+10\%$ 의 정확한 $\Delta K^{*}/K^{*}$ 는 $-17.4\%$ , $-37.9\%$ , $-85.1\%$ 이고 선형 예측은 $-20\%$ , $-50\%$ , $-200\%$ 다 — $F_x\to 0$ 일수록 도함수는 정확하되 유효한 근방이 사라진다. 계산은 항목 1의 근 찾기 함수를 재사용해 $\alpha$ 별로 $r=0.5$ 와 0.55의 근을 구하고 (eq-w06-7)의 도함수에 $\Delta r$ 을 곱한 것과 비교한다.

예측: ____

2. 계산¶

패키지 호출로 답을 내지 않는다. 1계 조건과 유계 헤시안을 원시 연산으로 직접 쌓고, 근은 이분법·뉴턴법으로 찾는다. 닫힌 해는 쓰지 않는다.

# 도구를 올린다 — 수치는 전부 아래에서 직접 쌓는다.
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 _c in _cands:
    if _c in _have:
        plt.rcParams['font.family'] = _c
        break
plt.rcParams['axes.unicode_minus'] = False
print('고른 글꼴:', plt.rcParams['font.family'][0])

고른 글꼴: Apple SD Gothic Neo

# (eq-w06-6)의 1계 조건 F(K,r)=p f'(K)-r 와 그 근을 이분법으로 쌓는다.
def F(K, r, alpha=0.5, p=1.0):
    return p * alpha * K ** (alpha - 1.0) - r


def root_K(r, alpha=0.5, p=1.0, lo=1e-12, hi=1e12):
    """F(K,r)=0 의 근. 닫힌 해를 쓰지 않는다 — 부호가 갈리는 구간을 반씩 줄인다."""
    a, b = lo, hi
    for _ in range(300):
        m = 0.5 * (a + b)
        if m <= a or m >= b:          # 인접한 부동소수에 닿으면 멈춘다
            break
        if F(m, r, alpha, p) > 0.0:   # F 는 K 에 대해 감소
            a = m
        else:
            b = m
    return 0.5 * (a + b)


K_star = root_K(0.5)
print(f'K*(r=0.5) = {K_star:.12f},  1계 조건 잔차 F = {F(K_star, 0.5):.3e}')

K*(r=0.5) = 1.000000000000,  1계 조건 잔차 F = 0.000e+00

# 예측 1 — 격자 h 마다 전진차분과 중심차분을 직접 쌓는다.
r0, exact = 0.5, -4.0
rows = []
for h in [0.1, 0.01, 0.001]:
    fwd = (root_K(r0 + h) - root_K(r0)) / h
    cen = (root_K(r0 + h) - root_K(r0 - h)) / (2.0 * h)
    rows.append((h, fwd, cen, abs(fwd - exact), abs(cen - exact)))

fpp = 0.5 * (0.5 - 1.0) * K_star ** (0.5 - 2.0)
print(f"(eq-w06-7) dK*/dr = 1/(p f''(K*)) = 1/({fpp:.4f}) = {1.0 / fpp:.4f}")
print(f"{'h':>7} | {'전진차분':>6} | {'중심차분':>6} | {'전진오차':>6} | {'중심오차':>6}")
for h, fwd, cen, ef, ec in rows:
    print(f'{h:>7.3f} | {fwd:>10.4f} | {cen:>10.4f} | {ef:>10.3e} | {ec:>10.3e}')

(eq-w06-7) dK*/dr = 1/(p f''(K*)) = 1/(-0.2500) = -4.0000
      h |   전진차분 |   중심차분 |   전진오차 |   중심오차
  0.100 |    -3.0556 |    -4.3403 |  9.444e-01 |  3.403e-01
  0.010 |    -3.8831 |    -4.0032 |  1.169e-01 |  3.202e-03
  0.001 |    -3.9880 |    -4.0000 |  1.197e-02 |  3.200e-05

# 오차의 차수를 log-log 기울기로 읽는다 — 1차와 2차.
h_arr = np.array([q[0] for q in rows])
err_f = np.array([q[3] for q in rows])
err_c = np.array([q[4] for q in rows])
slope_f = np.polyfit(np.log(h_arr), np.log(err_f), 1)[0]
slope_c = np.polyfit(np.log(h_arr), np.log(err_c), 1)[0]
print(f'log-log 기울기 — 전진 {slope_f:.3f}, 중심 {slope_c:.3f}')

fig, ax = plt.subplots(figsize=(5.4, 3.8))
ax.loglog(h_arr, err_f, 'o-', label=f'전진차분 (기울기 {slope_f:.2f})')
ax.loglog(h_arr, err_c, 's-', label=f'중심차분 (기울기 {slope_c:.2f})')
ax.set_xticks(h_arr, labels=['0.1', '0.01', '0.001'])   # 눈금을 직접 단다
ax.set_xticks([], minor=True)
_yt = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0]
ax.set_yticks(_yt, labels=['1e-5', '1e-4', '1e-3', '1e-2', '1e-1', '1'])
ax.set_yticks([], minor=True)
ax.set_xlabel('h')
ax.set_ylabel('|차분 - (-4)|')
ax.set_title('근사 차수 — 기울기 1과 2')
ax.grid(True, which='both', alpha=0.3)
ax.legend()
plt.show()

log-log 기울기 — 전진 0.949, 중심 2.013

# (eq-w06-8) 유계 헤시안과 (eq-w06-9) 1계 조건계를 직접 쌓는다. u = ln x1 + ln x2.
def foc(z, p1, p2, w):
    x1, x2, lam = z
    return np.array([1.0 / x1 - lam * p1,
                     1.0 / x2 - lam * p2,
                     w - p1 * x1 - p2 * x2])


def bordered(z, p1, p2):
    x1, x2, _ = z
    return np.array([[-1.0 / x1 ** 2, 0.0, -p1],
                     [0.0, -1.0 / x2 ** 2, -p2],
                     [-p1, -p2, 0.0]])


def solve_consumer(p1, p2, w, z0=(1.0, 1.0, 1.0)):
    """뉴턴법 — 야코비안은 위 유계 헤시안 그대로다."""
    z = np.array(z0, dtype=float)
    for _ in range(200):
        dz = np.linalg.solve(bordered(z, p1, p2), -foc(z, p1, p2, w))
        z = z + dz
        if np.max(np.abs(dz)) < 1e-14:
            break
    return z


print('뉴턴 해 (x1, x2, λ) =', np.round(solve_consumer(1.0, 1.0, 2.0), 12))

뉴턴 해 (x1, x2, λ) = [1. 1. 1.]

# 예측 2 — 해와 det H̄, 그리고 중심차분으로 얻은 마셜 도함수.
p1, p2, w = 1.0, 1.0, 2.0
z0 = solve_consumer(p1, p2, w)
x1_0, x2_0, lam0 = z0
Hbar = bordered(z0, p1, p2)
detH = float(np.linalg.det(Hbar))
print(f'해  x = ({x1_0:.10f}, {x2_0:.10f}),  λ = {lam0:.10f}')
print('유계 헤시안 H̄ =\n', Hbar)
print(f'(eq-w06-11) det H̄ = {detH:.10f}   [= 2p1p2u12 − p1²u22 − p2²u11]')

d = 0.01
dx1_dp1 = (solve_consumer(p1 + d, p2, w)[0] - solve_consumer(p1 - d, p2, w)[0]) / (2 * d)
dx1_dw = (solve_consumer(p1, p2, w + d)[0] - solve_consumer(p1, p2, w - d)[0]) / (2 * d)
print(f'∂x1/∂p1 (중심차분) = {dx1_dp1:.6f}')
print(f'∂x1/∂w  (중심차분) = {dx1_dw:.6f}')

해  x = (1.0000000000, 1.0000000000),  λ = 1.0000000000
유계 헤시안 H̄ =
 [[-1.  0. -1.]
 [ 0. -1. -1.]
 [-1. -1.  0.]]
(eq-w06-11) det H̄ = 2.0000000000   [= 2p1p2u12 − p1²u22 − p2²u11]
∂x1/∂p1 (중심차분) = -1.000100
∂x1/∂w  (중심차분) = 0.500000

# 보상소득을 이분법으로 찾아 (eq-w06-15)의 두 항으로 가른다.
u0 = np.log(x1_0) + np.log(x2_0)


def comp_income(p1n, p2n=1.0, lo=0.1, hi=10.0):
    """u(x(p',w')) = u0 가 되는 w'. 효용은 w' 에 증가하므로 이분법이 닿는다."""
    for _ in range(300):
        m = 0.5 * (lo + hi)
        if m <= lo or m >= hi:
            break
        zz = solve_consumer(p1n, p2n, m)
        if np.log(zz[0]) + np.log(zz[1]) < u0:
            lo = m
        else:
            hi = m
    return 0.5 * (lo + hi)


h_up = solve_consumer(p1 + d, p2, comp_income(p1 + d))[0]
h_dn = solve_consumer(p1 - d, p2, comp_income(p1 - d))[0]
subst = (h_up - h_dn) / (2 * d)
income = dx1_dp1 - subst
print(f'대체효과 ∂h1/∂p1   = {subst:>9.6f}   [−λp2²/det H̄ = {-lam0 * p2 ** 2 / detH:.6f}]')
print(f'소득효과 −x1·∂x1/∂w = {income:>9.6f}   [−x1·∂x1/∂w  = {-x1_0 * dx1_dw:.6f}]')
print(f'두 항의 합 = {subst + income:.6f}   vs   ∂x1/∂p1 = {dx1_dp1:.6f}')

대체효과 ∂h1/∂p1   = -0.500031   [−λp2²/det H̄ = -0.500000]
소득효과 −x1·∂x1/∂w = -0.500069   [−x1·∂x1/∂w  = -0.500000]
두 항의 합 = -1.000100   vs   ∂x1/∂p1 = -1.000100

# 마셜 수요와 보상 수요를 한 판에 올린다 — 기울기의 차이가 소득효과다.
grid = np.linspace(0.8, 1.3, 26)
marsh = np.array([solve_consumer(g, p2, w)[0] for g in grid])
hicks = np.array([solve_consumer(g, p2, comp_income(g))[0] for g in grid])

fig, ax = plt.subplots(figsize=(5.4, 3.8))
ax.plot(grid, marsh, label='마셜 수요 $x_1(p_1, w=2)$')
ax.plot(grid, hicks, label='보상 수요 $h_1(p_1, u_0)$')
ax.plot([p1], [x1_0], 'ko', ms=4)
ax.set_xlabel('$p_1$')
ax.set_ylabel('$x_1$')
ax.set_title(f'(eq-w06-15)  전체 {dx1_dp1:.3f} = 대체 {subst:.3f} + 소득 {income:.3f}')
ax.grid(alpha=0.3)
ax.legend()
plt.show()

# 예측 3 — α별로 근을 다시 찾아 정확한 변화와 선형 예측을 견준다.
alphas = [0.5, 0.8, 0.95]
exact_pct, lin_pct = [], []
print(f"{'α':>5} | {'탄력성':>4} | {'K*(0.5)':>12} | {'정확 ΔK*/K*':>10} | "
      f"{'선형 예측':>6} | {'격차(%p)':>7}")
for a in alphas:
    K0 = root_K(0.5, alpha=a)
    K1 = root_K(0.55, alpha=a)
    slope = 1.0 / (a * (a - 1.0) * K0 ** (a - 2.0))   # (eq-w06-7) 의 1/(p f''(K*))
    ex = 100.0 * (K1 - K0) / K0
    li = 100.0 * slope * 0.05 / K0
    exact_pct.append(ex)
    lin_pct.append(li)
    print(f'{a:>5.2f} | {-1.0 / (1.0 - a):>7.1f} | {K0:>12.4f} | {ex:>11.1f}% | '
          f'{li:>9.1f}% | {ex - li:>9.1f}')

    α |  탄력성 |      K*(0.5) |  정확 ΔK*/K* |  선형 예측 |  격차(%p)
 0.50 |    -2.0 |       1.0000 |       -17.4% |     -20.0% |       2.6
 0.80 |    -5.0 |      10.4858 |       -37.9% |     -50.0% |      12.1
 0.95 |   -20.0 |  375899.7346 |       -85.1% |    -200.0% |     114.9

# α가 1에 다가갈수록 선형 예측이 어디서 무너지는지 본다.
fig, ax = plt.subplots(figsize=(5.8, 3.9))
r_grid = np.linspace(0.5, 0.75, 26)
for a, c in zip(alphas, ['C0', 'C1', 'C2']):
    K0 = root_K(0.5, alpha=a)
    ex = np.array([100.0 * (root_K(x, alpha=a) - K0) / K0 for x in r_grid])
    slope = 1.0 / (a * (a - 1.0) * K0 ** (a - 2.0))
    ax.plot(r_grid, ex, color=c, label=f'α={a} 정확')
    ax.plot(r_grid, 100.0 * slope * (r_grid - 0.5) / K0, color=c, ls='--',
            label=f'α={a} 선형')
ax.axvline(0.55, color='k', lw=0.8)
ax.set_ylim(-250, 15)
ax.set_xlabel('$r$')
ax.set_ylabel('$\\Delta K^{*}/K^{*}$ (%)')
ax.set_title('도함수는 정확하되 유효한 근방이 좁아진다')
ax.grid(alpha=0.3)
ax.legend(fontsize=8, ncol=3)
plt.show()

3. 대조¶

예측과 계산이 어긋난 지점을 적는다. 어느 쪽이 틀렸는지 판정한다.

예측	결과	어긋남	원인

본문 확인: 어긋났으면 (eq-w06-7)·(eq-w06-11)·(eq-w06-15)로 돌아간다.