W02 실습 · 차원과 단위

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 F(K, L, A=1.0, a=0.3, b=0.7):
    return A * K ** a * L ** b

A0, ALPHA = 1.0, 0.3
print('글꼴:', plt.rcParams['font.family'][0])
print('F(100, 100) =', round(F(100.0, 100.0), 4))

# 예측 1 — L을 고정한 채 K를 2배 했을 때의 산출 배율을 쌓는다
L0 = 100.0
K1, K2 = 100.0, 200.0
Y1, Y2 = F(K1, L0), F(K2, L0)
ratio = Y2 / Y1
print(f'Y(K=100) = {Y1:.4f}     Y(K=200) = {Y2:.4f}')
print(f'배율  Y2/Y1    = {ratio:.4f}      2^0.3 = {2.0 ** ALPHA:.4f}')
print(f'증가율         = {100 * (ratio - 1):.2f}%     (30%가 아니다)')
print(f'Δln Y / Δln K  = {(np.log(Y2) - np.log(Y1)) / (np.log(K2) - np.log(K1)):.6f}')

# 같은 배율이 K의 수준을 모르는지 여러 자리에서 쌓아 본다
print(f"{'K':>9}{'Y(K)':>13}{'Y(2K)':>13}{'Y(2K)/Y(K)':>14}")
for K in [10.0, 100.0, 1000.0, 10000.0]:
    print(f'{K:9.0f}{F(K, L0):13.4f}{F(2 * K, L0):13.4f}{F(2 * K, L0) / F(K, L0):14.6f}')
print()
print('덧셈 충격 K→K+100의 배율 — 수준에 의존한다')
print(f'  K=100: {F(200.0, L0) / F(100.0, L0):.6f}     K=1000: {F(1100.0, L0) / F(1000.0, L0):.6f}')

# 전방 유한차분으로 한계생산을 쌓는다 — 미분 패키지를 쓰지 않는다
def FK_fd(K, L, h, A=1.0, a=0.3, b=0.7):
    return (F(K + h, L, A, a, b) - F(K, L, A, a, b)) / h


def FL_fd(K, L, h, A=1.0, a=0.3, b=0.7):
    return (F(K, L + h, A, a, b) - F(K, L, A, a, b)) / h


K0, Lb = 300.0, 100.0
print(f'F_K 전방차분(h=1e-3) = {FK_fd(K0, Lb, 1e-3):.10f}')
print(f'F_K 정확값 αY/K      = {ALPHA * F(K0, Lb) / K0:.10f}')
print(f'F_L 전방차분(h=1e-3) = {FL_fd(K0, Lb, 1e-3):.10f}')
print(f'F_L 정확값 βY/L      = {0.7 * F(K0, Lb) / Lb:.10f}')

# 예측 2 — ν별로 KF_K, LF_L, 합, νY, 잔여를 쌓아 3절 표와 대조한다
h = 1e-3
rows = []
print(f"{'ν':>5}{'β':>6}{'Y':>11}{'K·F_K':>11}{'L·F_L':>11}{'합':>11}{'νY':>11}{'잔여':>9}{'상대어긋남':>14}")
for b in [0.6, 0.7, 0.8]:
    nu = ALPHA + b
    Y = F(K0, Lb, A0, ALPHA, b)
    kfk = K0 * FK_fd(K0, Lb, h, A0, ALPHA, b)
    lfl = Lb * FL_fd(K0, Lb, h, A0, ALPHA, b)
    tot = kfk + lfl
    rel = (tot - nu * Y) / (nu * Y)
    rows.append((nu, Y, kfk, lfl, tot, rel))
    print(f'{nu:5.1f}{b:6.1f}{Y:11.4f}{kfk:11.4f}{lfl:11.4f}{tot:11.4f}'
          f'{nu * Y:11.4f}{(1 - nu) * Y:9.2f}{rel:14.3e}')
print()
print('자본의 몫 K·F_K/Y :', [round(r[2] / r[1], 6) for r in rows])

# h를 10배씩 줄이며 어긋남이 어떻게 주는지 쌓는다 (ν=1)
hs, rels = [1e-1, 1e-2, 1e-3, 1e-4], []
print(f"{'h':>9}{'합 - νY':>15}{'상대어긋남':>17}{'직전 대비':>13}")
for h_ in hs:
    Y = F(K0, Lb)
    tot = K0 * FK_fd(K0, Lb, h_) + Lb * FL_fd(K0, Lb, h_)
    rel = (tot - Y) / Y
    tail = '' if not rels else f'{rel / rels[-1]:13.3f}'
    rels.append(rel)
    print(f'{h_:9.0e}{tot - Y:15.3e}{rel:17.3e}{tail}')

# 생산물 대비 요소지급을 쌓아 그린다
sk = [r[2] / r[1] for r in rows]
sl = [r[3] / r[1] for r in rows]
x = np.arange(3)
fig, ax = plt.subplots(figsize=(5.4, 3.4))
ax.bar(x, sk, 0.5, color='#3b6ea5', label='자본 몫 $KF_K/Y$')
ax.bar(x, sl, 0.5, bottom=sk, color='#d98c3f', label='노동 몫 $LF_L/Y$')
ax.axhline(1.0, ls='--', lw=1.0, color='k')
for i, (u, v) in enumerate(zip(sk, sl)):
    ax.text(i, u + v + 0.03, f'{u + v:.2f}', ha='center', fontsize=9)
