# 1. 준모수 가법 분해 및 AR(1) 혁신항 추출 함수
def fit_semi_parametric(data_slice: pd.DataFrame, gamma: float = 0.25, t_half: float = 1460.0):
"""로그 가격을 거듭제곱 추세, 감쇠 조화 진동자, AR(1) 정상 잔차로 분해"""
t = data_slice["t"].values
p = data_slice["p_t"].values
ln_t = np.log(t)
omega = 2.0 * np.pi / t_half
sin_term = (t ** -gamma) * np.sin(omega * t)
cos_term = (t ** -gamma) * np.cos(omega * t)
X = np.column_stack([np.ones_like(t), ln_t, sin_term, cos_term])
beta, _, _, _ = np.linalg.lstsq(X, p, rcond=None)
p_trend = beta[0] + beta[1] * ln_t
c = beta[2] * sin_term + beta[3] * cos_term
residual = p - (p_trend + c)
# 잔차의 AR(1) 평균회귀 계수(rho) 및 혁신항(u_t) 추정
res_lag = residual[:-1]
res_curr = residual[1:]
rho = float(np.clip(np.sum(res_lag * res_curr) / np.sum(res_lag ** 2), 0.0, 0.999))
innovations = res_curr - rho * res_lag
innovations = innovations - np.mean(innovations)
return beta, p_trend, c, residual, rho, innovations
# 2. 준모수 경로 생성 함수 (외삽 추세 + 감쇠 사이클 + AR(1)-MBB 잔차 필터링)
def simulate_semi_parametric_paths(data_slice: pd.DataFrame, tf: int, gamma: float = 0.25,
t_half: float = 1460.0, block_len: int = 30,
m_paths: int = 10000, seed: int = 42):
"""구조적 성분 외삽과 AR(1) 평균회귀 잔차 MBB 합성 경로 산출"""
beta, _, _, residual, rho, innovations = fit_semi_parametric(data_slice, gamma, t_half)
t_last = data_slice.iloc[-1]["t"]
t_fut = np.arange(t_last + 1, t_last + tf + 1)
ln_tf = np.log(t_fut)
omega = 2.0 * np.pi / t_half
sin_f = (t_fut ** -gamma) * np.sin(omega * t_fut)
cos_f = (t_fut ** -gamma) * np.cos(omega * t_fut)
trend_future = beta[0] + beta[1] * ln_tf
c_future = beta[2] * sin_f + beta[3] * cos_f
# 혁신항(u_t)에 대한 이동 블록 부트스트랩
n_b = len(innovations) - block_len + 1
blocks = np.array([innovations[i : i + block_len] for i in range(n_b)])
k_blocks = int(np.ceil(tf / block_len))
np.random.seed(seed)
b_idx = np.random.randint(0, n_b, size=(m_paths, k_blocks))
sampled_u = blocks[b_idx].reshape(m_paths, -1)[:, :tf]
# AR(1) 평균회귀 필터링: eps_k = rho * eps_{k-1} + u_k^*
eps_paths = np.empty((m_paths, tf))
eps_prev = np.full(m_paths, residual[-1])
for k in range(tf):
eps_prev = rho * eps_prev + sampled_u[:, k]
eps_paths[:, k] = eps_prev
log_paths_future = trend_future + c_future + eps_paths
p0_log = data_slice.iloc[-1]["p_t"]
full_log_paths = np.hstack([np.full((m_paths, 1), p0_log), log_paths_future])
return np.exp(full_log_paths), beta, rho
# 3. 벤치마크 기준 시점(Cycle 2 Peak) 3단 구조 분해 시각화
target_slice = get_insample_data(df, ORIGINS["Cycle 2 Peak"], "All")
beta_est, p_tr, c_cyc, res_est, rho_est, _ = fit_semi_parametric(target_slice, gamma=0.25, t_half=1460.0)
adf_res_p = get_adf_pval(res_est, reg="c")
half_life_res = -np.log(2) / np.log(rho_est)
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 8.5), sharex=True)
ax1.plot(target_slice["date"], target_slice["p_t"], color="#2c3e50", lw=1.2, label="Actual Log Price ($p_t$)")
ax1.plot(target_slice["date"], p_tr, color="#e74c3c", lw=1.5, linestyle="--", label=rf"Power-Law Trend: $\hat{{p}}_{{\text{{trend}}}} = {beta_est[0]:.2f} + {beta_est[1]:.2f}\ln t$")
ax1.set_title("Structural Component 1: Power-Law Trend Extracted In-Sample (Cycle 2 Peak)", fontsize=10.5, fontweight="bold")
ax1.set_ylabel("Log Price", fontsize=9)
ax1.legend(loc="upper left", fontsize=8)
ax2.plot(target_slice["date"], c_cyc, color="#8e44ad", lw=1.2, label=rf"Halving Oscillator ($\gamma=0.25, T_{{\text{{half}}}}=1460$d)")
ax2.axhline(0, color="gray", lw=0.8, linestyle=":")
ax2.set_title("Structural Component 2: Damped Cyclical Harmonic Oscillator ($c(t)$)", fontsize=10.5, fontweight="bold")
ax2.set_ylabel("Cycle Component", fontsize=9)
ax2.legend(loc="upper left", fontsize=8)
ax3.plot(target_slice["date"], res_est, color="#16a085", lw=1.0, label=rf"Residual $\hat{{\epsilon}}_t$ (ADF $p={adf_res_p:.2e}$, $\hat{{\rho}}={rho_est:.4f}$, Half-life={half_life_res:.0f}d)")
ax3.axhline(0, color="black", lw=0.8, linestyle="--")
ax3.set_title("Structural Component 3: Stationary Residual ($\hat{\epsilon}_t$)", fontsize=10.5, fontweight="bold")
ax3.set_xlabel("Date", fontsize=9)
ax3.set_ylabel("Residual", fontsize=9)
ax3.legend(loc="upper left", fontsize=8)
plt.tight_layout()
plt.show()
# 4. 감쇠 계수(gamma) 및 반감 주기(T_half_life) 그리드 민감도 집계 및 2x2 팬 차트
semi_summary_rows = []
semi_cache = {}
for origin_name, origin_date in ORIGINS.items():
in_data = get_insample_data(df, origin_date, "All")
out_real = get_outsample_realized(df, origin_date, 1460)
tf_eval = len(out_real)
p_real_end = out_real.iloc[-1]["PriceUSD"]
for g in [0.0, 0.25, 0.50]:
paths_semi, b_vec, rho_val = simulate_semi_parametric_paths(in_data, tf_eval, gamma=g, t_half=1460.0, block_len=30)
semi_cache[(origin_name, g, 1460)] = paths_semi
q50 = np.median(paths_semi[:, -1])
bias_4y = compute_median_bias(paths_semi[:, -1], p_real_end)
cov_4y = compute_coverage_prob(paths_semi, out_real["PriceUSD"].values)
semi_summary_rows.append({
"Origin": origin_name,
"gamma": f"{g:.2f}",
"T_half": "1460d",
"Slope (b)": f"{b_vec[1]:.2f}",
"AR(1) rho": f"{rho_val:.4f}",
"Med 4Y*": f"${q50:,.0f}",
"Realized 4Y*": f"${p_real_end:,.0f}",
"Bias (4Y*)": f"{bias_4y:+.1%}",
"Cov (4Y*)": f"{cov_4y:.1%}"
})
# T_half_life in {1400, 1460, 1500} 민감도 추가 검증 (Cycle 3 Peak, gamma=0.25)
in_c3p = get_insample_data(df, ORIGINS["Cycle 3 Peak"], "All")
out_c3p = get_outsample_realized(df, ORIGINS["Cycle 3 Peak"], 1460)
for th in [1400, 1500]:
p_th, b_th, r_th = simulate_semi_parametric_paths(in_c3p, len(out_c3p), gamma=0.25, t_half=float(th), block_len=30)
semi_summary_rows.append({
"Origin": "Cycle 3 Peak (Period Test)",
"gamma": "0.25",
"T_half": f"{th}d",
"Slope (b)": f"{b_th[1]:.2f}",
"AR(1) rho": f"{r_th:.4f}",
"Med 4Y*": f"${np.median(p_th[:, -1]):,.0f}",
"Realized 4Y*": f"${out_c3p.iloc[-1]['PriceUSD']:,.0f}",
"Bias (4Y*)": f"{compute_median_bias(p_th[:, -1], out_c3p.iloc[-1]['PriceUSD']):+.1%}",
"Cov (4Y*)": f"{compute_coverage_prob(p_th, out_c3p['PriceUSD'].values):.1%}"
})
semi_table_df = pd.DataFrame(semi_summary_rows)
print(semi_table_df.to_string(index=False))
# 5. 4대 기준 시점별 준모수 시뮬레이션 2x2 팬 차트 (gamma in {0.0, 0.25, 0.50})
fig, axes = plt.subplots(2, 2, figsize=(14, 10), sharey=False)
gamma_Colors = {0.0: "#e67e22", 0.25: "#27ae60", 0.50: "#2980b9"}
for ax, (origin_name, origin_date) in zip(axes.flatten(), ORIGINS.items()):
realized_df = get_outsample_realized(df, origin_date, 1460)
t_days = np.arange(len(realized_df) + 1)
in_last = get_insample_data(df, origin_date, 1).iloc[-1]["PriceUSD"]
p_real_series = np.array([in_last] + realized_df["PriceUSD"].tolist())
ax.plot(t_days, p_real_series, color="black", lw=2.0, zorder=6, label="Realized Price")
for g in [0.0, 0.25, 0.50]:
paths = semi_cache[(origin_name, g, 1460)]
q05, q50, q95 = np.percentile(paths, [5, 50, 95], axis=0)
c = gamma_Colors[g]
ax.plot(t_days, q50, color=c, lw=1.4, label=rf"Median ($\gamma={g:.2f}$)")
if g == 0.25:
ax.fill_between(t_days, q05, q95, color=c, alpha=0.18, label=r"90% CI ($\gamma=0.25$)")
ax.set_yscale("log")
ax.set_title(f"Semi-Parametric MBB: {origin_name} ({origin_date.strftime('%Y-%m-%d')})", fontsize=11, fontweight="bold")
ax.set_xlabel("Days Elapsed from Cut-off ($k$)", fontsize=9)
ax.set_ylabel("Price in USD (Log Scale)", fontsize=9)
ax.legend(loc="lower right" if "Peak" in origin_name else "upper left", fontsize=7.5, frameon=True)
plt.tight_layout()
plt.show()
