W12 실습 · 벨만 방정식 — 연산자와 그 고정점

# 도구와 한국어 글꼴, 그리고 이 회차의 파라미터를 쌓는다
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager

CAND = ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic', 'Noto Sans CJK KR', 'Malgun Gothic']
AVAIL = {f.name for f in font_manager.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'

ALPHA = 0.4      # 자본분배율
BETA = 0.96      # 할인인자 = 압축계수
TOL = 1e-6       # 정지 규칙
NGRID = 300      # 격자점 수
print('글꼴:', FONT)
print(f'alpha={ALPHA}, beta={BETA}, 격자 {NGRID}점, 정지 기준 {TOL:g}')

# 격자 위 소비 c = k_i^alpha - k_j 와 그 효용을 행렬로 쌓는다 (delta=1)
def build(n, lo=0.01, hi=1.0, alpha=ALPHA):
    k = np.linspace(lo, hi, n)
    y = k ** alpha                  # 처분가능 자원 y(k) = k^alpha
    c = y[:, None] - k[None, :]     # (i,j)칸의 소비
    u = np.where(c > 0, np.log(np.where(c > 0, c, 1.0)), -np.inf)
    return k, u

k300, u300 = build(NGRID)
n_feas = np.isfinite(u300).sum(axis=1)
print(f'격자 간격 = {k300[1] - k300[0]:.5f}')
print(f'실현가능한 선택 수: k=0.010에서 {n_feas[0]}개, k=1.000에서 {n_feas[-1]}개')
print(f'효용의 범위: {u300[np.isfinite(u300)].min():.3f} ~ {u300[np.isfinite(u300)].max():.3f} (유한상태이므로 유계)')

# 벨만 연산자 한 번과 그 반복을 직접 쌓는다
def bellman(V, u, beta):
    M = u + beta * V[None, :]       # (i,j)칸의 후보값
    return M.max(axis=1), M.argmax(axis=1)

def vfi(u, beta, tol=TOL, maxit=20000):
    V = np.zeros(u.shape[0])        # V_0 = 0
    gaps = []
    for _ in range(maxit):
        Vn, pol = bellman(V, u, beta)
        gaps.append(float(np.max(np.abs(Vn - V))))   # 상한 노름 갱신량
        V = Vn
        if gaps[-1] < tol:
            break
    return V, pol, np.array(gaps)

_, _, g1 = vfi(u300, BETA, tol=np.inf, maxit=1)     # 첫 갱신 하나만 본다
n_apriori = np.log(TOL * (1 - BETA) / g1[0]) / np.log(BETA)
print(f'한 번 갱신한 뒤 ||TV_0 - V_0|| = {g1[0]:.4f}')
print(f'사전 상계 ||V_0 - V|| <= ||TV_0 - V_0||/(1-beta) = {g1[0] / (1 - BETA):.1f}')
print(f'오차를 {TOL:g} 아래로 보증하는 n = ln(TOL(1-beta)/||TV_0-V_0||)/ln(beta) '
      f'= {n_apriori:.1f} -> {np.ceil(n_apriori):.0f}회')

# beta=0.96에서 갱신량의 비·반감기·사후 상계를 쌓는다 (예측 1)
V96, pol96, gaps96 = vfi(u300, BETA)
ratio96 = gaps96[1:] / gaps96[:-1]
half96 = np.log(2) / -np.log(BETA)
post96 = BETA / (1 - BETA) * gaps96[-1]
print(f'[예측 1] 정지까지 반복 = {len(gaps96)}회')
print(f'[예측 1] 갱신비: n=1 {ratio96[0]:.4f} | n=5 {ratio96[4]:.4f} | '
      f'n=10 {ratio96[9]:.4f} | n=100 {ratio96[99]:.4f} | 마지막 {ratio96[-1]:.4f}')
print(f'[예측 1] 오차 반감기 = ln2/(-ln beta) = {half96:.2f}회')
print(f'[예측 1] 마지막 갱신량 {gaps96[-1]:.3e}, 사후 상계 계수 {BETA / (1 - BETA):.0f} '
      f'-> 진짜 오차 <= {post96:.2e}')
print(f'[예측 1] 사전 상계가 요구하는 n = {np.ceil(n_apriori):.0f}회 vs 실제 정지 {len(gaps96)}회 '
      f'— 첫 갱신 하나로 미리 잡은 횟수가 보수적으로 크다')

# 같은 격자에서 beta만 바꿔 반복 횟수를 다시 센다 (예측 2)
BETA2 = 0.90
V90, pol90, gaps90 = vfi(u300, BETA2)
half90 = np.log(2) / -np.log(BETA2)
print(f'[예측 2] beta=0.90 정지까지 반복 = {len(gaps90)}회 (beta=0.96은 {len(gaps96)}회)')
print(f'[예측 2] 횟수 비 = {len(gaps90) / len(gaps96):.3f}, '
      f'이론값 ln0.96/ln0.90 = {np.log(0.96) / np.log(0.90):.3f}')
print(f'[예측 2] 반감기 = {half90:.2f}회, 마지막 갱신비 = {gaps90[-1] / gaps90[-2]:.4f}')

# 닫힌 해 (eq-w12-17)을 쌓아 격자 해와 맞댄다 (예측 3-i, 3-ii)
def closed(k, beta, alpha=ALPHA):
    B = alpha / (1 - alpha * beta)
    A = (np.log(1 - alpha * beta)
         + alpha * beta / (1 - alpha * beta) * np.log(alpha * beta)) / (1 - beta)
    return A + B * np.log(k), (1 - alpha * beta) * k ** alpha, A, B

Vc, cc, A, B = closed(k300, BETA)
c96 = k300 ** ALPHA - k300[pol96]        # 격자 정책이 고른 소비
rel96 = np.abs(c96 - cc) / cc
kstar = (ALPHA * BETA) ** (1 / (1 - ALPHA))
print(f'[예측 3-i] 정책 최대 상대오차 = {rel96.max() * 100:.2f}% (격자 간격 {k300[1] - k300[0]:.5f})')
print(f'[예측 3-ii] max|V_n - V| = {np.max(np.abs(V96 - Vc)):.2e}  '
      f'(A={A:.2f}, B={B:.4f}, ||V||_inf={np.abs(Vc).max():.2f})')
print(f'정상상태 k* = {kstar:.4f}, 닫힌 해의 c*(k*) = {(1 - ALPHA * BETA) * kstar ** ALPHA:.4f}')

# 유한차분 V'와 u'(c) g_k 의 비를 쌓는다 (예측 3-iii)
h = k300[1] - k300[0]
dV = np.empty_like(V96)                            # 유한차분 V'를 직접 쌓는다
dV[1:-1] = (V96[2:] - V96[:-2]) / (k300[2:] - k300[:-2])   # 내부는 중심차분, 오차 O(h^2)
dV[0] = (V96[1] - V96[0]) / h                      # 양 끝은 한쪽 차분이라 정확도가
dV[-1] = (V96[-1] - V96[-2]) / h                   # 한 차수 낮다 — 오차 O(h)
env = (1.0 / c96) * ALPHA * k300 ** (ALPHA - 1)    # u'(c)*g_k, delta=1이므로 g_k=f'(k)
ratio = dV / env
inner = slice(5, -5)
print(f'[예측 3-iii] 내부(양 끝 5점 제외) 비 범위 = '
      f'{ratio[inner].min():.4f} ~ {ratio[inner].max():.4f}')
print(f'[예측 3-iii] 내부에서 |비-1|<0.01인 점의 비율 = {np.mean(np.abs(ratio[inner] - 1) < 0.01) * 100:.1f}%')
print(f'[예측 3-iii] 양 끝: k=0.010에서 {ratio[0]:.3f}, k=1.000에서 {ratio[-1]:.3f}')
off = np.flatnonzero(np.abs(ratio - 1) > 0.01)
print(f'[예측 3-iii] 격자 전체에서 |비-1|>0.01인 점 = {off.size}개, '
      f'k = {k300[off].min():.3f}~{k300[off].max():.3f}')
i_s = int(np.argmin(np.abs(k300 - kstar)))
print(f'정상상태 근처 k={k300[i_s]:.4f}: V\'={dV[i_s]:.3f}, u\'(c)g_k={env[i_s]:.3f}, '
      f'닫힌 해 B/k={B / k300[i_s]:.3f}')

# 격자를 600점으로 늘려 격자 오차와 반복 오차를 갈라 본다
k600, u600 = build(600)
V600, pol600, gaps600 = vfi(u600, BETA)
Vc600, cc600, _, _ = closed(k600, BETA)
c600 = k600 ** ALPHA - k600[pol600]
print('격자 | 반복 | max|V_n-V| | 정책 최대 상대오차 | 마지막 갱신량')
print(f'{300:4d} | {len(gaps96):4d} | {np.max(np.abs(V96 - Vc)):10.2e} | '
      f'{np.max(rel96) * 100:16.2f}% | {gaps96[-1]:.2e}')
print(f'{600:4d} | {len(gaps600):4d} | {np.max(np.abs(V600 - Vc600)):10.2e} | '
      f'{np.max(np.abs(c600 - cc600) / cc600) * 100:16.2f}% | {gaps600[-1]:.2e}')
print('반복 횟수는 격자가 아니라 beta의 것이고, 격자를 늘리면 격자 오차만 준다')

# 갱신량의 감쇠를 semi-log로 그린다
fig, ax = plt.subplots(figsize=(7.2, 4.2))
ax.semilogy(np.arange(1, len(gaps96) + 1), gaps96, color='C0', label=r'$\beta=0.96$')
ax.semilogy(np.arange(1, len(gaps90) + 1), gaps90, color='C1', label=r'$\beta=0.90$')
ax.axvline(half96, color='C0', ls=':', lw=1)
ax.axvline(half90, color='C1', ls=':', lw=1)
ax.axhline(TOL, color='gray', ls='--', lw=0.8)
ax.set_yticks([1e0, 1e-2, 1e-4, 1e-6])
ax.set_yticklabels([r'$10^{0}$', r'$10^{-2}$', r'$10^{-4}$', r'$10^{-6}$'])
ax.text(half96 + 4, 3e-6, f'반감기 {half96:.1f}회', color='C0', fontsize=9)
ax.text(half90 + 4, 1e-6, f'반감기 {half90:.1f}회', color='C1', fontsize=9)
ax.set_xlabel('반복 횟수 n')
ax.set_ylabel(r'$\|V_{n+1}-V_n\|_\infty$')
ax.set_title(r'갱신량의 기하급수 감쇠 — semi-log 기울기는 $\ln\beta$')
ax.legend()
fig.tight_layout()
plt.show()

# 가치함수·정책·포락선 비를 한 줄에 그린다
fig, axes = plt.subplots(1, 3, figsize=(12.6, 3.8))
axes[0].plot(k300, V96, color='C0', lw=2, label=r'격자 해 $V_n$')
axes[0].plot(k300, Vc, 'r--', lw=1, label=r'닫힌 해 $A+B\ln k$')
axes[0].set_xlabel('$k$')
axes[0].set_ylabel('가치 [효용]')
axes[0].legend()
axes[1].plot(k300, c96, color='C0', lw=2, label=r'격자 정책 $c_n$')
axes[1].plot(k300, cc, 'r--', lw=1, label=r'$(1-\alpha\beta)k^{\alpha}$')
axes[1].set_xlabel('$k$')
axes[1].set_ylabel('소비 [수량]')
axes[1].legend()
axes[2].plot(k300, ratio, color='C0', lw=1)
axes[2].axhline(1.0, color='r', ls='--', lw=1, label='비 = 1')
axes[2].axhspan(0.99, 1.01, color='gray', alpha=0.2, label=r'$\pm1\%$')
axes[2].set_ylim(0.85, 1.06)
axes[2].set_xlabel('$k$')
axes[2].set_ylabel(r"$V'/(u'(c)\,g_k)$")
axes[2].set_title('포락선 비', fontsize=10)
axes[2].legend(fontsize=8)
fig.tight_layout()
plt.show()

# 예측 3항목에 대응하는 수치를 한 곳에 모은다
print(f'예측 1 | 갱신비 -> {ratio96[-1]:.4f} | 반감기 {half96:.2f}회 | '
      f'정지 {len(gaps96)}회 | 진짜 오차 <= {post96:.2e}')
print(f'예측 2 | beta=0.90 정지 {len(gaps90)}회 | 횟수 비 {len(gaps90) / len(gaps96):.3f} | '
      f'반감기 {half90:.2f}회')
print(f'예측 3 | 정책 최대 상대오차 {rel96.max() * 100:.2f}% | '
      f'max|V_n-V| {np.max(np.abs(V96 - Vc)):.2e} | '
      f'내부 포락선 비 {ratio[inner].min():.3f}~{ratio[inner].max():.3f} | '
      f'양 끝 왼쪽 {ratio[0]:.3f} / 오른쪽 {ratio[-1]:.3f}')