W03 실습 · 듀레이션과 컨벡시티

# 계산 도구와 한국어 글꼴만 쌓는다.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm

CANDIDATES = ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic',
              'Noto Sans CJK KR', 'Malgun Gothic']
AVAILABLE = {f.name for f in fm.fontManager.ttflist}
FONT = next((f for f in CANDIDATES if f in AVAILABLE), 'DejaVu Sans')
plt.rcParams.update({'font.family': FONT, 'font.size': 11, 'axes.unicode_minus': False,
                     'mathtext.fontset': 'cm', 'figure.dpi': 110})
print('한국어 글꼴:', FONT)

# 현금흐름 벡터와 시점을 만들고 할인계수를 직접 쌓아 가격을 낸다 (eq-w03-3).
KAPPA, T0, Y0, FACE = 0.05, 10, 0.05, 100.0


def cashflows(kappa, T, face=FACE):
    t = np.arange(1, T + 1, dtype=float)
    c = np.full(T, kappa * face)
    c[-1] += face                      # 마지막 시점에 액면을 더한다
    return t, c


def price(kappa, T, y):
    t, c = cashflows(kappa, T)
    return float(np.sum(c * (1.0 + y) ** (-t)))


P0 = price(KAPPA, T0, Y0)
print(f'10년 · 연 1회 5% 이표 · y = 5%  →  P = {P0:.6f}  (액면 {FACE:.0f}, 액면가 채권)')

# eq-w03-5의 합으로 가중치 → D_M → D → C 를 쌓는다. 유한차분이 아니다.
def dur_conv(kappa, T, y):
    t, c = cashflows(kappa, T)
    pv = c * (1.0 + y) ** (-t)
    P = float(pv.sum())
    w = pv / P                                              # 현재가치 가중치
    DM = float(np.sum(t * w))
    D = DM / (1.0 + y)
    Cx = float(np.sum(t * (t + 1.0) * w)) / (1.0 + y) ** 2
    K3 = float(np.sum(t * (t + 1.0) * (t + 2.0) * w)) / (1.0 + y) ** 3
    return P, DM, D, Cx, K3


P0, DM0, D0, C0, K30 = dur_conv(KAPPA, T0, Y0)
print(f'D_M = {DM0:.4f} 년   D = {D0:.4f} 년   '
      f'C = {C0:.3f} 년^2   K3^0 = {K30:.2f} 년^3')

# 검산 — 중심차분 [P(y+h)-P(y-h)]/(2hP) 이 -D 와 1e-6 안에서 맞는가 (h = 1e-5).
h = 1e-5
fd = (price(KAPPA, T0, Y0 + h) - price(KAPPA, T0, Y0 - h)) / (2 * h * P0)
print(f'중심차분 {fd:+.8f}   -D {-D0:+.8f}   차이 {abs(fd + D0):.2e}  '
      f'({"통과" if abs(fd + D0) < 1e-6 else "실패"})')

# 같은 가격 함수에 y+dy 를 넣은 정확 재평가와 1차·2차 전개를 나란히 쌓는다 (본문 3절 표).
def expand(dy, T=T0, P=None, D=None, C=None):
    P = P0 if P is None else P
    D = D0 if D is None else D
    C = C0 if C is None else C
    exact = price(KAPPA, T, Y0 + dy) / P - 1.0            # 정확 재평가
    first = -D * dy                                       # 1차 (듀레이션)
    second = first + 0.5 * C * dy ** 2                    # 2차 (컨벡시티까지)
    return exact, first, second, exact - first, exact - second


GRID = [1e-4, 1e-3, 5e-3, 1e-2, 2e-2, 4e-2, -4e-2, 5e-2]
head = f'{"dy(bp)":>8}{"정확(%)":>11}{"1차(%)":>11}{"2차(%)":>11}'
print(head + f'{"E1(%p)":>12}{"E2(%p)":>13}{"rho":>9}')
for dy in GRID:
    ex, o1, o2, E1, E2 = expand(dy)
    rho = C0 * abs(dy) / (2.0 * D0)
    print(f'{dy * 1e4:+8.0f}{ex * 100:+11.4f}{o1 * 100:+11.4f}{o2 * 100:+11.4f}'
          f'{E1 * 100:+12.6f}{E2 * 100:+13.7f}{rho * 100:8.2f}%')

# 예측 1 — +50bp 한 줄만 따로 뽑는다.
ex50, o1_50, o2_50, E1_50, E2_50 = expand(0.005)
sign = '양(+)' if E1_50 > 0 else '음(-)'
print(f'[예측 1] +50bp  정확 {ex50 * 100:+.4f}%   1차 {o1_50 * 100:+.4f}%   '
      f'2차 {o2_50 * 100:+.4f}%')
print(f'[예측 1] E1 = 정확 - 1차 = {E1_50 * 100:+.4f}%p   부호 {sign}')
print(f'[예측 1] 1차 근사는 손실을 {abs(E1_50) * 100:.4f}%p 과대평가한다 — 접선이 가격 아래에 있다.')

