W10 실습 · 위상선과 수렴속도

# 배열·그림 도구와 한국어 글꼴을 얹는다
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm

CAND = ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic', 'Noto Sans CJK KR', 'Malgun Gothic']
AVAIL = {f.name for f in fm.fontManager.ttflist}
FONT = next((f for f in CAND if f in AVAIL), 'DejaVu Sans')
plt.rcParams['font.family'] = FONT
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['mathtext.fontset'] = 'cm'
print('글꼴:', FONT)

# 파라미터에서 고정점과 수렴속도를 쌓는다
alpha, n, delta, s, A = 1 / 3, 0.01, 0.05, 0.2, 1.0
nd = n + delta                                  # 자본을 깎는 속도 [1/년]
kstar = (s * A / nd) ** (1 / (1 - alpha))       # eq-w10-10
lam = (1 - alpha) * nd                          # eq-w10-10
t_half = np.log(2) / lam                        # eq-w10-9

print(f'k*            = {kstar:.4f} [재화/노동]')
print(f'lambda        = {lam:.4f} [1/년]')
print(f'선형화 반감기 = {t_half:.2f} 년')

def euler(k0, s_, T, dt=0.01):
    """오일러 격자로 dk/dt = s A k^alpha - (n+delta) k 를 한 칸씩 직접 쌓는다."""
    N = int(round(T / dt))
    t = np.arange(N + 1) * dt
    k = np.empty(N + 1)
    k[0] = k0
    for i in range(N):
        k[i + 1] = k[i] + dt * (s_ * A * k[i] ** alpha - nd * k[i])
    return t, k


def crossing(t, k, target):
    """목표 수준을 처음 지나는 두 격자점 사이를 선형 보간한다."""
    j = int(np.argmax(k >= target)) if k[0] < target else int(np.argmax(k <= target))
    return t[j - 1] + (target - k[j - 1]) / (k[j] - k[j - 1]) * (t[j] - t[j - 1])


def exact_time(kap0, kap1):
    """정확해 eq-w10-13을 kappa로 옮겨 간극 절반 도달 시간을 얻는다."""
    return -np.log((kap1 ** (1 - alpha) - 1) / (kap0 ** (1 - alpha) - 1)) / lam

# 항목 1 — 두 출발점에서 간극 절반까지 걸린 연수를 센다
item1 = {}
print(f'{"출발":>8s}{"목표":>8s}{"오일러[년]":>12s}{"정확해[년]":>12s}{"선형화[년]":>12s}{"정확/선형":>11s}')
for kap0, kap1 in [(0.5, 0.75), (2.0, 1.5)]:
    t, k = euler(kap0 * kstar, s, 80.0)
    tg = crossing(t, k, kap1 * kstar)
    tx = exact_time(kap0, kap1)
    item1[kap0] = (tg, tx)
    print(f'{kap0:7.2f}k*{kap1:7.2f}k*{tg:12.2f}{tx:12.2f}{t_half:12.2f}{(tx / t_half - 1) * 100:+10.1f}%')

def ls_slope(x, y):
    """정규방정식 (X'X)b = X'y 를 직접 풀어 1차 최소제곱 기울기를 얻는다."""
    X = np.column_stack([np.ones_like(x), x])
    return float(np.linalg.solve(X.T @ X, X.T @ y)[1])


# 항목 2 — 두 저축률에서 후반 40년의 ln|k-k*| 기울기를 잰다
item2 = {}
k_base = None
print(f'{"s":>6s}{"k*":>10s}{"k* 비":>8s}{"후반40년 기울기":>18s}{"-lambda":>10s}')
for s_ in (0.2, 0.3):
    ks_ = (s_ * A / nd) ** (1 / (1 - alpha))
    t, k = euler(0.5 * ks_, s_, 200.0)
    m = t >= 160.0
    b1 = ls_slope(t[m], np.log(np.abs(k[m] - ks_)))
    k_base = ks_ if k_base is None else k_base
    item2[s_] = (ks_, b1)
    print(f'{s_:6.2f}{ks_:10.4f}{ks_ / k_base:8.2f}{b1:18.4f}{-lam:10.4f}')

def G(k_):
    """연 단위 사상 eq-w10-15를 그대로 쓴다."""
    return (s * A * k_ ** alpha + (1 - delta) * k_) / (1 + n)


# 항목 3 (앞) — 사상을 반복해 간극 절반까지의 연수를 정수와 보간 둘로 센다
path = [0.5 * kstar]
for _ in range(40):
    path.append(G(path[-1]))
path = np.array(path)
years = np.arange(path.size)
j = int(np.argmax(path >= 0.75 * kstar))
t_int = j
t_lin = (j - 1) + (0.75 * kstar - path[j - 1]) / (path[j] - path[j - 1])

