wehub-resource-sync
642 lines
23 KiB
Python
642 lines
23 KiB
Python
# ---
|
||
# jupyter:
|
||
# jupytext:
|
||
# cell_metadata_filter: tags,-all
|
||
# text_representation:
|
||
# extension: .py
|
||
# format_name: percent
|
||
# format_version: '1.3'
|
||
# jupytext_version: 1.19.3
|
||
# kernelspec:
|
||
# display_name: Python 3 (ipykernel)
|
||
# language: python
|
||
# name: python3
|
||
# ---
|
||
|
||
# %% [markdown]
|
||
# # Uncertainty Features
|
||
#
|
||
# **Docker image**: `ml4t`
|
||
#
|
||
# This notebook demonstrates Bayesian and frequentist approaches to extracting
|
||
# **uncertainty features** — posterior distributions and prediction intervals
|
||
# become ML inputs, not just diagnostics.
|
||
#
|
||
# **Learning Objectives**:
|
||
# - Understand the uncertainty-as-feature principle
|
||
# - Implement a walk-forward stochastic volatility (SV) model with filtered posteriors
|
||
# - Extract `vol_posterior_std`, `vol_ci_width`, `vol_of_vol`, and `vol_persistence`
|
||
# - Compute ARIMA forecast uncertainty on log-volatility with proper back-transform
|
||
#
|
||
# **Book Reference**: Chapter 9, Section 9.4 (Uncertainty Features)
|
||
#
|
||
# **Prerequisites**: `08_garch_volatility` for frequentist volatility modeling
|
||
# (GARCH provides point estimates; this notebook adds uncertainty quantification).
|
||
# `09_har_rough_volatility` for realized volatility estimators.
|
||
|
||
# %%
|
||
"""Uncertainty Features — walk-forward SV posteriors and ARIMA forecast intervals."""
|
||
|
||
from datetime import datetime
|
||
|
||
import arviz as az
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import pandas as pd
|
||
import polars as pl
|
||
import pymc as pm
|
||
from IPython.display import display
|
||
from ml4t.diagnostic.evaluation.autocorrelation import analyze_autocorrelation
|
||
from ml4t.diagnostic.evaluation.stationarity import analyze_stationarity
|
||
from ml4t.engineer.features.volatility import garman_klass_volatility
|
||
from statsforecast import StatsForecast
|
||
from statsforecast.models import AutoARIMA
|
||
|
||
from data import load_etfs
|
||
from utils.paths import get_case_study_dir
|
||
from utils.reproducibility import set_global_seeds
|
||
|
||
# %% tags=["parameters"]
|
||
# Production defaults — Papermill injects overrides for CI
|
||
N_DRAWS = 1000
|
||
N_TUNE = 2000
|
||
N_CHAINS = 2
|
||
REFIT_INTERVAL = 63 # Quarterly (~252/4 trading days)
|
||
ARIMA_WINDOW = 252
|
||
TRAIN_DAYS = 252
|
||
SEED = 42
|
||
|
||
# %%
|
||
set_global_seeds(SEED)
|
||
|
||
# %% [markdown]
|
||
# ## Load Data
|
||
|
||
# %%
|
||
etfs = load_etfs(symbols=["SPY"])
|
||
sp500 = etfs.select(["timestamp", "open", "high", "low", "close"]).sort("timestamp")
|
||
|
||
start_date = datetime(2015, 1, 1)
|
||
end_date = datetime(2024, 12, 31)
|
||
|
||
sp500 = sp500.filter((pl.col("timestamp") >= start_date) & (pl.col("timestamp") <= end_date))
|
||
|
||
sp500 = sp500.with_columns(
|
||
returns=pl.col("close").pct_change() * 100,
|
||
).drop_nulls()
|
||
|
||
returns = sp500.to_pandas().set_index("timestamp")["returns"]
|
||
|
||
print(f"Returns: {len(returns)} observations")
|
||
print(f"Date range: {returns.index.min()} to {returns.index.max()}")
|
||
print(f"Sample mean: {returns.mean():.4f}%, Sample std: {returns.std():.4f}%")
|
||
|
||
# %% [markdown]
|
||
# ## Garman-Klass Realized Volatility
|
||
#
|
||
# ARIMA forecast uncertainty requires a volatility target series. We use the
|
||
# Garman-Klass estimator (Garman and Klass, 1980) rather than close-to-close
|
||
# rolling standard deviation. GK incorporates open, high, low, and close prices,
|
||
# producing a more efficient estimator — see `09_har_rough_volatility` for the
|
||
# full comparison of range-based estimators.
|
||
|
||
# %%
|
||
sp500_gk = sp500.with_columns(
|
||
rv_gk_21=garman_klass_volatility("open", "high", "low", "close", period=21),
|
||
).drop_nulls(subset=["rv_gk_21"])
|
||
|
||
rv_gk = sp500_gk.select(["timestamp", "rv_gk_21"]).to_pandas().set_index("timestamp")["rv_gk_21"]
|
||
|
||
print(f"Garman-Klass RV (21-day): {len(rv_gk)} observations")
|
||
print(f"Mean: {rv_gk.mean():.4f}, Std: {rv_gk.std():.4f}")
|
||
|
||
# %%
|
||
fig, ax = plt.subplots(figsize=(14, 4))
|
||
ax.plot(rv_gk.index, rv_gk.values, linewidth=0.5, alpha=0.8)
|
||
ax.set_title("21-Day Garman-Klass Realized Volatility (SPY)")
|
||
ax.set_ylabel("Annualized Volatility")
|
||
ax.set_xlabel("Date")
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
# %% [markdown]
|
||
# ## The Uncertainty-as-Feature Principle
|
||
#
|
||
# Any point estimate becomes a distributional feature under Bayesian inference.
|
||
# A brief example: the Sharpe ratio. The frequentist estimate is a single number;
|
||
# the Bayesian posterior reveals how uncertain that number is. The same principle
|
||
# applies to volatility, hedge ratios, and any model parameter — the posterior
|
||
# width is itself a feature. Chapter 17 develops Bayesian Sharpe estimation fully;
|
||
# here we apply the principle to volatility and forecasts.
|
||
|
||
# %%
|
||
freq_sharpe = (returns.mean() / returns.std()) * np.sqrt(252)
|
||
print(f"Frequentist Sharpe (annualized): {freq_sharpe:.3f}")
|
||
print("A single number with no uncertainty quantification.")
|
||
print("The Bayesian version (Chapter 17) produces a full posterior distribution,")
|
||
print("yielding credible intervals and probability statements like P(SR > 0).")
|
||
|
||
# %% [markdown]
|
||
# ## Stochastic Volatility Model — Walk-Forward
|
||
#
|
||
# The canonical SV model specifies an AR(1) process in log-volatility:
|
||
#
|
||
# $$\log(\sigma_t) = \mu_h + \phi \cdot (\log(\sigma_{t-1}) - \mu_h) + \sigma_\eta \cdot \eta_t$$
|
||
#
|
||
# We implement this using `pm.AR` with the centered parameterization:
|
||
# $\tilde{h}_t = \phi \cdot \tilde{h}_{t-1} + \sigma_\eta \cdot \eta_t$
|
||
# where $\tilde{h}_t = h_t - \mu_h$. The observation model uses Student-t
|
||
# rather than Gaussian to handle fat tails in financial returns.
|
||
#
|
||
# ### Walk-Forward Protocol
|
||
#
|
||
# To maintain point-in-time integrity, we refit the SV model at quarterly
|
||
# boundaries on a trailing training window. At each refit, we extract the
|
||
# **filtered final state** — the posterior at the last time step, conditioned
|
||
# only on data up to that point. Between refits, features carry forward the
|
||
# most recent posterior state. This mirrors the GARCH walk-forward in
|
||
# `08_garch_volatility` but with quarterly cadence due to MCMC cost.
|
||
|
||
# %% [markdown]
|
||
# ### SV Model Specification
|
||
#
|
||
# Define the model as a function for reuse across walk-forward folds.
|
||
|
||
|
||
# %%
|
||
def fit_sv_model(returns_array, n_draws, n_tune, n_chains):
|
||
"""Fit AR(1) stochastic volatility with Student-t observations.
|
||
|
||
Returns the MCMC trace (InferenceData object).
|
||
"""
|
||
with pm.Model() as sv_model:
|
||
mu_h = pm.Normal("mu_h", mu=0, sigma=5)
|
||
phi = pm.Uniform("phi", lower=0, upper=1)
|
||
sigma_eta = pm.Exponential("sigma_eta", lam=2)
|
||
nu = pm.Deterministic("nu", pm.Gamma("nu_minus2", alpha=2, beta=0.1) + 2)
|
||
|
||
h_centered = pm.AR(
|
||
"h_centered",
|
||
rho=[phi],
|
||
sigma=sigma_eta,
|
||
init_dist=pm.Normal.dist(0, 1),
|
||
shape=len(returns_array),
|
||
)
|
||
h = pm.Deterministic("h", h_centered + mu_h)
|
||
volatility = pm.Deterministic("volatility", pm.math.exp(h / 2))
|
||
|
||
pm.StudentT("obs", nu=nu, mu=0, sigma=volatility, observed=returns_array)
|
||
|
||
with sv_model:
|
||
trace = pm.sample(
|
||
n_draws,
|
||
tune=n_tune,
|
||
chains=n_chains,
|
||
cores=1,
|
||
progressbar=False,
|
||
random_seed=SEED,
|
||
target_accept=0.99,
|
||
)
|
||
return trace
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Walk-Forward Execution
|
||
#
|
||
# Refit every `REFIT_INTERVAL` trading days (~quarterly). Each refit trains
|
||
# on the preceding `TRAIN_DAYS` observations. Between refits, we carry forward
|
||
# the filtered posterior from the training window's final state.
|
||
|
||
# %%
|
||
forecast_start = len(returns) - TRAIN_DAYS
|
||
refit_points = list(range(forecast_start, len(returns), REFIT_INTERVAL))
|
||
n_refits = len(refit_points)
|
||
|
||
print(f"Walk-forward: {n_refits} refits, interval={REFIT_INTERVAL} days")
|
||
print(f"Training window: {TRAIN_DAYS} days, forecast horizon: {len(returns) - forecast_start} days")
|
||
|
||
sv_features_list = []
|
||
last_trace = None
|
||
|
||
# %%
|
||
for refit_idx, t in enumerate(refit_points):
|
||
train_start = max(0, t - TRAIN_DAYS)
|
||
train_data = returns.iloc[train_start:t].values
|
||
|
||
carry_end = min(t + REFIT_INTERVAL, len(returns))
|
||
carry_dates = returns.index[t:carry_end]
|
||
|
||
print(
|
||
f" Refit {refit_idx + 1}/{n_refits}: train [{train_start}:{t}] ({len(train_data)} obs), "
|
||
f"carry [{t}:{carry_end}] ({len(carry_dates)} days)"
|
||
)
|
||
|
||
trace = fit_sv_model(train_data, N_DRAWS, N_TUNE, N_CHAINS)
|
||
last_trace = trace
|
||
|
||
# Extract FILTERED final state — posterior at the last training time step
|
||
vol_final = trace.posterior["volatility"].values[:, :, -1].flatten()
|
||
phi_post = trace.posterior["phi"].values.flatten()
|
||
sigma_eta_post = trace.posterior["sigma_eta"].values.flatten()
|
||
nu_post = trace.posterior["nu"].values.flatten()
|
||
|
||
features = {
|
||
"vol_posterior_mean": vol_final.mean(),
|
||
"vol_posterior_std": vol_final.std(),
|
||
"vol_ci_width": np.percentile(vol_final, 97.5) - np.percentile(vol_final, 2.5),
|
||
"vol_of_vol": sigma_eta_post.mean(),
|
||
"vol_persistence": phi_post.mean(),
|
||
}
|
||
|
||
# Carry forward to each day until next refit
|
||
for dt in carry_dates:
|
||
row = {"timestamp": dt, **features}
|
||
sv_features_list.append(row)
|
||
|
||
# Per-refit diagnostics
|
||
n_div = int(trace.sample_stats["diverging"].values.sum())
|
||
print(
|
||
f" phi={phi_post.mean():.3f}, sigma_eta={sigma_eta_post.mean():.3f}, "
|
||
f"nu={nu_post.mean():.1f}, divergences={n_div}"
|
||
)
|
||
|
||
sv_features_df = pd.DataFrame(sv_features_list).set_index("timestamp")
|
||
print(f"\nSV features: {len(sv_features_df)} daily observations ({n_refits} refit points)")
|
||
|
||
# %% [markdown]
|
||
# ### MCMC Diagnostics
|
||
#
|
||
# Beyond divergence counts, we check $\hat{R}$ (convergence across chains) and
|
||
# effective sample size (ESS). $\hat{R} > 1.05$ or ESS below 100 indicate the
|
||
# sampler has not explored the posterior adequately.
|
||
#
|
||
# SV models are notoriously difficult for NUTS because the latent AR(1) path
|
||
# creates a funnel geometry between $\phi$, $\sigma_\eta$, and the 252 latent
|
||
# states. Expect lower ESS for these parameters than for $\mu_h$ and $\nu$,
|
||
# which are better identified. This is a genuine operational consideration:
|
||
# production SV implementations often use particle MCMC or sequential Monte
|
||
# Carlo rather than NUTS to handle this geometry.
|
||
|
||
# %%
|
||
# Diagnostics from the last refit
|
||
diag_vars = ["phi", "sigma_eta", "mu_h", "nu"]
|
||
summary = az.summary(last_trace, var_names=diag_vars)
|
||
display(summary)
|
||
|
||
n_divergences = int(last_trace.sample_stats["diverging"].values.sum())
|
||
print(f"\nDivergences: {n_divergences}")
|
||
|
||
# %% [markdown]
|
||
# ### Walk-Forward SV Features
|
||
#
|
||
# The plot shows piecewise-constant features that update at each quarterly
|
||
# refit. This honestly reflects the production protocol: features carry the
|
||
# last filtered state forward between re-estimations.
|
||
|
||
# %%
|
||
fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)
|
||
|
||
ax = axes[0]
|
||
wf_dates = sv_features_df.index
|
||
ax.plot(
|
||
returns.loc[wf_dates].index,
|
||
returns.loc[wf_dates].values,
|
||
linewidth=0.4,
|
||
alpha=0.6,
|
||
label="Returns",
|
||
)
|
||
ax.fill_between(
|
||
wf_dates,
|
||
-2 * sv_features_df["vol_posterior_mean"],
|
||
2 * sv_features_df["vol_posterior_mean"],
|
||
alpha=0.3,
|
||
color="red",
|
||
label="±2σ (filtered)",
|
||
)
|
||
ax.set_title("Returns with Walk-Forward Volatility Bands")
|
||
ax.set_ylabel("Return (%)")
|
||
ax.legend()
|
||
|
||
ax = axes[1]
|
||
ax.plot(
|
||
wf_dates,
|
||
sv_features_df["vol_posterior_mean"],
|
||
linewidth=1.5,
|
||
label="Posterior mean (filtered)",
|
||
drawstyle="steps-post",
|
||
)
|
||
ax.fill_between(
|
||
wf_dates,
|
||
sv_features_df["vol_posterior_mean"] - sv_features_df["vol_ci_width"] / 2,
|
||
sv_features_df["vol_posterior_mean"] + sv_features_df["vol_ci_width"] / 2,
|
||
alpha=0.3,
|
||
label="95% CI width",
|
||
)
|
||
for rp in refit_points:
|
||
if rp < len(returns):
|
||
ax.axvline(returns.index[rp], color="gray", linestyle=":", linewidth=0.5)
|
||
ax.set_title("Filtered Posterior Volatility (Step-Function at Refit Points)")
|
||
ax.set_ylabel("Volatility (%)")
|
||
ax.set_xlabel("Date")
|
||
ax.legend()
|
||
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
# %% [markdown]
|
||
# **Interpretation**: The step-function updates show how SV features evolve at
|
||
# quarterly refits. Between refits, features are constant — reflecting the
|
||
# production reality that MCMC is too expensive to run daily. The CI width
|
||
# varies across refits, widening during uncertain periods.
|
||
|
||
# %% [markdown]
|
||
# ### Parameter Posteriors (Last Refit)
|
||
|
||
# %%
|
||
phi_samples = last_trace.posterior["phi"].values.flatten()
|
||
sigma_eta_samples = last_trace.posterior["sigma_eta"].values.flatten()
|
||
mu_h_samples = last_trace.posterior["mu_h"].values.flatten()
|
||
nu_samples = last_trace.posterior["nu"].values.flatten()
|
||
|
||
print("=== SV Parameter Posteriors (Last Refit) ===")
|
||
print(f"phi (persistence): mean={phi_samples.mean():.3f}, std={phi_samples.std():.3f}")
|
||
print(
|
||
f"sigma_eta (vol-of-vol): mean={sigma_eta_samples.mean():.3f}, std={sigma_eta_samples.std():.3f}"
|
||
)
|
||
print(f"mu_h (mean log-vol): mean={mu_h_samples.mean():.3f}, std={mu_h_samples.std():.3f}")
|
||
print(f"nu (Student-t df): mean={nu_samples.mean():.1f}, std={nu_samples.std():.1f}")
|
||
|
||
fig, axes = plt.subplots(1, 4, figsize=(16, 4))
|
||
|
||
for ax, samples, label, title in zip(
|
||
axes,
|
||
[phi_samples, sigma_eta_samples, mu_h_samples, nu_samples],
|
||
["φ (persistence)", "σ_η (vol-of-vol)", "μ_h (mean log-vol)", "ν (Student-t df)"],
|
||
["Volatility Persistence", "Volatility of Volatility", "Mean Log-Volatility", "Tail Thickness"],
|
||
strict=False,
|
||
):
|
||
ax.hist(samples, bins=30, density=True, alpha=0.7, edgecolor="white")
|
||
ax.axvline(samples.mean(), color="red", linestyle="--", label=f"Mean: {samples.mean():.2f}")
|
||
ax.set_xlabel(label)
|
||
ax.set_title(title)
|
||
ax.legend(fontsize=8)
|
||
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
# %% [markdown]
|
||
# **Interpretation**: The $\phi$ posterior captures volatility persistence — values
|
||
# near 1 mean shocks decay slowly (clustering), lower values indicate faster
|
||
# mean-reversion. The $\nu$ posterior indicates tail thickness: values near 4–6
|
||
# confirm that Student-t is substantially better than Gaussian ($\nu \to \infty$);
|
||
# values above 30 suggest Gaussian is adequate. Both $\sigma_\eta$ (vol-of-vol)
|
||
# and $\phi$ are re-estimated quarterly and serve as slowly varying features.
|
||
|
||
# %% [markdown]
|
||
# ### Scalability
|
||
#
|
||
# Full MCMC-based SV is practical for small universes (tens of assets). For
|
||
# large panels (1000+ assets), pragmatic alternatives include:
|
||
# - **Sector aggregates**: Run SV on sector ETFs, broadcast uncertainty features
|
||
# - **Variational inference**: `pm.fit()` gives faster but less accurate posteriors
|
||
# - **GARCH baseline with SV overlay**: Reserve SV for portfolio-level risk
|
||
#
|
||
# See the chapter text (§9.4) for universe-size guidance.
|
||
|
||
# %% [markdown]
|
||
# ## ARIMA Forecast Uncertainty on Log-Volatility
|
||
#
|
||
# ARIMA produces prediction intervals alongside point forecasts. We apply it
|
||
# to **log-volatility** rather than raw volatility for two reasons:
|
||
#
|
||
# 1. Log-transform ensures forecasts stay positive (volatility cannot be negative)
|
||
# 2. ARIMA residuals on log-vol are more Gaussian, improving interval coverage
|
||
#
|
||
# The back-transform uses the log-normal property: $\exp(\mu_f)$ gives the
|
||
# **median** forecast on the original scale (not the mean, which is
|
||
# $\exp(\mu_f + \sigma_f^2/2)$). The prediction interval is
|
||
# $[\exp(\mu_f - z \cdot \sigma_f),\; \exp(\mu_f + z \cdot \sigma_f)]$.
|
||
# Both the median and the interval are exact, not approximations.
|
||
|
||
# %% [markdown]
|
||
# ### Pre-ARIMA Diagnostics
|
||
#
|
||
# Before fitting, verify that log-volatility is stationary (confirming $d=0$)
|
||
# and examine autocorrelation structure.
|
||
|
||
# %%
|
||
log_rv = np.log(rv_gk.clip(lower=1e-8))
|
||
|
||
stat_check = analyze_stationarity(log_rv.dropna().values)
|
||
acf_analysis = analyze_autocorrelation(log_rv.dropna().values)
|
||
|
||
print("=== Pre-ARIMA Diagnostics (log Garman-Klass RV) ===")
|
||
print(f"Stationarity: {stat_check.consensus} (agreement: {stat_check.agreement_score:.2f})")
|
||
print(f"Suggested ARIMA order: {acf_analysis.suggested_arima_order}")
|
||
|
||
# %% [markdown]
|
||
# ### Order Selection and Rolling Forecast via Nixtla
|
||
#
|
||
# We use Nixtla's `AutoARIMA` (compiled C, AIC-based order selection) on the
|
||
# first training window to choose (p, d, q), then roll forward with that fixed
|
||
# order using Nixtla's `ARIMA` for speed. The model operates on
|
||
# $\log(\text{RV}_{GK})$; we back-transform to the original scale afterward.
|
||
#
|
||
# $\exp(\hat{\mu}_f)$ gives the **median** forecast on the original scale
|
||
# (not the mean $\exp(\hat{\mu}_f + \hat{\sigma}_f^2/2)$), which is more
|
||
# robust for position sizing.
|
||
|
||
# %%
|
||
from statsforecast.models import ARIMA as NixtlaARIMA
|
||
|
||
log_rv_clean = log_rv.dropna()
|
||
|
||
# Step 1: AutoARIMA order selection on first training window
|
||
first_window_df = pd.DataFrame(
|
||
{
|
||
"unique_id": "SPY",
|
||
"ds": log_rv_clean.index[:ARIMA_WINDOW],
|
||
"y": log_rv_clean.values[:ARIMA_WINDOW],
|
||
}
|
||
)
|
||
|
||
sf_auto = StatsForecast(models=[AutoARIMA(season_length=1)], freq="B", n_jobs=1)
|
||
sf_auto.fit(first_window_df)
|
||
selected_order = sf_auto.fitted_[0, 0].model_["arma"]
|
||
p, q, P, Q, s, d, D = selected_order
|
||
print(f"AutoARIMA selected: ARIMA({p},{d},{q})")
|
||
|
||
# %%
|
||
# Step 2: Rolling forecast with fixed order (fast)
|
||
n_total = len(log_rv_clean)
|
||
n_windows = n_total - ARIMA_WINDOW
|
||
|
||
arima_input = pd.DataFrame(
|
||
{
|
||
"unique_id": "SPY",
|
||
"ds": log_rv_clean.index,
|
||
"y": log_rv_clean.values,
|
||
}
|
||
)
|
||
|
||
sf = StatsForecast(
|
||
models=[NixtlaARIMA(order=(p, d, q), season_length=1)],
|
||
freq="B",
|
||
n_jobs=1,
|
||
)
|
||
cv_result = sf.cross_validation(
|
||
df=arima_input,
|
||
h=1,
|
||
step_size=1,
|
||
n_windows=n_windows,
|
||
level=[95],
|
||
)
|
||
|
||
# Back-transform from log scale to original scale
|
||
col = "ARIMA"
|
||
arima_df = pd.DataFrame(
|
||
{
|
||
"timestamp": cv_result["ds"].values,
|
||
"arima_forecast": np.exp(cv_result[col].values),
|
||
"arima_lo_95": np.exp(cv_result[f"{col}-lo-95"].values),
|
||
"arima_hi_95": np.exp(cv_result[f"{col}-hi-95"].values),
|
||
"actual": np.exp(cv_result["y"].values),
|
||
}
|
||
)
|
||
arima_df["arima_ci_width"] = arima_df["arima_hi_95"] - arima_df["arima_lo_95"]
|
||
arima_df["forecast_uncertainty_ratio"] = arima_df["arima_ci_width"] / arima_df["arima_forecast"]
|
||
|
||
# Forecast std from the CI width on original scale (robust to negative log-CIs)
|
||
arima_df["arima_forecast_std"] = arima_df["arima_ci_width"] / (2 * 1.96)
|
||
|
||
arima_df = arima_df.set_index("timestamp")
|
||
|
||
print(f"ARIMA({p},{d},{q}) forecast features: {len(arima_df)} observations")
|
||
|
||
# %%
|
||
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)
|
||
|
||
ax = axes[0]
|
||
ax.plot(
|
||
arima_df.index,
|
||
arima_df["arima_forecast"],
|
||
linewidth=0.8,
|
||
label=f"ARIMA({p},{d},{q}) forecast (median)",
|
||
)
|
||
ax.plot(arima_df.index, arima_df["actual"], linewidth=0.5, alpha=0.5, label="Actual GK RV")
|
||
ax.fill_between(
|
||
arima_df.index,
|
||
arima_df["arima_forecast"] - 1.96 * arima_df["arima_forecast_std"],
|
||
arima_df["arima_forecast"] + 1.96 * arima_df["arima_forecast_std"],
|
||
alpha=0.2,
|
||
label="95% PI",
|
||
)
|
||
ax.set_title(f"ARIMA({p},{d},{q}) Forecast on Garman-Klass RV (Log-Normal Back-Transform)")
|
||
ax.set_ylabel("Annualized Volatility")
|
||
ax.legend()
|
||
|
||
ax = axes[1]
|
||
ax.plot(arima_df.index, arima_df["arima_forecast_std"], linewidth=0.8, color="orange")
|
||
ax.set_title("Forecast Standard Error on Log Scale (arima_forecast_std)")
|
||
ax.set_ylabel("Log-Scale Std Error")
|
||
|
||
ax = axes[2]
|
||
ax.plot(arima_df.index, arima_df["forecast_uncertainty_ratio"], linewidth=0.8, color="purple")
|
||
ax.axhline(
|
||
arima_df["forecast_uncertainty_ratio"].median(),
|
||
color="red",
|
||
linestyle="--",
|
||
linewidth=0.5,
|
||
label=f"Median: {arima_df['forecast_uncertainty_ratio'].median():.2f}",
|
||
)
|
||
ax.set_title("Forecast Uncertainty Ratio (CI Width / Forecast Level)")
|
||
ax.set_ylabel("Ratio")
|
||
ax.set_xlabel("Date")
|
||
ax.legend()
|
||
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
# %% [markdown]
|
||
# **Interpretation**: The forecast standard error spikes during volatile periods
|
||
# (COVID, 2022 rate hikes) — precisely when the model is least confident.
|
||
# The log-normal back-transform produces asymmetric prediction intervals on the
|
||
# original scale: wider on the upside, reflecting that volatility spikes are
|
||
# unbounded while the floor is zero. The order was selected by AutoARIMA on
|
||
# the first window; periodic reselection would guard against structural change.
|
||
|
||
# %% [markdown]
|
||
# ## Feature Summary
|
||
|
||
# %%
|
||
print("SV Uncertainty Features (Walk-Forward):")
|
||
display(sv_features_df.describe().round(4))
|
||
|
||
print("ARIMA Uncertainty Features (AutoARIMA on Log-Vol):")
|
||
display(
|
||
arima_df[["arima_forecast_std", "arima_ci_width", "forecast_uncertainty_ratio"]]
|
||
.describe()
|
||
.round(4)
|
||
)
|
||
|
||
# %% [markdown]
|
||
# ### Production Notes
|
||
#
|
||
# The walk-forward protocol implemented above is the correct production pattern:
|
||
# SV features are re-estimated at periodic boundaries and carried forward between
|
||
# refits. ARIMA features update daily with negligible cost (refit on a rolling
|
||
# window and call `get_forecast()`).
|
||
#
|
||
# For production hardening beyond this notebook:
|
||
# - **Particle filtering** provides online posterior updates without full MCMC
|
||
# re-estimation — useful for daily SV updates between quarterly refits
|
||
# - AutoARIMA reselects the order at each window automatically; for even
|
||
# faster rolling forecasts, fix the order after an initial selection pass
|
||
|
||
# %% [markdown]
|
||
# ## Save Uncertainty Features
|
||
|
||
# %%
|
||
MODEL_DIR = get_case_study_dir("etfs") / "models" / "time_series"
|
||
MODEL_DIR.mkdir(parents=True, exist_ok=True)
|
||
|
||
arima_save_cols = [
|
||
"arima_forecast",
|
||
"arima_forecast_std",
|
||
"arima_ci_width",
|
||
"forecast_uncertainty_ratio",
|
||
]
|
||
arima_output = pl.from_pandas(arima_df[arima_save_cols].reset_index())
|
||
arima_path = MODEL_DIR / "arima_uncertainty.parquet"
|
||
arima_output.write_parquet(arima_path)
|
||
print(f"Saved ARIMA uncertainty features to {arima_path}")
|
||
print(f" Shape: {arima_output.shape}")
|
||
|
||
sv_output = pl.from_pandas(sv_features_df.reset_index())
|
||
sv_path = MODEL_DIR / "sv_uncertainty.parquet"
|
||
sv_output.write_parquet(sv_path)
|
||
print(f"Saved SV uncertainty features to {sv_path}")
|
||
print(f" Shape: {sv_output.shape}")
|
||
|
||
# %% [markdown]
|
||
# ## Key Takeaways
|
||
#
|
||
# 1. **Walk-forward SV preserves point-in-time integrity** — the filtered final
|
||
# state from each training window produces features conditioned only on past
|
||
# data, unlike full-sample smoothed posteriors
|
||
# 2. **Student-t observation model handles fat tails** — the estimated $\nu$
|
||
# posterior confirms whether Gaussian is adequate; values near 4–6 show
|
||
# Student-t materially improves the fit
|
||
# 3. **Vol persistence ($\phi$) and vol-of-vol ($\sigma_\eta$)** are complementary:
|
||
# $\phi$ near 1 means shocks decay slowly; high $\sigma_\eta$ means the
|
||
# volatility process itself is unstable
|
||
# 4. **Log-vol ARIMA ensures positive forecasts** — the log-normal back-transform
|
||
# produces exact, asymmetric prediction intervals on the original scale
|
||
# 5. **MCMC diagnostics go beyond divergences** — $\hat{R}$, ESS, and the $\nu$
|
||
# posterior all provide information about posterior reliability
|
||
# 6. **Computational cost scales with universe** — full MCMC is practical for
|
||
# tens of assets; large panels require variational inference or sector aggregates
|
||
#
|
||
# **Previous**: `09_har_rough_volatility` for multi-horizon volatility decomposition.
|
||
# **Next**: `11_hmm_regimes` for regime detection via Hidden Markov Models.
|