# 예측 2 — 만기 T=1..30 을 돌며 +50bp 의 E1 과 컨벡시티 항을 쌓는다.
DY = 0.005
Ts = np.arange(1, 31)
E1_T = np.empty(len(Ts))
halfC_T = np.empty(len(Ts))
for i, T in enumerate(Ts):
    P, DM, D, Cx, _ = dur_conv(KAPPA, int(T), Y0)
    ex = price(KAPPA, int(T), Y0 + DY) / P - 1.0
    E1_T[i] = ex + D * DY
    halfC_T[i] = 0.5 * Cx * DY ** 2

print(f'{"T(년)":>6}{"E1(+50bp, %p)":>16}{"(1/2)C dy^2 (%p)":>20}{"C(년^2)":>12}')
for T in (1, 5, 10, 20, 30):
    i = T - 1
    Cx = dur_conv(KAPPA, int(T), Y0)[3]
    print(f'{T:6d}{E1_T[i] * 100:+16.4f}{halfC_T[i] * 100:+20.4f}{Cx:12.2f}')
print(f'[예측 2] 30년 / 10년 = {E1_T[29] / E1_T[9]:.2f}배, '
      f'10년 / 1년 = {E1_T[9] / E1_T[0]:.2f}배 — 만기가 길수록 커진다.')

# 예측 3 — 1bp 와 100bp 의 오차를 견주어 배율을 쌓는다.
_, _, _, E1_1, E2_1 = expand(1e-4)
_, _, _, E1_100, E2_100 = expand(1e-2)
print(f'[예측 3] E1: {E1_1 * 100:.4e}%p → {E1_100 * 100:.6f}%p   '
      f'배율 {E1_100 / E1_1:,.0f}   (10^4 = 10,000)')
print(f'[예측 3] E2: {E2_1 * 100:.4e}%p → {E2_100 * 100:.6f}%p   '
      f'배율 {E2_100 / E2_1:,.0f}   (10^6 = 1,000,000)')

# |dy| 를 로그 간격으로 훑어 오차를 쌓고 로그-로그 기울기로 차수를 재다.
dys = np.logspace(-4, np.log10(0.05), 40)
E1s = np.array([expand(d)[3] for d in dys])
E2s = np.array([expand(d)[4] for d in dys])
slope1 = float(np.polyfit(np.log(dys), np.log(np.abs(E1s)), 1)[0])
slope2 = float(np.polyfit(np.log(dys), np.log(np.abs(E2s)), 1)[0])
print(f'[예측 3] 로그-로그 기울기: |E1| {slope1:.3f} (2 기대)   '
      f'|E2| {slope2:.3f} (3 기대)')

# 그림 1 — 오차의 크기 순위를 로그-로그로 쌓는다.
from matplotlib.ticker import FuncFormatter, NullFormatter

fig, ax = plt.subplots(figsize=(6.6, 4.2))
ax.loglog(dys * 1e4, np.abs(E1s), color='#1f5fbf', label='$|E_1|$ 1차 근사 오차')
ax.loglog(dys * 1e4, np.abs(E2s), color='#c0392b', label='$|E_2|$ 2차 근사 오차')
ax.loglog(dys * 1e4, np.abs(E1s[0]) * (dys / dys[0]) ** 2,
          color='#8a8a8a', lw=0.9, ls='--', label='기울기 2')
ax.loglog(dys * 1e4, np.abs(E2s[0]) * (dys / dys[0]) ** 3,
          color='#8a8a8a', lw=0.9, ls=':', label='기울기 3')
# 눈금은 직접 적는다 — 기본 눈금표의 유니코드 마이너스가 한국어 글꼴에 없다.
ax.xaxis.set_major_formatter(FuncFormatter(lambda v, _: f'{v:g}'))
ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: '$10^{%d}$' % round(np.log10(v))))
ax.xaxis.set_minor_formatter(NullFormatter())
ax.yaxis.set_minor_formatter(NullFormatter())
ax.set_xlabel('$|\\Delta y|$ (bp)')
ax.set_ylabel('오차 (가격 대비, 절댓값)')
ax.set_title(f'추정 기울기 {slope1:.2f} 와 {slope2:.2f} — 2차와 3차')
ax.legend(loc='upper left')
plt.show()

# 그림 2 — +50bp 오차가 만기를 따라 자라는 모습을 쌓는다.
fig, ax = plt.subplots(figsize=(6.6, 4.2))
ax.plot(Ts, halfC_T * 100, color='#1f5fbf', label='$\\frac{1}{2}C(T)(\\Delta y)^2$')
ax.plot(Ts, E1_T * 100, 'o', ms=4, color='#222222', label='$E_1(+50\\,\\mathrm{bp})$')
for T, dx, dy_ in ((1, 0.6, 0.050), (10, 0.9, -0.033), (30, 0.7, -0.033)):
    ax.annotate(f'{E1_T[T - 1] * 100:.3f}%p', xy=(T, E1_T[T - 1] * 100),
                xytext=(T + dx, E1_T[T - 1] * 100 + dy_), fontsize=9)
ax.set_xlim(0, 34)
ax.set_xlabel('만기 $T$ (년)')
ax.set_ylabel('1차 근사 오차 (%p)')
ax.set_title('$+50$bp에서 만기별 $E_1$')
ax.legend(loc='upper left')
plt.show()