노트북은 예측 → 계산 → 대조 세 부분으로 고정한다. 새 개념은 도입하지 않는다. 본문에서 이미 유도한 것을 수치로 확인할 뿐이다.
1. 예측 (코드를 쓰기 전에)¶
아래 셋을 먼저 채운다. 계산하지 않는다.
기준 파라미터에서 , , , 반감기는 얼마인가. 같은 의 솔로우(W10, )보다 빠른가 느린가.
예측: ____
에서 을 안장경로 값의 로 잡고 정방향 적분하면 각각 어디로 가는가. 안장경로 값 그대로 정방향으로 쏘면 상대오차 10-6이 언제 이 되는가.
예측: ____
를 로 영구히 내리면 옛 정상상태에서 소비는 어느 방향으로 얼마나 뛰는가.
예측: ____
2. 계산¶
패키지 호출로 답을 내지 않는다. 정상상태·야코비안·적분기·역행 사격을 차례로 직접 쌓는다. 결과는 표와 그림 둘로 낸다.
# 기준 파라미터와 생산함수를 쌓는다.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import font_manager
from matplotlib.ticker import FixedLocator, FuncFormatter, LogLocator, NullFormatter, NullLocator
# 한국어 글꼴 설정
_have = {fo.name for fo in font_manager.fontManager.ttflist}
for _cand in ['Apple SD Gothic Neo', 'AppleGothic', 'NanumGothic',
'Noto Sans CJK KR', 'Malgun Gothic']:
if _cand in _have:
plt.rcParams['font.family'] = _cand
break
plt.rcParams['axes.unicode_minus'] = False
# 로그축 눈금은 10의 거듭제곱으로 직접 적는다. 기본 로그 포매터는 U+2212(−)를 쓰는데
# 위에서 고른 한글 글꼴에 그 글리프가 없어 10^-2 의 부호가 네모로 깨진다.
pow10 = FuncFormatter(lambda v, _: '$10^{%d}$' % round(np.log10(v)))
ALPHA, RHO, DELTA, SIGMA = 1 / 3, 0.03, 0.05, 2.0 # 연 단위 기준 파라미터
def f(k): # 집약형 생산함수
return k ** ALPHA
def fp(k): # 자본의 한계생산
return ALPHA * k ** (ALPHA - 1)
def fpp(k): # 한계생산의 도함수
return ALPHA * (ALPHA - 1) * k ** (ALPHA - 2)
def phi(k): # 유지가능 소비 — 자본이 멈추는 곡선
return f(k) - DELTA * k
print('alpha=%.4f rho=%.2f delta=%.2f sigma=%.1f' % (ALPHA, RHO, DELTA, SIGMA))alpha=0.3333 rho=0.03 delta=0.05 sigma=2.0
# 정상상태를 두 길로 구해 대조한다 — 닫힌 형과 뉴턴 반복.
def steady_closed(rho=RHO):
ks = (ALPHA / (rho + DELTA)) ** (1 / (1 - ALPHA))
return ks, phi(ks)
def steady_newton(rho=RHO, k_init=1.0, tol=1e-13, itmax=60):
k, steps = k_init, 0
for _ in range(itmax):
step = (fp(k) - rho - DELTA) / fpp(k) # g(k)=f'(k)-rho-delta 의 뉴턴 걸음
k, steps = k - step, steps + 1
if abs(step) < tol:
break
return k, phi(k), steps
KS, CS = steady_closed()
kn, cn, nstep = steady_newton()
K_GOLD = (ALPHA / DELTA) ** (1 / (1 - ALPHA))
K_BAR = DELTA ** (-1 / (1 - ALPHA))
S_STAR = DELTA * KS / f(KS)
print('닫힌 형 k*=%.6f c*=%.6f' % (KS, CS))
print('뉴턴 반복 k*=%.6f c*=%.6f (반복 %d회)' % (kn, cn, nstep))
print('두 길의 차이 |dk*|=%.2e' % abs(KS - kn))
print('k_gold=%.3f k_bar=%.3f y*=%.4f s*=delta k*/y*=%.4f (1/s*=%.2f)'
% (K_GOLD, K_BAR, f(KS), S_STAR, 1 / S_STAR))닫힌 형 k*=8.505173 c*=1.615983
뉴턴 반복 k*=8.505173 c*=1.615983 (반복 9회)
두 길의 차이 |dk*|=1.78e-15
k_gold=17.213 k_bar=89.443 y*=2.0412 s*=delta k*/y*=0.2083 (1/s*=4.80)
# 야코비안과 고유값을 (eq-w11-13)·(eq-w11-14) 공식으로 손으로 쌓는다.
def jacobian(ks, cs, sigma=SIGMA, rho=RHO):
return np.array([[fp(ks) - DELTA, -1.0],
[cs * fpp(ks) / sigma, (fp(ks) - rho - DELTA) / sigma]])
def eigen_hand(ks, cs, sigma=SIGMA, rho=RHO):
tr, det = rho, cs * fpp(ks) / sigma
disc = np.sqrt(tr ** 2 - 4 * det)
return tr, det, (tr - disc) / 2, (tr + disc) / 2
J = jacobian(KS, CS)
TR, DET, LAM1, LAM2 = eigen_hand(KS, CS)
chk = np.sort(np.linalg.eigvals(J).real) # np.linalg.eig 는 검산에만
HALF = np.log(2) / abs(LAM1)
HALF_SOLOW = np.log(2) / ((1 - ALPHA) * DELTA) # W10 솔로우(n=0, s 고정)
print('J =\n%s' % np.array2string(J, precision=5))
print('trJ=%.5f detJ=%.6f 판별식=%.6f f\'\'(k*)=%.6f' % (TR, DET, TR ** 2 - 4 * DET, fpp(KS)))
print('손 계산 lambda1=%.5f lambda2=%.5f (lambda1+lambda2=%.5f=rho)' % (LAM1, LAM2, LAM1 + LAM2))
print('eig 검산 lambda1=%.5f lambda2=%.5f 차이 %.1e' % (chk[0], chk[1], abs(chk[0] - LAM1)))
print('[예측 1] k*=%.3f c*=%.3f lambda1=%.4f 반감기=%.2f년'
% (KS, CS, LAM1, HALF))
print('[예측 1] 솔로우 반감기=%.2f년 -> 램지가 %s' % (HALF_SOLOW, '빠르다' if HALF < HALF_SOLOW else '느리다'))J =
[[ 3.00000e-02 -1.00000e+00]
[-5.06667e-03 6.93889e-18]]
trJ=0.03000 detJ=-0.005067 판별식=0.021167 f''(k*)=-0.006271
손 계산 lambda1=-0.05774 lambda2=0.08774 (lambda1+lambda2=0.03000=rho)
eig 검산 lambda1=-0.05774 lambda2=0.08774 차이 6.9e-18
[예측 1] k*=8.505 c*=1.616 lambda1=-0.0577 반감기=12.00년
[예측 1] 솔로우 반감기=20.79년 -> 램지가 빠르다
# RK4 적분기를 직접 쓴다 — 벡터장은 (eq-w11-1).
def field(x, rho=RHO, sigma=SIGMA):
k, c = x
return np.array([f(k) - c - DELTA * k, c * (fp(k) - rho - DELTA) / sigma])
def rk4_step(x, h, rho=RHO, sigma=SIGMA):
a = field(x, rho, sigma)
b = field(x + 0.5 * h * a, rho, sigma)
c = field(x + 0.5 * h * b, rho, sigma)
d = field(x + h * c, rho, sigma)
return x + h / 6 * (a + 2 * b + 2 * c + d)
def integrate(x0, h, tmax, rho=RHO, sigma=SIGMA, stop=None):
x = np.array(x0, float) # h<0 이면 시간을 거꾸로 간다
ts, xs, t = [0.0], [x.copy()], 0.0
while t < tmax - 1e-12:
x = rk4_step(x, h, rho, sigma)
t += abs(h)
if not np.all(np.isfinite(x)) or x[0] <= 0 or x[1] <= 0:
break
ts.append(t)
xs.append(x.copy())
if stop is not None and stop(x):
break
return np.array(ts), np.array(xs)
print('RK4 한 걸음 검산: field(k*,c*)=%s' % np.array2string(field([KS, CS]), precision=8))RK4 한 걸음 검산: field(k*,c*)=[5.55111512e-17 1.12131333e-17]
# 역행 사격으로 안장경로를 쌓는다 — 시간을 뒤집으면 안정 다양체가 불안정 다양체가 된다.
def arm_branches(rho=RHO, sigma=SIGMA, eps=1e-3, h=0.05, kmin=0.3, kmax=40.0):
"""왼쪽·오른쪽 두 가지를 각각 (뒤로 간 시간, 궤적) 으로 돌려준다."""
ks, cs = steady_closed(rho)
lam2 = eigen_hand(ks, cs, sigma, rho)[3]
v = np.array([1.0, lam2]) / np.hypot(1.0, lam2) # 안정 고유벡터 방향 (eq-w11-15)
return [integrate(np.array([ks, cs]) + s * eps * v, -h, 4000.0, rho, sigma,
stop=lambda x: not (kmin < x[0] < kmax))
for s in (-1.0, +1.0)]
def arm_table(branches):
"""두 가지를 k 순으로 이어 붙여 표 c_arm(k) 로 만든다."""
xs = np.vstack([b[1] for b in branches])
order = np.argsort(xs[:, 0])
return xs[order, 0], xs[order, 1]
ARM = arm_branches()
KA, CA = arm_table(ARM)
K0 = 3.0
C_ARM0 = float(np.interp(K0, KA, CA))
C_LIN0 = CS + LAM2 * (K0 - KS)
print('안장경로 표: k in [%.2f, %.2f], 점 %d개' % (KA[0], KA[-1], KA.size))
print('k0=%.0f 에서 c_arm=%.4f 선형 근사 c*+lambda2(k-k*)=%.4f 차이 %+.1f%%'
% (K0, C_ARM0, C_LIN0, 100 * (C_LIN0 / C_ARM0 - 1)))
print('안장경로 위 저축률 1-c_arm/f(k0)=%.4f > 정상상태 s*=%.4f' % (1 - C_ARM0 / f(K0), S_STAR))안장경로 표: k in [0.29, 40.04], 점 6728개
k0=3 에서 c_arm=1.0142 선형 근사 c*+lambda2(k-k*)=1.1329 차이 +11.7%
안장경로 위 저축률 1-c_arm/f(k0)=0.2968 > 정상상태 s*=0.2083
# 그림 1 — 두 널클라인, 안장경로, k0=3 에서 ±3% 로 출발한 두 경로.
up = integrate([K0, C_ARM0 * 1.03], 0.05, 400.0, stop=lambda x: x[0] < 0.05)
dn = integrate([K0, C_ARM0 * 0.97], 0.05, 600.0)
kk = np.linspace(0.05, 30.0, 400)
fig, ax = plt.subplots(figsize=(6.6, 4.3))
ax.plot(kk, phi(kk), color='0.45', lw=1.2, label=r'$\dot k=0$ (유지가능 소비)')
ax.axvline(KS, color='0.45', lw=1.0, ls='--', label=r'$\dot c=0$ ($k=k^*$)')
ax.plot(KA, CA, color='C0', lw=2.2, label='안장경로 (역행 사격)')
ax.plot(up[1][:, 0], up[1][:, 1], color='C3', lw=1.4, label='+3% 경로')
ax.plot(dn[1][:, 0], dn[1][:, 1], color='C1', lw=1.4, label='-3% 경로')
ax.plot([KS], [CS], 'ko', ms=5)
ax.plot([K0, K0], [C_ARM0 * 0.97, C_ARM0 * 1.03], color='0.2', lw=0.9, ls=':')
ax.set_xlim(0, 30)
ax.set_ylim(0, 2.6)
ax.set_xlabel('1인당 자본 k')
ax.set_ylabel('1인당 소비 c')
ax.set_title('램지 위상평면 — 안장경로와 두 발산 경로')
ax.legend(fontsize=8, loc='lower right')
plt.show()
# [예측 2] 세 경로의 종착지와 횡단조건 양 e^{-rho t} mu_t k_t 를 쌓는다.
def tvc(ts, xs, rho=RHO, sigma=SIGMA):
return np.exp(-rho * ts) * xs[:, 1] ** (-sigma) * xs[:, 0] # mu = c^{-sigma}
def arm_path(branches, k0, ks=KS):
"""안장경로의 시간 경로. 안장경로 위에서 정방향으로 쏘면 다음 셀이 재는 속도로
다양체를 이탈하므로, 역행 사격 궤적에서 k0 쪽 구간을 잘라 t -> t_max - t 로
시간만 뒤집어 쓴다."""
ts, xs = branches[0] if k0 < ks else branches[1]
m = (xs[:, 0] >= k0) if k0 < ks else (xs[:, 0] <= k0)
return (ts[m].max() - ts[m])[::-1], xs[m][::-1]
arm = arm_path(ARM, K0)
print('경로 종료 t(년) k_T c_T e^{-rho t}mu k (T)')
for name, (ts, xs) in [('안장경로 ', arm), ('+3% ', up), ('-3% ', dn)]:
print('%s %8.1f %10.4f %10.5f %14.3e' % (name, ts[-1], xs[-1, 0], xs[-1, 1], tvc(ts, xs)[-1]))
print(' (안장경로 행의 종료 t 는 역행 사격이 k* 의 1e-3 근방에 닿는 시점 — k_T -> k*=%.4f)' % KS)
print()
print('[예측 2] +3%%: k=0 도달 t=%.1f년 (실행 불가능)' % up[0][-1])
print('[예측 2] -3%%: (k,c) -> (%.3f, %.5f), k_bar=%.3f, 횡단조건 양 %.3e 로 발산'
% (dn[1][-1, 0], dn[1][-1, 1], K_BAR, tvc(dn[0], dn[1])[-1]))
print(' -3%% 경로의 발산 기울기 (1-alpha)delta=%.4f, 안장경로는 기울기 -rho=%.4f 로 0 에 간다'
% ((1 - ALPHA) * DELTA, -RHO))
late = arm[0] > 100
print(' 검산: 안장경로 후반(t>100년) 로그 기울기 적합 %.5f vs -rho=%.5f'
% (np.polyfit(arm[0][late], np.log(tvc(*arm)[late]), 1)[0], -RHO))경로 종료 t(년) k_T c_T e^{-rho t}mu k (T)
안장경로 150.9 8.5042 1.61590 3.527e-02
+3% 24.9 0.0498 2.52824 3.694e-03
-3% 600.0 89.4427 0.00000 1.163e+09
(안장경로 행의 종료 t 는 역행 사격이 k* 의 1e-3 근방에 닿는 시점 — k_T -> k*=8.5052)
[예측 2] +3%: k=0 도달 t=24.9년 (실행 불가능)
[예측 2] -3%: (k,c) -> (89.443, 0.00000), k_bar=89.443, 횡단조건 양 1.163e+09 로 발산
-3% 경로의 발산 기울기 (1-alpha)delta=0.0333, 안장경로는 기울기 -rho=-0.0300 로 0 에 간다
검산: 안장경로 후반(t>100년) 로그 기울기 적합 -0.03000 vs -rho=-0.03000
# [예측 2] 정방향 사격 — 상대오차 1e-6 이 언제 O(1) 이 되는가.
EPS = 1e-6
t_b, x_b = integrate([K0, C_ARM0], 0.01, 400.0, stop=lambda x: x[0] < 0.05)
t_p, x_p = integrate([K0, C_ARM0 * (1 + EPS)], 0.01, 400.0, stop=lambda x: x[0] < 0.05)
n = min(t_b.size, t_p.size)
dev = np.abs(x_p[:n, 1] - x_b[:n, 1]) / x_b[:n, 1]
band = (t_b[:n] > 60) & (t_b[:n] < 130)
rate = np.polyfit(t_b[:n][band], np.log(dev[band]), 1)[0]
print('상대오차 dev(t) 표')
for lvl in (1e-6, 1e-4, 1e-2, 1e-1):
idx = int(np.argmax(dev >= lvl))
print(' dev >= %7.0e 에서 t=%6.1f년' % (lvl, t_b[idx]))
t10 = t_b[int(np.argmax(dev >= 0.1))]
print('[예측 2] 외삽한 dev=1 도달 t=%.1f년 (t(dev=0.1)=%.1f 에 ln10/lambda2 를 더한다),'
' 이론값 ln(1/eps)/lambda2=%.1f년'
% (t10 + np.log(10) / LAM2, t10, np.log(1 / EPS) / LAM2))
print(' 검산: 증폭률 적합(60~130년) %.5f vs 손 계산 lambda2 %.5f (외삽값 %.1f년)'
% (rate, LAM2, t10 + np.log(10) / rate))
print(' 정방향 사격은 이 속도로 다양체를 이탈한다 — 그래서 그림 2 의 안장경로는'
' 역행 사격 궤적을 시간만 뒤집어 그린다.')상대오차 dev(t) 표
dev >= 1e-06 에서 t= 0.0년
dev >= 1e-04 에서 t= 55.2년
dev >= 1e-02 에서 t= 107.4년
dev >= 1e-01 에서 t= 132.9년
[예측 2] 외삽한 dev=1 도달 t=159.2년 (t(dev=0.1)=132.9 에 ln10/lambda2 를 더한다), 이론값 ln(1/eps)/lambda2=157.5년
검산: 증폭률 적합(60~130년) 0.08830 vs 손 계산 lambda2 0.08774 (외삽값 159.0년)
정방향 사격은 이 속도로 다양체를 이탈한다 — 그래서 그림 2 의 안장경로는 역행 사격 궤적을 시간만 뒤집어 그린다.
# 그림 2 — 세 경로의 소비와 횡단조건 양. 지평 T=150년까지 잘라 그린다.
T_FIG = 150.0
fig, (a1, a2) = plt.subplots(2, 1, figsize=(6.6, 5.4), sharex=True)
for name, (ts, xs), col in [('안장경로', arm, 'C0'), ('+3%', up, 'C3'), ('-3%', dn, 'C1')]:
m = ts <= T_FIG
a1.plot(ts[m], xs[m, 1], color=col, lw=1.6, label=name)
a2.semilogy(ts[m], tvc(ts, xs)[m], color=col, lw=1.6, label=name)
a1.axhline(CS, color='0.45', lw=0.9, ls='--')
a1.set_ylabel('소비 c_t')
a1.set_ylim(0, 2.8)
a1.set_title('횡단조건이 경로를 고른다')
a1.legend(fontsize=8, loc='upper left')
a2.yaxis.set_major_locator(LogLocator(base=10.0))
a2.yaxis.set_major_formatter(pow10)
a2.yaxis.set_minor_formatter(NullFormatter())
a2.set_ylabel(r'$e^{-\rho t}\mu_t k_t$')
a2.set_xlabel('시간(년)')
a2.set_xlim(0, T_FIG)
plt.show()
# [예측 3] rho 를 0.03 -> 0.02 로 영구히 내린다 — 옛 정상상태에서의 점프.
RHO2 = 0.02
KS2, CS2 = steady_closed(RHO2)
_, DET2, LAM1b, LAM2b = eigen_hand(KS2, CS2, SIGMA, RHO2)
KA2, CA2 = arm_table(arm_branches(rho=RHO2))
c_lin = CS2 + LAM2b * (KS - KS2) # (eq-w11-16) 의 선형 근사
c_nl = float(np.interp(KS, KA2, CA2)) # 역행 사격으로 얻은 비선형 값
print('새 정상상태 k*=%.4f c*=%.4f lambda1=%.5f lambda2=%.5f 반감기=%.2f년'
% (KS2, CS2, LAM1b, LAM2b, np.log(2) / abs(LAM1b)))
print('옛 정상상태 (k0,c0)=(%.4f, %.4f) 은 새 phi 곡선 위에 그대로 있다 — phi(k0)=%.4f'
% (KS, CS, phi(KS)))
print('[예측 3] 선형 근사 c0=%.4f (%+.2f%%), 역행 사격 c0=%.4f (%+.2f%%)'
% (c_lin, 100 * (c_lin / CS - 1), c_nl, 100 * (c_nl / CS - 1)))
print('[예측 3] 두 값의 차이 %.2f%%p — 소비는 아래로 뛴다. 자본은 그 자리(스톡)'
% abs(100 * (c_lin / CS - 1) - 100 * (c_nl / CS - 1)))
print('점프 직후 저축률 1-c0/y0=%.4f (옛 s*=%.4f), 소비 성장률 (f\'(k0)-rho2-delta)/sigma=%.4f'
% (1 - c_nl / f(KS), S_STAR, (fp(KS) - RHO2 - DELTA) / SIGMA))새 정상상태 k*=10.3913 c*=1.6626 lambda1=-0.05191 lambda2=0.07191 반감기=13.35년
옛 정상상태 (k0,c0)=(8.5052, 1.6160) 은 새 phi 곡선 위에 그대로 있다 — phi(k0)=1.6160
[예측 3] 선형 근사 c0=1.5270 (-5.51%), 역행 사격 c0=1.5203 (-5.92%)
[예측 3] 두 값의 차이 0.42%p — 소비는 아래로 뛴다. 자본은 그 자리(스톡)
점프 직후 저축률 1-c0/y0=0.2552 (옛 s*=0.2083), 소비 성장률 (f'(k0)-rho2-delta)/sigma=0.0050
# 그림 3 — 선형 근사와 안장경로의 차이를 |k-k*| 에 대해 쌓는다.
kg = np.concatenate([np.linspace(0.35 * KS, 0.99 * KS, 60),
np.linspace(1.01 * KS, 2.5 * KS, 60)])
c_true = np.interp(kg, KA, CA)
c_lin_g = CS + LAM2 * (kg - KS)
gap = np.abs(c_lin_g - c_true) / c_true
d = np.abs(kg - KS)
ref = np.argmin(np.abs(d - 1.0))
near = d < 0.5 # 정상상태 근방의 국소 기울기
slope = np.polyfit(np.log(d[near]), np.log(gap[near]), 1)[0]
fig, ax = plt.subplots(figsize=(6.4, 4.2))
ax.loglog(d, gap, 'o', ms=3, color='C0', label='선형 근사의 상대오차')
ax.loglog(d, gap[ref] * (d / d[ref]) ** 2, color='0.45', lw=1.0, ls='--', label='기울기 2 기준선')
ax.axvline(abs(K0 - KS), color='C3', lw=1.0, ls=':')
for axis in (ax.xaxis, ax.yaxis):
axis.set_major_locator(LogLocator(base=10.0))
axis.set_major_formatter(pow10)
axis.set_minor_formatter(NullFormatter())
ax.set_xlabel(r'$|k-k^*|$')
ax.set_ylabel('상대오차')
ax.set_title('버린 항은 간극의 제곱으로 되살아난다')
ax.legend(fontsize=8)
plt.show()
print('k=%.0f (|k-k*|=%.2f): 선형 근사 %.4f 대 안장경로 %.4f — 오차 %.1f%%'
% (K0, abs(K0 - KS), C_LIN0, C_ARM0, 100 * (C_LIN0 / C_ARM0 - 1)))
print('근방 %d개 점(|k-k*|<0.5)의 로그-로그 기울기 적합 %.2f — 버린 항은 (k-k*)^2 이다'
% (near.sum(), slope))
k=3 (|k-k*|=5.51): 선형 근사 1.1329 대 안장경로 1.0142 — 오차 11.7%
근방 7개 점(|k-k*|<0.5)의 로그-로그 기울기 적합 2.02 — 버린 항은 (k-k*)^2 이다
# 그림 4 — sigma 스윕. 반감기 곡선이 솔로우 기준선을 어디서 지나는가.
SIG_CROSS = (RHO + DELTA) / (ALPHA * DELTA)
def halflife(sigma):
return np.log(2) / abs(eigen_hand(KS, CS, sigma)[2])
print('sigma 반감기(년) 솔로우 대비')
for s in (0.5, 1.0, 2.0, 3.0, SIG_CROSS, 6.0, 10.0):
hl = halflife(s)
tag = ('같다' if abs(hl - HALF_SOLOW) < 1e-9 * HALF_SOLOW # 허용오차를 먼저 본다
else ('빠르다' if hl < HALF_SOLOW else '느리다'))
print('%6.2f %10.2f %s' % (s, hl, tag))
print('교차 sigma=(rho+delta)/(alpha delta)=%.4f, 1/s*=%.4f, 솔로우 반감기=%.4f년'
% (SIG_CROSS, 1 / S_STAR, HALF_SOLOW))
print('교차점에서 반감기 차이 %.2e년 — 부동소수점 잔차다' % (halflife(SIG_CROSS) - HALF_SOLOW))
sg = np.logspace(np.log10(0.2), np.log10(50.0), 300)
hg = np.array([halflife(s) for s in sg])
fig, ax = plt.subplots(figsize=(6.4, 4.0))
ax.plot(sg, hg, color='C0', lw=2.0, label=r'램지 $\ln 2/|\lambda_1|$')
ax.axhline(HALF_SOLOW, color='0.45', lw=1.2, ls='--', label='솔로우 %.1f년' % HALF_SOLOW)
ax.plot([SIG_CROSS], [halflife(SIG_CROSS)], 'o', color='C3', ms=6, zorder=5)
ax.plot([SIGMA], [HALF], 'ko', ms=5, zorder=5)
ax.annotate(r'$\sigma=1/s^*=%.1f$' % SIG_CROSS, (SIG_CROSS, HALF_SOLOW),
textcoords='offset points', xytext=(8, -16), color='C3', fontsize=9)
ax.annotate(r'$\sigma=2$: %.1f년' % HALF, (SIGMA, HALF),
textcoords='offset points', xytext=(8, -16), fontsize=9)
ax.set_xscale('log')
ax.set_xlim(0.2, 50)
ax.set_ylim(0, 130)
ax.xaxis.set_major_locator(FixedLocator([0.2, 0.5, 1, 2, 5, 10, 20, 50]))
ax.xaxis.set_minor_locator(NullLocator())
ax.xaxis.set_major_formatter(FuncFormatter(lambda v, _: '%g' % v))
ax.set_xlabel(r'상대적 위험회피도 $\sigma$ (로그축)')
ax.set_ylabel('반감기(년)')
ax.set_title(r'반감기 곡선과 솔로우 기준선 — 교차점은 $\sigma=1/s^*$')
ax.legend(fontsize=8, loc='upper left')
plt.show()sigma 반감기(년) 솔로우 대비
0.50 5.41 빠르다
1.00 7.99 빠르다
2.00 12.00 빠르다
3.00 15.40 빠르다
4.80 20.79 같다
6.00 24.11 느리다
10.00 34.33 느리다
교차 sigma=(rho+delta)/(alpha delta)=4.8000, 1/s*=4.8000, 솔로우 반감기=20.7944년
교차점에서 반감기 차이 -1.07e-14년 — 부동소수점 잔차다
# 교차 sigma 에서 안장경로 위 저축률 — k 에 무관한가.
KAc, CAc = arm_table(arm_branches(sigma=SIG_CROSS))
print('sigma=%.1f 의 안장경로 위 저축률' % SIG_CROSS)
print(' k c_arm(k) 1-c_arm/f(k)')
for k in (1.0, 2.0, 3.0, 5.0, KS, 12.0, 20.0):
c = float(np.interp(k, KAc, CAc))
print('%6.2f %10.4f %13.4f' % (k, c, 1 - c / f(k)))
print('s* = alpha delta/(rho+delta) = %.4f — k 에 무관하다' % S_STAR)sigma=4.8 의 안장경로 위 저축률
k c_arm(k) 1-c_arm/f(k)
1.00 0.7917 0.2083
2.00 0.9974 0.2083
3.00 1.1418 0.2083
5.00 1.3537 0.2083
8.51 1.6160 0.2083
12.00 1.8125 0.2083
20.00 2.1489 0.2083
s* = alpha delta/(rho+delta) = 0.2083 — k 에 무관하다