print(f'18번째 반복 후 kappa = {path[18] / kstar:.4f}')
print(f'19번째 반복 후 kappa = {path[19] / kstar:.4f}')
print(f'정수 연수            = {t_int} 년')
print(f'교차점 보간          = {t_lin:.2f} 년')
print(f'연속형 정확해        = {exact_time(0.5, 0.75):.2f} 년')

# 항목 3 (뒤) — 고정점에서의 사상 기울기를 중심차분으로 잰다
h = 1e-5
Gp = (G(kstar + h) - G(kstar - h)) / (2 * h)

print(f"G'(k*) 수치미분  = {Gp:.6f}")
print(f"1 - lambda/(1+n) = {1 - lam / (1 + n):.6f}")
print(f"exp(-lambda)     = {np.exp(-lam):.6f}")
print(f"1 - lambda       = {1 - lam:.6f}")
print(f"차 G'(k*) - exp(-lambda) = {Gp - np.exp(-lam):+.6f}")

# 위상선과 고정점에서의 접선을 한 장에 그린다
kk = np.linspace(0.05, 2.3 * kstar, 400)
phi = s * A * kk ** alpha - nd * kk

fig, ax = plt.subplots(figsize=(6.4, 4.0))
ax.plot(kk, phi, color='#1f5fbf', label=r'$\phi(k)=sAk^{\alpha}-(n+\delta)k$')
ax.plot(kk, -lam * (kk - kstar), '--', color='#c0392b', label=r'접선: 기울기 $-\lambda$')
ax.axhline(0.0, color='#222222', lw=0.8)
ax.plot([kstar], [0.0], 'o', color='#222222')
ax.annotate(f'$k^*$ = {kstar:.2f}', (kstar, 0.0), textcoords='offset points', xytext=(8, 10))
ax.set_xlim(0, 2.3 * kstar)
ax.set_ylim(-0.30, 0.30)
ax.set_xlabel('$k$ [재화/노동]')
ax.set_ylabel(r'$\dot k$ [재화/노동/년]')
ax.set_title(f'위상선과 고정점의 기울기 ($-\\lambda$ = {-lam:.2f}/년)')
ax.legend(loc='lower left')
plt.show()

# 간극의 로그가 직선으로 주는 것을 두 출발점과 사상에서 겹쳐 본다
fig, ax = plt.subplots(figsize=(6.4, 4.0))
for kap0, col in [(0.5, '#1f5fbf'), (2.0, '#e67e22')]:
    t, k = euler(kap0 * kstar, s, 80.0)
    ax.plot(t, np.log(np.abs(k / kstar - 1.0)), color=col, label=f'오일러 $\\kappa_0$ = {kap0}')
ax.plot(years, np.log(np.abs(path / kstar - 1.0)), 'o', ms=3, color='#8a8a8a', label='연 단위 사상')
tt = np.linspace(0.0, 80.0, 2)
ax.plot(tt, np.log(0.5) - lam * tt, 'k--', lw=1.0, label=r'기준선 기울기 $-\lambda$')
ax.axvline(t_half, color='#c0392b', ls=':', lw=1.0)
ax.annotate(f'선형화 반감기 {t_half:.1f}년', (t_half, -3.6), textcoords='offset points', xytext=(6, 0))
ax.set_xlim(0, 80)
ax.set_xlabel('$t$ [년]')
ax.set_ylabel(r'$\ln|\kappa-1|$')
ax.set_title('간극의 로그는 기울기 $-\\lambda$로 준다')
ax.legend(loc='upper right')
plt.show()

# 예측 항목 세 개에 대응하는 수치를 한자리에 모은다
g_below, x_below = item1[0.5]
g_above, x_above = item1[2.0]
print('항목 1 · 간극 절반 도달 시간 [년]')
print(f'  아래 출발 0.5k* -> 0.75k* : 오일러 {g_below:.2f} · 정확해 {x_below:.2f} · 선형화 {t_half:.2f}')
print(f'  위   출발 2.0k* -> 1.5k*  : 오일러 {g_above:.2f} · 정확해 {x_above:.2f} · 선형화 {t_half:.2f}')
print('항목 2 · 후반 40년 기울기 [1/년]')
print(f'  s=0.20: k*={item2[0.2][0]:.4f}, 기울기 {item2[0.2][1]:.4f}')
print(f'  s=0.30: k*={item2[0.3][0]:.4f}, 기울기 {item2[0.3][1]:.4f}  (k* 비 {item2[0.3][0] / item2[0.2][0]:.2f}배)')
print('항목 3 · 연 단위 사상')
print(f"  정수 연수 {t_int}년 · 교차점 보간 {t_lin:.2f}년 · 연속형 {x_below:.2f}년")
print(f"  G'(k*) = {Gp:.4f} · exp(-lambda) = {np.exp(-lam):.4f}")