ax.set_xticks(x)
ax.set_xticklabels([f'ν={r[0]:.1f}' for r in rows])
ax.set_ylabel('생산물 대비 요소지급')
ax.set_ylim(0, 1.35)
ax.legend(loc='lower right', fontsize=8)
ax.set_title('합이 생산물과 만나는 것은 ν=1뿐이다')
plt.tight_layout()
plt.show()

# 예측 3 — 아홉 점 격자 자료와 정규방정식을 직접 세운다
Ks, Ls = np.array([100.0, 300.0, 1000.0]), np.array([50.0, 100.0, 200.0])
KK, LL = np.meshgrid(Ks, Ls, indexing='ij')
Kg, Lg = KK.ravel(), LL.ravel()
Yg = F(Kg, Lg)


def ols(K_col, L_col, y):
    X = np.column_stack([np.ones_like(K_col), np.log(K_col), np.log(L_col)])
    return np.linalg.solve(X.T @ X, X.T @ np.log(y))


c1, a1, b1 = ols(Kg, Lg, Yg)
print('격자 아홉 점 (K, L, Y)')
for K, L, Y in zip(Kg, Lg, Yg):
    print(f'  {K:7.0f}{L:7.0f}{Y:11.4f}')
print(f'원 단위: ĉ = {c1:+.6f}  α̂ = {a1:.12f}  β̂ = {b1:.12f}  Â = {np.exp(c1):.6f}')

# 같은 자료를 K/1000으로 다시 재서 추정한다
c2, a2, b2 = ols(Kg / 1000.0, Lg, Yg)
print(f'천 단위: ĉ = {c2:+.6f}  α̂ = {a2:.12f}  β̂ = {b2:.12f}  Â = {np.exp(c2):.6f}')
print(f'ĉ 이동량 = {c2 - c1:.6f}      α·ln1000 = {ALPHA * np.log(1000.0):.6f}')
print(f'Â 배율   = {np.exp(c2) / np.exp(c1):.4f}         1000^0.3 = {1000.0 ** ALPHA:.4f}')
print(f'α̂ 변화   = {a2 - a1:+.3e}      β̂ 변화 = {b2 - b1:+.3e}')

# 분배율은 단위에 비례하는 증분으로 쌓는다
sK = {}
for nm, K_, A_ in [('원 단위', 300.0, 1.0), ('천 단위', 0.3, 1000.0 ** ALPHA)]:
    sK[nm] = K_ * FK_fd(K_, Lb, 1e-6 * K_, A_) / F(K_, Lb, A_)
    print(f'{nm}: K = {K_:8.4f}  Y = {F(K_, Lb, A_):.4f}  s_K = {sK[nm]:.6f}')

# 같은 자료를 semi-log와 log-log 두 축에 찍는다
Kc = np.logspace(1, 3, 200)
m = Lg == 100.0
fig, axes = plt.subplots(1, 2, figsize=(9.6, 3.6))
axes[0].plot(Kc, np.log(F(Kc, 100.0)), color='k')
axes[0].scatter(Kg[m], np.log(Yg[m]), s=20, color='k', zorder=3)
axes[0].set_xlabel('자본 $K$')
axes[0].set_ylabel('$\\ln Y$')
axes[0].set_title('(a) semi-log — 오목 곡선')
axes[1].plot(np.log(Kc), np.log(F(Kc, 100.0)), color='#3b6ea5', label='원 단위')
axes[1].plot(np.log(Kc / 1000.0), np.log(F(Kc, 100.0)), ls='--', color='#c0392b', label='천 단위')
axes[1].scatter(np.log(Kg[m]), np.log(Yg[m]), s=20, color='k', zorder=3)
axes[1].annotate('', xy=(np.log(100.0), np.log(F(100.0, 100.0))),
                 xytext=(np.log(0.1), np.log(F(100.0, 100.0))),
                 arrowprops=dict(arrowstyle='<->', lw=0.9))
axes[1].text(np.log(3.0), np.log(F(100.0, 100.0)) + 0.08, 'ln1000 = 6.91', ha='center', fontsize=8)
axes[1].set_xlabel('$\\ln K$')
axes[1].set_ylabel('$\\ln Y$')
axes[1].set_title('(b) log-log — 기울기 0.3, 단위 변경은 평행이동')
axes[1].legend(loc='lower right', fontsize=8)
plt.tight_layout()
plt.show()

# 예측 항목 세 개에 대응하는 수치를 한자리에 모은다
d_ln = (np.log(Y2) - np.log(Y1)) / (np.log(K2) - np.log(K1))
print(f'[예측 1] Y 배율 {ratio:.4f} (= 2^0.3), 증가율 {100 * (ratio - 1):.2f}%, Δln Y/Δln K = {d_ln:.6f}')
print('[예측 2] h=1e-3 상대 어긋남  ' + ',  '.join(f'ν={r[0]:.1f}: {r[5]:+.3e}' for r in rows))
print(f'         h 10배 축소 시 오차 비 {rels[1] / rels[0]:.3f}, {rels[2] / rels[1]:.3f}, '
      f'{rels[3] / rels[2]:.3f}')
print(f'[예측 3] Â {np.exp(c1):.4f} → {np.exp(c2):.4f} ({np.exp(c2) / np.exp(c1):.4f}배), '
      f'α̂ 변화 {a2 - a1:+.2e}, s_K {sK["원 단위"]:.6f} → {sK["천 단위"]:.6f}')