{ "cells": [ { "cell_type": "markdown", "id": "56613d1a", "metadata": { "papermill": { "duration": 0.004031, "end_time": "2026-03-24T02:03:45.938079+00:00", "exception": false, "start_time": "2026-03-24T02:03:45.934048+00:00", "status": "completed" } }, "source": [ "# ETFs: Model-Based Features (Per-Fold Temporal Models)\n", "\n", "This notebook fits temporal models and extracts features that capture\n", "latent market dynamics. All models are fit **per CV fold** on training\n", "data only, eliminating parameter-level look-ahead bias. It produces\n", "features for three model families:\n", "\n", "1. **HMM Regime Detection**: 2-state Gaussian HMM on aggregate market (SPY)\n", " with filtered (causal) probabilities, regime transition indicators, and\n", " regime duration features.\n", "2. **Fractional Differencing**: Memory-preserving stationarity transforms\n", " on 10 reference ETFs spanning all major asset classes.\n", "3. **GARCH(1,1) Conditional Volatility**: Per-ETF volatility forecasts that\n", " provide each asset's own risk dynamics.\n", "\n", "Each row in the output carries a `fold` column identifying which fold's\n", "model produced it. This enables downstream CV to use the correct\n", "(non-leaked) features for each fold.\n", "\n", "## Learning Objectives\n", "\n", "- Fit temporal models per CV fold to avoid parameter-level look-ahead\n", "- Fit a 2-state HMM with k-means initialization and multiple restarts\n", "- Compute filtered (not smoothed) probabilities for production use\n", "- Derive regime transition and duration features from filtered probs\n", "- Apply fractional differencing with fixed $d$ values by asset class\n", "- Fit per-asset GARCH(1,1) with frozen-parameter filtering\n", "- Combine date-level and per-asset features into a per-fold panel\n", "\n", "## Book Reference\n", "\n", "Chapter 9, Sections 9.1 (Fractional Differencing), 9.3 (GARCH), and\n", "9.5 (Regime Features)\n", "\n", "## Prerequisites\n", "\n", "- [`02_labels`](02_labels.ipynb) (produces label parquet files)\n", "- `03_financial_features.py` (produces `features/financial.parquet`)" ] }, { "cell_type": "code", "execution_count": 2, "id": "5617b4d5", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:45.945084Z", "iopub.status.busy": "2026-03-24T02:03:45.944992Z", "iopub.status.idle": "2026-03-24T02:03:47.559470Z", "shell.execute_reply": "2026-03-24T02:03:47.559098Z" }, "lines_to_next_cell": 2, "papermill": { "duration": 1.618165, "end_time": "2026-03-24T02:03:47.560107+00:00", "exception": false, "start_time": "2026-03-24T02:03:45.941942+00:00", "status": "completed" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "~/ml4t/libraries/ml4t-engineer/src/ml4t/engineer/features/ml/__init__.py:9: UserWarning: Feature 'cyclical_encode': lookback=0 but has period/window parameter. Consider using lookback='period' or specifying the actual lookback.\n", " from ml4t.engineer.features.ml.cyclical_encode import * # noqa: F403\n" ] } ], "source": [ "\"\"\"ETFs: Model-Based Features (per-fold HMM + FFD + GARCH).\"\"\"\n", "\n", "import warnings\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import polars as pl\n", "import yaml\n", "from arch import arch_model\n", "from hmmlearn.hmm import GaussianHMM\n", "from ml4t.diagnostic.evaluation.stats import robust_ic\n", "from ml4t.engineer.features.fdiff import ffdiff\n", "from sklearn.cluster import KMeans\n", "\n", "from data import load_etfs\n", "from utils.cv_splits import generate_cv_splits, load_evaluation_config\n", "from utils.paths import get_case_study_dir\n", "\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": 3, "id": "5c611ee9", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:47.564876Z", "iopub.status.busy": "2026-03-24T02:03:47.564598Z", "iopub.status.idle": "2026-03-24T02:03:47.566450Z", "shell.execute_reply": "2026-03-24T02:03:47.566172Z" }, "papermill": { "duration": 0.004592, "end_time": "2026-03-24T02:03:47.566726+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.562134+00:00", "status": "completed" }, "tags": [ "parameters" ] }, "outputs": [], "source": [ "# Production defaults — Papermill injects overrides for CI\n", "START_DATE = None # None = use full dataset\n", "N_RESTARTS = 10\n", "GARCH_MIN_OBS = 504 # Minimum observations for GARCH fit (~2 years)\n", "MAX_SYMBOLS = 0 # 0 = all symbols (production)" ] }, { "cell_type": "code", "execution_count": 4, "id": "37ed4132", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:47.571043Z", "iopub.status.busy": "2026-03-24T02:03:47.570957Z", "iopub.status.idle": "2026-03-24T02:03:47.595274Z", "shell.execute_reply": "2026-03-24T02:03:47.594719Z" }, "papermill": { "duration": 0.027018, "end_time": "2026-03-24T02:03:47.595605+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.568587+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Prices: 470,662 rows, 100 assets\n" ] } ], "source": [ "CASE_DIR = get_case_study_dir(\"etfs\")\n", "\n", "prices = load_etfs()\n", "print(f\"Prices: {len(prices):,} rows, {prices['symbol'].n_unique()} assets\")" ] }, { "cell_type": "markdown", "id": "f50f04e9", "metadata": { "papermill": { "duration": 0.001858, "end_time": "2026-03-24T02:03:47.599637+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.597779+00:00", "status": "completed" } }, "source": [ "## CV Fold Setup\n", "\n", "We load the walk-forward CV splits from `setup.yaml` and add a holdout\n", "fold. All temporal models are fit per fold on training data only." ] }, { "cell_type": "code", "execution_count": 5, "id": "e0fc5cf0", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:47.604300Z", "iopub.status.busy": "2026-03-24T02:03:47.604142Z", "iopub.status.idle": "2026-03-24T02:03:47.611460Z", "shell.execute_reply": "2026-03-24T02:03:47.610856Z" }, "papermill": { "duration": 0.010304, "end_time": "2026-03-24T02:03:47.611785+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.601481+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SPY: 5,010 observations (2006-02-02 to 2025-12-31)\n" ] } ], "source": [ "# Generate CV splits\n", "cv_splits = generate_cv_splits(prices, case_study_id=\"etfs\", label_buffer=\"21D\")\n", "eval_config = load_evaluation_config(\"etfs\")\n", "\n", "# Add holdout fold: fit on all pre-holdout data, extract for holdout period\n", "holdout_start = str(eval_config[\"holdout_start\"])\n", "holdout_end = str(eval_config.get(\"holdout_end\", prices[\"timestamp\"].max()))\n", "\n", "# The last CV fold's val_end is the boundary before holdout\n", "# For holdout, train on everything up to holdout_start\n", "holdout_fold = {\n", " \"fold\": len(cv_splits),\n", " \"train_start\": cv_splits[0][\"train_start\"],\n", " \"train_end\": holdout_start,\n", " \"val_start\": holdout_start,\n", " \"val_end\": str(holdout_end),\n", "}\n", "all_folds = cv_splits + [holdout_fold]\n", "\n", "n_cv = len(cv_splits)\n", "n_total = len(all_folds)\n", "print(f\"CV folds: {n_cv}, plus holdout fold (fold {n_cv})\")\n", "for fold in all_folds:\n", " label = \"HOLDOUT\" if fold[\"fold\"] == n_cv else f\"Fold {fold['fold']}\"\n", " print(\n", " f\" {label}: train {fold['train_start']}..{fold['train_end']}, \"\n", " f\"val {fold['val_start']}..{fold['val_end']}\"\n", " )" ] }, { "cell_type": "markdown", "id": "2e887a92", "metadata": {}, "source": [ "## Part 1: HMM Regime Detection\n", "\n", "We fit a 2-state Gaussian HMM on SPY returns + volatility **per fold**.\n", "The aggregate market drives regime classification; individual ETFs\n", "inherit it.\n", "\n", "### Why SPY Only\n", "\n", "Using a single aggregate proxy (SPY) rather than per-asset HMMs:\n", "- Avoids overfitting 100 independent HMMs\n", "- Regime is a market-level phenomenon (risk-on/risk-off)\n", "- All ETFs inherit the same regime state, ensuring cross-sectional consistency" ] }, { "cell_type": "code", "execution_count": null, "id": "43a97cfd", "metadata": {}, "outputs": [], "source": [ "spy_full = (\n", " prices.filter(pl.col(\"symbol\") == \"SPY\")\n", " .sort(\"timestamp\")\n", " .with_columns(\n", " log_ret=(pl.col(\"close\").log().diff() * 100),\n", " vol_21d=(pl.col(\"close\").log().diff().rolling_std(window_size=21) * 100 * np.sqrt(252)),\n", " )\n", " .drop_nulls()\n", ")\n", "\n", "print(\n", " f\"SPY: {len(spy_full):,} observations ({spy_full['timestamp'].min()} to {spy_full['timestamp'].max()})\"\n", ")" ] }, { "cell_type": "markdown", "id": "8e19853a", "metadata": { "lines_to_next_cell": 2, "papermill": { "duration": 0.00182, "end_time": "2026-03-24T02:03:47.615604+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.613784+00:00", "status": "completed" } }, "source": [ "### K-Means-Seeded HMM Fitting\n", "\n", "K-means clustering provides better initial emission parameters than random\n", "initialization, reducing sensitivity to local optima." ] }, { "cell_type": "code", "execution_count": 6, "id": "e6d8c796", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:47.620042Z", "iopub.status.busy": "2026-03-24T02:03:47.619854Z", "iopub.status.idle": "2026-03-24T02:03:47.622581Z", "shell.execute_reply": "2026-03-24T02:03:47.622155Z" }, "lines_to_next_cell": 2, "papermill": { "duration": 0.005436, "end_time": "2026-03-24T02:03:47.622855+00:00", "exception": false, "start_time": "2026-03-24T02:03:47.617419+00:00", "status": "completed" } }, "outputs": [], "source": [ "def fit_hmm_kmeans_init(X: np.ndarray, n_states: int, random_state: int = 42) -> GaussianHMM:\n", " \"\"\"Fit HMM with k-means-seeded initialization.\"\"\"\n", " kmeans = KMeans(n_clusters=n_states, random_state=random_state, n_init=10)\n", " kmeans.fit(X)\n", "\n", " model = GaussianHMM(\n", " n_components=n_states,\n", " covariance_type=\"full\",\n", " n_iter=200,\n", " random_state=random_state,\n", " init_params=\"st\", # Only init startprob and transmat\n", " )\n", "\n", " model.means_ = kmeans.cluster_centers_\n", " model.covars_ = np.array(\n", " [np.cov(X[kmeans.labels_ == k].T) + np.eye(X.shape[1]) * 1e-6 for k in range(n_states)]\n", " )\n", "\n", " model.fit(X)\n", " return model" ] }, { "cell_type": "markdown", "id": "e79e549c", "metadata": { "lines_to_next_cell": 2, "papermill": { "duration": 0.006122, "end_time": "2026-03-24T02:03:50.033381+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.027259+00:00", "status": "completed" } }, "source": [ "### Label Switching Prevention\n", "\n", "Sort states by variance (ascending) so State 0 is always \"low volatility\"\n", "(calm) and State 1 is always \"high volatility\" (stressed)." ] }, { "cell_type": "code", "execution_count": 8, "id": "5ec3601c", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.046736Z", "iopub.status.busy": "2026-03-24T02:03:50.046556Z", "iopub.status.idle": "2026-03-24T02:03:50.054166Z", "shell.execute_reply": "2026-03-24T02:03:50.053698Z" }, "lines_to_next_cell": 2, "papermill": { "duration": 0.017851, "end_time": "2026-03-24T02:03:50.057408+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.039557+00:00", "status": "completed" } }, "outputs": [], "source": [ "def sort_states_by_variance(model: GaussianHMM) -> np.ndarray:\n", " \"\"\"Sort HMM states by variance (ascending) for consistent labeling.\"\"\"\n", " variances = np.array([np.trace(model.covars_[k]) for k in range(model.n_components)])\n", " return np.argsort(variances) # Low vol first\n", "\n", "\n", "def relabel_states(states: np.ndarray, probs: np.ndarray, order: np.ndarray) -> tuple:\n", " \"\"\"Relabel states according to the given order.\"\"\"\n", " inv_order = np.argsort(order)\n", " new_states = inv_order[states]\n", " new_probs = probs[:, order]\n", " return new_states, new_probs" ] }, { "cell_type": "markdown", "id": "4e78ae17", "metadata": { "lines_to_next_cell": 2, "papermill": { "duration": 0.006123, "end_time": "2026-03-24T02:03:50.068749+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.062626+00:00", "status": "completed" } }, "source": [ "### Filtered Probabilities (No Look-Ahead)\n", "\n", "In production, we must use **filtered** probabilities $P(z_t | x_{1:t})$\n", "which condition only on past and present observations. hmmlearn's\n", "`predict_proba()` returns **smoothed** probabilities $P(z_t | x_{1:T})$\n", "which use future data and would introduce look-ahead bias." ] }, { "cell_type": "code", "execution_count": 9, "id": "6fe164e7", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.082404Z", "iopub.status.busy": "2026-03-24T02:03:50.082220Z", "iopub.status.idle": "2026-03-24T02:03:50.089465Z", "shell.execute_reply": "2026-03-24T02:03:50.088796Z" }, "lines_to_next_cell": 2, "papermill": { "duration": 0.014982, "end_time": "2026-03-24T02:03:50.089814+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.074832+00:00", "status": "completed" } }, "outputs": [], "source": [ "def compute_filtered_probs(model: GaussianHMM, X: np.ndarray) -> np.ndarray:\n", " \"\"\"Compute filtered probabilities P(state_t | obs_{1:t}).\n", "\n", " Uses the forward algorithm, then normalizes.\n", " \"\"\"\n", " framelogprob = model._compute_log_likelihood(X)\n", "\n", " n_samples = X.shape[0]\n", " n_components = model.n_components\n", "\n", " log_startprob = np.log(model.startprob_ + 1e-300)\n", " log_transmat = np.log(model.transmat_ + 1e-300)\n", "\n", " # Forward pass (log-domain for numerical stability)\n", " fwdlattice = np.zeros((n_samples, n_components))\n", " fwdlattice[0] = log_startprob + framelogprob[0]\n", "\n", " for t in range(1, n_samples):\n", " for j in range(n_components):\n", " fwdlattice[t, j] = framelogprob[t, j] + np.logaddexp.reduce(\n", " fwdlattice[t - 1] + log_transmat[:, j]\n", " )\n", "\n", " # Normalize to get probabilities\n", " log_normalizer = np.logaddexp.reduce(fwdlattice, axis=1, keepdims=True)\n", " filtered = np.exp(fwdlattice - log_normalizer)\n", "\n", " return filtered" ] }, { "cell_type": "markdown", "id": "007de7e7", "metadata": { "lines_to_next_cell": 2 }, "source": [ "### Derive Regime Features from Filtered Probabilities\n", "\n", "From filtered probabilities, derive three feature types:\n", "1. **regime_prob_stress**: Filtered probability of being in the high-vol state\n", "2. **regime_transition**: Absolute change in stress probability (detects regime shifts)\n", "3. **regime_duration**: Days since last regime change (persistence indicator)" ] }, { "cell_type": "code", "execution_count": 10, "id": "1f4b9131", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.096057Z", "iopub.status.busy": "2026-03-24T02:03:50.095887Z", "iopub.status.idle": "2026-03-24T02:03:50.175498Z", "shell.execute_reply": "2026-03-24T02:03:50.175099Z" }, "papermill": { "duration": 0.084884, "end_time": "2026-03-24T02:03:50.177857+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.092973+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "State 0 (Low-Vol): 69.0% of time, mean ret=0.072%, vol=10.9%\n", "State 1 (High-Vol): 31.0% of time, mean ret=-0.030%, vol=27.3%\n" ] } ], "source": [ "def derive_regime_features(\n", " timestamps: pl.Series,\n", " filtered_probs: np.ndarray,\n", " states: np.ndarray,\n", " order: np.ndarray,\n", ") -> pl.DataFrame:\n", " \"\"\"Derive regime features from HMM output for a single fold window.\"\"\"\n", " states_sorted, filtered_sorted = relabel_states(states, filtered_probs, order)\n", "\n", " regime_prob_stress = filtered_sorted[:, 1] # P(high-vol state)\n", "\n", " # Transition: absolute 1-day change in stress probability\n", " regime_transition = np.abs(np.diff(regime_prob_stress, prepend=regime_prob_stress[0]))\n", "\n", " # Duration: days since last regime change\n", " regime_duration = np.zeros(len(states_sorted))\n", " current_run = 0\n", " for i in range(len(states_sorted)):\n", " if i == 0 or states_sorted[i] != states_sorted[i - 1]:\n", " current_run = 1\n", " else:\n", " current_run += 1\n", " regime_duration[i] = current_run\n", "\n", " return pl.DataFrame(\n", " {\n", " \"timestamp\": timestamps,\n", " \"regime_prob_stress\": regime_prob_stress,\n", " \"regime_transition\": regime_transition,\n", " \"regime_log_duration\": np.log1p(regime_duration),\n", " }\n", " )" ] }, { "cell_type": "markdown", "id": "2b5896fe", "metadata": {}, "source": [ "### Illustrative Full-Sample HMM Fit\n", "\n", "Before the per-fold loop, we fit one full-sample HMM for visualization\n", "purposes only. This shows readers the regime overlay on SPY and validates\n", "the methodology. **These features are NOT used in the saved output.**" ] }, { "cell_type": "code", "execution_count": 11, "id": "2d0e2a06", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.194752Z", "iopub.status.busy": "2026-03-24T02:03:50.194620Z", "iopub.status.idle": "2026-03-24T02:03:50.201758Z", "shell.execute_reply": "2026-03-24T02:03:50.200960Z" }, "papermill": { "duration": 0.018249, "end_time": "2026-03-24T02:03:50.202188+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.183939+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Transition matrix (rows=from, cols=to):\n", " Calm Stress\n", " Calm 0.9925 0.0075\n", " Stress 0.0166 0.9834\n", "\n", "Mean regime duration: calm=133d, stress=60d\n" ] } ], "source": [ "N_STATES = 2\n", "X_full = spy_full.select([\"log_ret\", \"vol_21d\"]).to_numpy()\n", "\n", "best_model_illustrative = None\n", "best_ll = -np.inf\n", "\n", "for seed in range(N_RESTARTS):\n", " try:\n", " model = fit_hmm_kmeans_init(X_full, n_states=N_STATES, random_state=seed)\n", " ll = model.score(X_full)\n", " if ll > best_ll:\n", " best_ll = ll\n", " best_model_illustrative = model\n", " except Exception:\n", " continue\n", "\n", "print(f\"Illustrative full-sample HMM: best log-likelihood = {best_ll:.1f}\")\n", "\n", "order_illustrative = sort_states_by_variance(best_model_illustrative)\n", "filtered_illustrative = compute_filtered_probs(best_model_illustrative, X_full)\n", "states_illustrative = best_model_illustrative.predict(X_full)\n", "states_sorted_ill, _ = relabel_states(\n", " states_illustrative, filtered_illustrative, order_illustrative\n", ")\n", "\n", "for k in range(N_STATES):\n", " mask = states_sorted_ill == k\n", " label = \"Low-Vol\" if k == 0 else \"High-Vol\"\n", " mean_ret = X_full[mask, 0].mean()\n", " mean_vol = X_full[mask, 1].mean()\n", " pct = mask.mean()\n", " print(f\"State {k} ({label}): {pct:.1%} of time, mean ret={mean_ret:.3f}%, vol={mean_vol:.1f}%\")" ] }, { "cell_type": "markdown", "id": "c9028209", "metadata": { "papermill": { "duration": 0.005464, "end_time": "2026-03-24T02:03:50.215275+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.209811+00:00", "status": "completed" } }, "source": [ "### Regime Overlay on SPY (Illustrative)\n", "\n", "Shading stressed periods on SPY cumulative returns shows whether the HMM\n", "captures known market episodes (GFC, COVID, 2022 rate shock)." ] }, { "cell_type": "code", "execution_count": 12, "id": "f1e8a49b", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.233293Z", "iopub.status.busy": "2026-03-24T02:03:50.233169Z", "iopub.status.idle": "2026-03-24T02:03:50.301461Z", "shell.execute_reply": "2026-03-24T02:03:50.300716Z" }, "papermill": { "duration": 0.077028, "end_time": "2026-03-24T02:03:50.302153+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.225125+00:00", "status": "completed" } }, "outputs": [], "source": [ "spy_cum = spy_full.with_columns(cum_ret=(pl.col(\"close\") / pl.col(\"close\").first() - 1) * 100)\n", "\n", "fig_regime, ax = plt.subplots(figsize=(12, 5))\n", "dates = spy_cum[\"timestamp\"].to_numpy()\n", "cum_ret = spy_cum[\"cum_ret\"].to_numpy()\n", "ax.plot(dates, cum_ret, linewidth=0.8, color=\"0.3\")\n", "\n", "# Shade stressed periods\n", "stressed = states_sorted_ill == 1\n", "in_stress = False\n", "start = None\n", "for i in range(len(stressed)):\n", " if stressed[i] and not in_stress:\n", " start = dates[i]\n", " in_stress = True\n", " elif not stressed[i] and in_stress:\n", " ax.axvspan(start, dates[i], alpha=0.15, color=\"red\", linewidth=0)\n", " in_stress = False\n", "if in_stress:\n", " ax.axvspan(start, dates[-1], alpha=0.15, color=\"red\", linewidth=0)\n", "\n", "ax.set_ylabel(\"Cumulative Return (%)\")\n", "ax.set_title(\"SPY with HMM Stress Regimes (full-sample illustrative)\")\n", "fig_regime.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6a292625", "metadata": { "papermill": { "duration": 0.002015, "end_time": "2026-03-24T02:03:50.307431+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.305416+00:00", "status": "completed" } }, "source": [ "### Per-Fold HMM Fitting\n", "\n", "For each fold, we fit the HMM on **training data only**, then apply the\n", "forward algorithm to the full fold window (train_start through val_end)\n", "for filtered probabilities. The model parameters $\\theta$ are estimated\n", "exclusively from training observations, eliminating parameter-level\n", "look-ahead." ] }, { "cell_type": "code", "execution_count": null, "id": "c847797c", "metadata": {}, "outputs": [], "source": [ "hmm_fold_results = []\n", "\n", "for fold in all_folds:\n", " fold_idx = fold[\"fold\"]\n", " train_start, train_end = fold[\"train_start\"], fold[\"train_end\"]\n", " val_end = fold[\"val_end\"]\n", "\n", " # Training data: fit HMM parameters\n", " spy_train = spy_full.filter(\n", " (pl.col(\"timestamp\") >= pl.lit(train_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") < pl.lit(train_end).cast(pl.Date))\n", " )\n", " X_train = spy_train.select([\"log_ret\", \"vol_21d\"]).to_numpy()\n", "\n", " if len(X_train) < 252:\n", " print(\n", " f\" Fold {fold_idx}: insufficient SPY training data ({len(X_train)} obs), skipping HMM\"\n", " )\n", " continue\n", "\n", " # Fit HMM with multiple restarts on training data\n", " best_fold_model = None\n", " best_fold_ll = -np.inf\n", " for seed in range(N_RESTARTS):\n", " try:\n", " m = fit_hmm_kmeans_init(X_train, n_states=N_STATES, random_state=seed)\n", " ll = m.score(X_train)\n", " if ll > best_fold_ll:\n", " best_fold_ll = ll\n", " best_fold_model = m\n", " except Exception:\n", " continue\n", "\n", " if best_fold_model is None:\n", " print(f\" Fold {fold_idx}: HMM fitting failed\")\n", " continue\n", "\n", " order = sort_states_by_variance(best_fold_model)\n", "\n", " # Full fold window (train_start through val_end) for filtered probs\n", " spy_fold = spy_full.filter(\n", " (pl.col(\"timestamp\") >= pl.lit(train_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") <= pl.lit(val_end).cast(pl.Date))\n", " )\n", " X_fold = spy_fold.select([\"log_ret\", \"vol_21d\"]).to_numpy()\n", "\n", " filtered = compute_filtered_probs(best_fold_model, X_fold)\n", " states = best_fold_model.predict(X_fold)\n", "\n", " regime_df = derive_regime_features(spy_fold[\"timestamp\"], filtered, states, order)\n", " regime_df = regime_df.with_columns(pl.lit(fold_idx).alias(\"fold\"))\n", "\n", " hmm_fold_results.append(regime_df)\n", " print(\n", " f\" Fold {fold_idx}: HMM LL={best_fold_ll:.1f}, {len(regime_df)} dates, \"\n", " f\"stress={regime_df['regime_prob_stress'].mean():.3f}\"\n", " )\n", "\n", "hmm_features = (\n", " pl.concat(hmm_fold_results)\n", " if hmm_fold_results\n", " else pl.DataFrame(schema={\"timestamp\": pl.Date, \"fold\": pl.Int64})\n", ")\n", "n_hmm_folds = hmm_features[\"fold\"].n_unique() if len(hmm_features) > 0 else 0\n", "print(f\"\\nHMM features: {len(hmm_features):,} rows across {n_hmm_folds} folds\")" ] }, { "cell_type": "markdown", "id": "0b496f76", "metadata": { "papermill": { "duration": 0.001928, "end_time": "2026-03-24T02:03:50.311483+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.309555+00:00", "status": "completed" } }, "source": [ "## Part 2: Fractional Differencing (Per Fold)\n", "\n", "Fractional differencing preserves long-range memory while achieving\n", "stationarity. We apply fixed $d$ values by asset class to 10 reference\n", "ETFs spanning all major asset classes in our universe.\n", "\n", "**Why fixed $d$?** Using pre-specified $d$ values avoids parameter estimation\n", "lookahead entirely -- no data-dependent optimization, so the transform is\n", "purely mechanical. We still compute per fold so each fold window gets a\n", "clean series starting from its own `train_start`.\n", "\n", "| Asset Class | Reference ETFs | $d$ |\n", "|-----------------|---------------|-----|\n", "| US Equities | SPY, QQQ, IWM | 0.4 |\n", "| Int'l Equities | EFA, EEM | 0.4 |\n", "| Fixed Income | TLT | 0.5 |\n", "| Gold | GLD | 0.4 |\n", "| Real Estate | VNQ | 0.4 |\n", "| High Yield | HYG | 0.5 |\n", "| Inv. Grade | LQD | 0.5 |" ] }, { "cell_type": "code", "execution_count": 14, "id": "7c5aae31", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.333922Z", "iopub.status.busy": "2026-03-24T02:03:50.333758Z", "iopub.status.idle": "2026-03-24T02:03:50.336571Z", "shell.execute_reply": "2026-03-24T02:03:50.336064Z" }, "papermill": { "duration": 0.005888, "end_time": "2026-03-24T02:03:50.336883+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.330995+00:00", "status": "completed" } }, "outputs": [], "source": [ "REFERENCE_ETFS = {\n", " \"SPY\": 0.4, # US large-cap equities\n", " \"QQQ\": 0.4, # US tech equities\n", " \"IWM\": 0.4, # US small-cap equities\n", " \"EFA\": 0.4, # International developed\n", " \"EEM\": 0.4, # Emerging markets\n", " \"TLT\": 0.5, # Long-term treasuries\n", " \"GLD\": 0.4, # Gold\n", " \"VNQ\": 0.4, # Real estate\n", " \"HYG\": 0.5, # High yield bonds\n", " \"LQD\": 0.5, # Investment grade bonds\n", "}" ] }, { "cell_type": "markdown", "id": "fdeb7b24", "metadata": { "papermill": { "duration": 0.001931, "end_time": "2026-03-24T02:03:50.340958+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.339027+00:00", "status": "completed" } }, "source": [ "### Per-Fold FFD Application\n", "\n", "For each fold, apply fractional differencing to the fold window\n", "(train_start through val_end). The FFD filter uses a fixed-width window\n", "of historical weights, so it only needs a warmup period at the start --\n", "no fitting is involved. Computing per fold ensures the warmup loss\n", "does not leak information across fold boundaries." ] }, { "cell_type": "code", "execution_count": 15, "id": "fb2891a0", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.345376Z", "iopub.status.busy": "2026-03-24T02:03:50.345289Z", "iopub.status.idle": "2026-03-24T02:03:50.515022Z", "shell.execute_reply": "2026-03-24T02:03:50.514577Z" }, "papermill": { "duration": 0.172458, "end_time": "2026-03-24T02:03:50.515343+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.342885+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " SPY (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " QQQ (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " IWM (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " EFA (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " EEM (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " TLT (d=0.5): 5,031 valid / 5,031 total (warmup=0)\n", " GLD (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " VNQ (d=0.4): 5,031 valid / 5,031 total (warmup=0)\n", " HYG (d=0.5): 4,713 valid / 4,713 total (warmup=0)\n", " LQD (d=0.5): 5,031 valid / 5,031 total (warmup=0)\n" ] } ], "source": [ "ffd_fold_results = []\n", "\n", "for fold in all_folds:\n", " fold_idx = fold[\"fold\"]\n", " train_start = fold[\"train_start\"]\n", " val_end = fold[\"val_end\"]\n", "\n", " fold_frames = []\n", " for symbol, d in REFERENCE_ETFS.items():\n", " etf = (\n", " prices.filter(\n", " (pl.col(\"symbol\") == symbol)\n", " & (pl.col(\"timestamp\") >= pl.lit(train_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") <= pl.lit(val_end).cast(pl.Date))\n", " )\n", " .sort(\"timestamp\")\n", " .select([\"timestamp\", \"close\"])\n", " )\n", "\n", " if len(etf) == 0:\n", " continue\n", "\n", " log_close = etf[\"close\"].log()\n", " ffd_series = ffdiff(log_close, d=d)\n", "\n", " ffd_df = pl.DataFrame(\n", " {\n", " \"timestamp\": etf[\"timestamp\"],\n", " f\"ffd_{symbol.lower()}\": ffd_series,\n", " }\n", " ).drop_nulls()\n", "\n", " fold_frames.append(ffd_df)\n", "\n", " if fold_frames:\n", " ffd_fold = fold_frames[0]\n", " for df in fold_frames[1:]:\n", " ffd_fold = ffd_fold.join(df, on=\"timestamp\", how=\"outer_coalesce\")\n", " ffd_fold = ffd_fold.sort(\"timestamp\").with_columns(pl.lit(fold_idx).alias(\"fold\"))\n", " ffd_fold_results.append(ffd_fold)\n", " ffd_cols = [c for c in ffd_fold.columns if c.startswith(\"ffd_\")]\n", " print(f\" Fold {fold_idx}: FFD {len(ffd_cols)} series, {len(ffd_fold):,} dates\")\n", " else:\n", " print(f\" Fold {fold_idx}: no FFD series produced (no qualifying ETFs)\")\n", "\n", "ffd_features = pl.concat(ffd_fold_results) if ffd_fold_results else pl.DataFrame()\n", "ffd_col_names = [c for c in ffd_features.columns if c.startswith(\"ffd_\")]\n", "print(f\"\\nFFD features: {len(ffd_col_names)} series, {len(ffd_features):,} rows across folds\")" ] }, { "cell_type": "markdown", "id": "dae312b1", "metadata": { "papermill": { "duration": 0.00256, "end_time": "2026-03-24T02:03:50.520153+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.517593+00:00", "status": "completed" } }, "source": [ "## Part 3: Per-ETF GARCH(1,1) (Per Fold)\n", "\n", "For each fold, fit GARCH(1,1) on each ETF's **training returns**, then\n", "use `model.fix(params)` to run the variance recursion on the full fold\n", "window (train through val_end) without re-estimating parameters. This\n", "is the **fit-then-filter** paradigm: parameters come from training only,\n", "but conditional volatility is computed for every date in the window.\n", "\n", "The `fix()` method applies frozen parameters to a new data series,\n", "producing a causal conditional volatility path $\\sigma_t$ that depends\n", "only on past returns given the frozen parameters." ] }, { "cell_type": "code", "execution_count": 16, "id": "8e987fff", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.527379Z", "iopub.status.busy": "2026-03-24T02:03:50.527269Z", "iopub.status.idle": "2026-03-24T02:03:50.546354Z", "shell.execute_reply": "2026-03-24T02:03:50.545327Z" }, "papermill": { "duration": 0.02357, "end_time": "2026-03-24T02:03:50.547016+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.523446+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "FFD features: 10 series, 5,031 rows\n" ] } ], "source": [ "all_symbols = sorted(prices[\"symbol\"].unique().to_list())\n", "if MAX_SYMBOLS > 0:\n", " all_symbols = all_symbols[:MAX_SYMBOLS]\n", "\n", "print(f\"Fitting GARCH(1,1) on {len(all_symbols)} ETFs across {len(all_folds)} folds...\")" ] }, { "cell_type": "code", "execution_count": null, "id": "252c615f", "metadata": {}, "outputs": [], "source": [ "def fit_garch_fold(\n", " prices_df: pl.DataFrame,\n", " symbols: list[str],\n", " fold: dict,\n", " min_obs: int,\n", ") -> pl.DataFrame:\n", " \"\"\"Fit GARCH(1,1) per symbol for one fold.\n", "\n", " Parameters\n", " ----------\n", " prices_df : pl.DataFrame\n", " Full prices panel with timestamp, symbol, close columns.\n", " symbols : list[str]\n", " Symbols to fit.\n", " fold : dict\n", " Fold dict with train_start, train_end, val_end keys.\n", " min_obs : int\n", " Minimum training observations for GARCH fit.\n", "\n", " Returns\n", " -------\n", " pl.DataFrame\n", " Conditional volatility for the full fold window with fold column.\n", " \"\"\"\n", " fold_idx = fold[\"fold\"]\n", " train_start = fold[\"train_start\"]\n", " train_end = fold[\"train_end\"]\n", " val_end = fold[\"val_end\"]\n", "\n", " results = []\n", " n_success = 0\n", " n_fail = 0\n", "\n", " for sym in symbols:\n", " # Get full fold window data (train_start through val_end)\n", " sym_data = (\n", " prices_df.filter(\n", " (pl.col(\"symbol\") == sym)\n", " & (pl.col(\"timestamp\") >= pl.lit(train_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") <= pl.lit(val_end).cast(pl.Date))\n", " )\n", " .sort(\"timestamp\")\n", " .with_columns(ret=pl.col(\"close\").pct_change())\n", " .drop_nulls(subset=[\"ret\"])\n", " )\n", "\n", " # Training returns only\n", " train_data = sym_data.filter(pl.col(\"timestamp\") < pl.lit(train_end).cast(pl.Date))\n", "\n", " if len(train_data) < min_obs:\n", " n_fail += 1\n", " continue\n", "\n", " train_returns_pct = (train_data[\"ret\"] * 100).to_numpy()\n", "\n", " try:\n", " # Fit on training data\n", " train_model = arch_model(\n", " train_returns_pct,\n", " mean=\"Constant\",\n", " vol=\"GARCH\",\n", " p=1,\n", " q=1,\n", " dist=\"Normal\",\n", " )\n", " train_result = train_model.fit(disp=\"off\", show_warning=False)\n", "\n", " # Apply frozen parameters to full fold window\n", " full_returns_pct = (sym_data[\"ret\"] * 100).to_numpy()\n", " full_model = arch_model(\n", " full_returns_pct,\n", " mean=\"Constant\",\n", " vol=\"GARCH\",\n", " p=1,\n", " q=1,\n", " dist=\"Normal\",\n", " )\n", " filtered = full_model.fix(train_result.params)\n", "\n", " # Annualized conditional vol (input is in % daily)\n", " cond_vol_ann = filtered.conditional_volatility * np.sqrt(252) / 100\n", "\n", " sym_result = pl.DataFrame(\n", " {\n", " \"timestamp\": sym_data[\"timestamp\"].to_list(),\n", " \"symbol\": [sym] * len(sym_data),\n", " \"garch_cond_vol\": cond_vol_ann,\n", " \"fold\": [fold_idx] * len(sym_data),\n", " }\n", " ).drop_nulls()\n", "\n", " if len(sym_result) > 0:\n", " results.append(sym_result)\n", " n_success += 1\n", " except Exception:\n", " n_fail += 1\n", "\n", " print(f\" Fold {fold_idx} GARCH: {n_success}/{len(symbols)} fitted, {n_fail} failed/skipped\")\n", " return pl.concat(results) if results else pl.DataFrame()" ] }, { "cell_type": "code", "execution_count": null, "id": "eb6992ce", "metadata": {}, "outputs": [], "source": [ "garch_fold_results = []\n", "\n", "for fold in all_folds:\n", " garch_fold = fit_garch_fold(prices, all_symbols, fold, GARCH_MIN_OBS)\n", " if len(garch_fold) > 0:\n", " garch_fold_results.append(garch_fold)\n", "\n", "garch_df = pl.concat(garch_fold_results) if garch_fold_results else pl.DataFrame()\n", "garch_cols = [\"garch_cond_vol\"]\n", "\n", "if len(garch_df) > 0:\n", " n_syms = garch_df[\"symbol\"].n_unique()\n", " print(\n", " f\"\\nGARCH features: {len(garch_df):,} rows, {n_syms} assets, {garch_df['fold'].n_unique()} folds\"\n", " )\n", " print(\n", " f\" Conditional vol: mean={garch_df['garch_cond_vol'].mean():.3f}, \"\n", " f\"std={garch_df['garch_cond_vol'].std():.3f}\"\n", " )" ] }, { "cell_type": "markdown", "id": "f3becf69", "metadata": { "papermill": { "duration": 0.002708, "end_time": "2026-03-24T02:03:50.553562+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.550854+00:00", "status": "completed" } }, "source": [ "## Part 4: Combine and Broadcast to Per-Fold Panel\n", "\n", "HMM regime features and FFD features are date-level (one value per day\n", "shared by all ETFs). GARCH features are per-asset. For each fold, we\n", "broadcast date-level features to all symbols and join with per-asset\n", "GARCH features, producing a `(fold, timestamp, symbol)` panel." ] }, { "cell_type": "code", "execution_count": 17, "id": "1e898785", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:50.558106Z", "iopub.status.busy": "2026-03-24T02:03:50.557983Z", "iopub.status.idle": "2026-03-24T02:03:50.562863Z", "shell.execute_reply": "2026-03-24T02:03:50.562518Z" }, "papermill": { "duration": 0.007487, "end_time": "2026-03-24T02:03:50.563099+00:00", "exception": false, "start_time": "2026-03-24T02:03:50.555612+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fitting GARCH(1,1) on 100 ETFs...\n" ] } ], "source": [ "fold_panels = []\n", "\n", "for fold in all_folds:\n", " fold_idx = fold[\"fold\"]\n", " train_start = fold[\"train_start\"]\n", " val_end = fold[\"val_end\"]\n", "\n", " # Get panel skeleton for this fold's date range\n", " fold_skeleton = (\n", " prices.filter(\n", " (pl.col(\"timestamp\") >= pl.lit(train_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") <= pl.lit(val_end).cast(pl.Date))\n", " )\n", " .select([\"timestamp\", \"symbol\"])\n", " .unique()\n", " .sort([\"timestamp\", \"symbol\"])\n", " .with_columns(pl.lit(fold_idx).alias(\"fold\"))\n", " )\n", "\n", " # Get date-level features for this fold\n", " fold_hmm = (\n", " hmm_features.filter(pl.col(\"fold\") == fold_idx).drop(\"fold\")\n", " if len(hmm_features) > 0\n", " else pl.DataFrame()\n", " )\n", " fold_ffd = (\n", " ffd_features.filter(pl.col(\"fold\") == fold_idx).drop(\"fold\")\n", " if len(ffd_features) > 0\n", " else pl.DataFrame()\n", " )\n", "\n", " # Combine date-level features\n", " if len(fold_hmm) > 0 and len(fold_ffd) > 0:\n", " date_level = fold_hmm.join(fold_ffd, on=\"timestamp\", how=\"outer_coalesce\")\n", " elif len(fold_hmm) > 0:\n", " date_level = fold_hmm\n", " elif len(fold_ffd) > 0:\n", " date_level = fold_ffd\n", " else:\n", " date_level = pl.DataFrame()\n", "\n", " # Broadcast date-level features to all assets\n", " if len(date_level) > 0:\n", " panel = fold_skeleton.join(date_level, on=\"timestamp\", how=\"left\")\n", " else:\n", " panel = fold_skeleton\n", "\n", " # Join per-asset GARCH features\n", " if len(garch_df) > 0:\n", " fold_garch = garch_df.filter(pl.col(\"fold\") == fold_idx).drop(\"fold\")\n", " panel = panel.join(fold_garch, on=[\"timestamp\", \"symbol\"], how=\"left\")\n", "\n", " fold_panels.append(panel)\n", "\n", "temporal = pl.concat(fold_panels).sort([\"fold\", \"timestamp\", \"symbol\"])\n", "\n", "temporal_cols = [c for c in temporal.columns if c not in (\"timestamp\", \"symbol\", \"fold\")]\n", "print(f\"Combined model-based features: {len(temporal_cols)} columns, {len(temporal):,} rows\")\n", "print(f\" Assets: {temporal['symbol'].n_unique()}, Folds: {temporal['fold'].n_unique()}\")" ] }, { "cell_type": "markdown", "id": "4fec8bad", "metadata": { "papermill": { "duration": 0.002072, "end_time": "2026-03-24T02:03:52.465075+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.463003+00:00", "status": "completed" } }, "source": [ "### Quality Check\n", "\n", "Verify that temporal features have reasonable coverage and distributions." ] }, { "cell_type": "code", "execution_count": 21, "id": "5e87a8dc", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:52.469693Z", "iopub.status.busy": "2026-03-24T02:03:52.469612Z", "iopub.status.idle": "2026-03-24T02:03:52.480274Z", "shell.execute_reply": "2026-03-24T02:03:52.479949Z" }, "papermill": { "duration": 0.013348, "end_time": "2026-03-24T02:03:52.480584+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.467236+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " regime_prob_stress : 469,423/470,662 valid, mean=0.3186, std=0.4508\n", " regime_transition : 469,423/470,662 valid, mean=0.0199, std=0.0803\n", " regime_log_duration : 469,423/470,662 valid, mean=4.2489, std=1.3018\n", " ffd_spy : 470,662/470,662 valid, mean=0.2299, std=0.0967\n", " ffd_qqq : 470,662/470,662 valid, mean=0.2126, std=0.0785\n", " ffd_iwm : 470,662/470,662 valid, mean=0.2014, std=0.0876\n", " ffd_efa : 470,662/470,662 valid, mean=0.1664, std=0.0807\n", " ffd_eem : 470,662/470,662 valid, mean=0.1490, std=0.0720\n", " ffd_tlt : 470,662/470,662 valid, mean=0.0936, std=0.0655\n", " ffd_gld : 470,662/470,662 valid, mean=0.2111, std=0.0878\n", " ffd_vnq : 470,662/470,662 valid, mean=0.1685, std=0.0757\n", " ffd_hyg : 449,008/470,662 valid, mean=0.0867, std=0.0669\n", " ffd_lqd : 470,662/470,662 valid, mean=0.0934, std=0.0646\n", " garch_cond_vol : 470,562/470,662 valid, mean=0.1891, std=0.1325\n" ] } ], "source": [ "for col in temporal_cols:\n", " valid = temporal[col].drop_nulls().len()\n", " total = len(temporal)\n", " mean = temporal[col].drop_nulls().mean()\n", " std = temporal[col].drop_nulls().std()\n", " print(f\" {col:30s}: {valid:,}/{total:,} valid, mean={mean:.4f}, std={std:.4f}\")" ] }, { "cell_type": "markdown", "id": "6342476d", "metadata": {}, "source": [ "### Per-Fold Feature Stability\n", "\n", "Check that HMM and GARCH features are stable across folds (model\n", "parameters should not drift drastically with overlapping training windows)." ] }, { "cell_type": "code", "execution_count": null, "id": "9bf3d031", "metadata": {}, "outputs": [], "source": [ "print(\"\\nPer-fold feature means:\")\n", "fold_summary = (\n", " temporal.group_by(\"fold\")\n", " .agg([pl.col(c).mean().alias(f\"mean_{c}\") for c in temporal_cols])\n", " .sort(\"fold\")\n", ")\n", "print(fold_summary)" ] }, { "cell_type": "markdown", "id": "b2cbdaa8", "metadata": { "papermill": { "duration": 0.002126, "end_time": "2026-03-24T02:03:52.484919+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.482793+00:00", "status": "completed" } }, "source": [ "## Save Artifacts" ] }, { "cell_type": "code", "execution_count": 22, "id": "1e44e1fe", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:52.489483Z", "iopub.status.busy": "2026-03-24T02:03:52.489397Z", "iopub.status.idle": "2026-03-24T02:03:52.517320Z", "shell.execute_reply": "2026-03-24T02:03:52.516978Z" }, "papermill": { "duration": 0.03081, "end_time": "2026-03-24T02:03:52.517722+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.486912+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved: /tmp/ml4t-test-output/etfs/features/model_based.parquet (470,662 rows, 14 features)\n" ] } ], "source": [ "FEATURES_DIR = CASE_DIR / \"features\"\n", "\n", "FEATURES_DIR.mkdir(parents=True, exist_ok=True)\n", "temporal.write_parquet(FEATURES_DIR / \"model_based.parquet\")\n", "print(\n", " f\"Saved: {FEATURES_DIR / 'model_based.parquet'} \"\n", " f\"({len(temporal):,} rows, {len(temporal_cols)} features + fold column)\"\n", ")" ] }, { "cell_type": "markdown", "id": "eb8549b5", "metadata": { "papermill": { "duration": 0.002027, "end_time": "2026-03-24T02:03:52.521952+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.519925+00:00", "status": "completed" } }, "source": [ "## Incremental Evaluation\n", "\n", "Evaluate feature quality using **validation-period data only** from each\n", "fold. We compute cross-sectional Spearman IC for per-asset features and\n", "time-series IC for date-level features." ] }, { "cell_type": "code", "execution_count": 23, "id": "74abfbab", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:52.526441Z", "iopub.status.busy": "2026-03-24T02:03:52.526364Z", "iopub.status.idle": "2026-03-24T02:03:52.562802Z", "shell.execute_reply": "2026-03-24T02:03:52.562474Z" }, "papermill": { "duration": 0.039265, "end_time": "2026-03-24T02:03:52.563204+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.523939+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Evaluation panel: 468,562 rows, 100 assets\n" ] } ], "source": [ "labels = pl.read_parquet(CASE_DIR / \"labels\" / \"fwd_ret_21d.parquet\")\n", "label_col = \"fwd_ret_21d\"\n", "\n", "# Only evaluate on validation periods (not training periods)\n", "val_rows = []\n", "for fold in all_folds:\n", " fold_idx = fold[\"fold\"]\n", " val_start = fold[\"val_start\"]\n", " val_end = fold[\"val_end\"]\n", " fold_val = temporal.filter(\n", " (pl.col(\"fold\") == fold_idx)\n", " & (pl.col(\"timestamp\") >= pl.lit(val_start).cast(pl.Date))\n", " & (pl.col(\"timestamp\") <= pl.lit(val_end).cast(pl.Date))\n", " )\n", " val_rows.append(fold_val)\n", "\n", "val_temporal = pl.concat(val_rows) if val_rows else temporal\n", "\n", "eval_df = val_temporal.join(labels, on=[\"timestamp\", \"symbol\"], how=\"inner\").drop_nulls(\n", " subset=[label_col]\n", ")\n", "print(\n", " f\"Evaluation panel (val periods only): {len(eval_df):,} rows, {eval_df['symbol'].n_unique()} assets\"\n", ")" ] }, { "cell_type": "markdown", "id": "61a3a148", "metadata": { "papermill": { "duration": 0.002063, "end_time": "2026-03-24T02:03:52.567487+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.565424+00:00", "status": "completed" } }, "source": [ "### Cross-Sectional IC for Per-Asset Features" ] }, { "cell_type": "code", "execution_count": 24, "id": "0e0428ec", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:52.572018Z", "iopub.status.busy": "2026-03-24T02:03:52.571930Z", "iopub.status.idle": "2026-03-24T02:03:53.663229Z", "shell.execute_reply": "2026-03-24T02:03:53.662789Z" }, "papermill": { "duration": 1.094134, "end_time": "2026-03-24T02:03:53.663620+00:00", "exception": false, "start_time": "2026-03-24T02:03:52.569486+00:00", "status": "completed" } }, "outputs": [], "source": [ "temporal_ic = {}\n", "\n", "# Per-asset features: cross-sectional IC (rank correlation per date)\n", "per_asset_features = [c for c in temporal_cols if c in garch_cols]\n", "date_level_features = [c for c in temporal_cols if c not in garch_cols]\n", "\n", "for feat in per_asset_features:\n", " ic_series = (\n", " eval_df.filter(pl.col(feat).is_not_null())\n", " .with_columns(\n", " pl.col(feat).rank(method=\"average\").over(\"timestamp\").alias(\"_feat_rank\"),\n", " pl.col(label_col).rank(method=\"average\").over(\"timestamp\").alias(\"_label_rank\"),\n", " )\n", " .group_by(\"timestamp\")\n", " .agg(pl.corr(\"_feat_rank\", \"_label_rank\").alias(\"ic\"))\n", " .sort(\"timestamp\")\n", " )\n", " ics = ic_series[\"ic\"].drop_nulls().to_numpy()\n", " if len(ics) > 50:\n", " temporal_ic[feat] = {\n", " \"ic\": float(np.nanmean(ics)),\n", " \"t_stat\": float(np.nanmean(ics) / (np.nanstd(ics) / np.sqrt(len(ics)))),\n", " \"p_value\": None,\n", " \"bootstrap_std\": float(np.nanstd(ics)),\n", " }" ] }, { "cell_type": "markdown", "id": "4d788e06", "metadata": { "papermill": { "duration": 0.002077, "end_time": "2026-03-24T02:03:53.667948+00:00", "exception": false, "start_time": "2026-03-24T02:03:53.665871+00:00", "status": "completed" } }, "source": [ "### Time-Series IC for Date-Level Features\n", "\n", "Date-level features (HMM, FFD) are identical across symbols on each date,\n", "so cross-sectional IC is zero by construction. We evaluate them via\n", "time-series correlation with the cross-sectional average return." ] }, { "cell_type": "code", "execution_count": 25, "id": "b3a443e1", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:03:53.672483Z", "iopub.status.busy": "2026-03-24T02:03:53.672395Z", "iopub.status.idle": "2026-03-24T02:04:29.183354Z", "shell.execute_reply": "2026-03-24T02:04:29.182944Z" }, "papermill": { "duration": 35.514016, "end_time": "2026-03-24T02:04:29.183974+00:00", "exception": false, "start_time": "2026-03-24T02:03:53.669958+00:00", "status": "completed" } }, "outputs": [], "source": [ "avg_ret = (\n", " eval_df.group_by(\"timestamp\")\n", " .agg(pl.col(label_col).mean().alias(\"avg_fwd_ret\"))\n", " .sort(\"timestamp\")\n", ")\n", "\n", "date_features = (\n", " val_temporal.select([\"timestamp\"] + date_level_features)\n", " .unique(subset=[\"timestamp\"])\n", " .sort(\"timestamp\")\n", ")\n", "eval_ts = date_features.join(avg_ret, on=\"timestamp\", how=\"inner\").drop_nulls()\n", "\n", "for feat in date_level_features:\n", " x = eval_ts[feat].to_numpy()\n", " y = eval_ts[\"avg_fwd_ret\"].to_numpy()\n", " valid = ~(np.isnan(x) | np.isnan(y))\n", " if valid.sum() < 50:\n", " continue\n", " result = robust_ic(x[valid], y[valid], return_details=True)\n", " temporal_ic[feat] = result" ] }, { "cell_type": "code", "execution_count": 26, "id": "a68bb556", "metadata": { "execution": { "iopub.execute_input": "2026-03-24T02:04:29.188873Z", "iopub.status.busy": "2026-03-24T02:04:29.188784Z", "iopub.status.idle": "2026-03-24T02:04:29.191283Z", "shell.execute_reply": "2026-03-24T02:04:29.191004Z" }, "papermill": { "duration": 0.00524, "end_time": "2026-03-24T02:04:29.191502+00:00", "exception": false, "start_time": "2026-03-24T02:04:29.186262+00:00", "status": "completed" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Model-Based Feature Evaluation:\n", "shape: (14, 5)\n", "┌────────────────────┬───────────┬───────────┬─────────┬──────────────┐\n", "│ feature ┆ ic ┆ t_stat ┆ p_value ┆ bootstrap_se │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ f64 ┆ f64 ┆ f64 ┆ f64 │\n", "╞════════════════════╪═══════════╪═══════════╪═════════╪══════════════╡\n", "│ regime_prob_stress ┆ 0.169398 ┆ 2.807827 ┆ 0.0 ┆ 0.060331 │\n", "│ ffd_tlt ┆ 0.070389 ┆ 1.541184 ┆ 0.114 ┆ 0.045672 │\n", "│ ffd_lqd ┆ 0.066419 ┆ 1.324522 ┆ 0.139 ┆ 0.050146 │\n", "│ garch_cond_vol ┆ 0.051639 ┆ 9.646657 ┆ null ┆ 0.37886 │\n", "│ ffd_gld ┆ 0.044882 ┆ 0.797033 ┆ 0.386 ┆ 0.056311 │\n", "│ … ┆ … ┆ … ┆ … ┆ … │\n", "│ ffd_qqq ┆ -0.084412 ┆ -1.434737 ┆ 0.09 ┆ 0.058834 │\n", "│ ffd_efa ┆ -0.110375 ┆ -2.185295 ┆ 0.022 ┆ 0.050508 │\n", "│ ffd_spy ┆ -0.114751 ┆ -2.089767 ┆ 0.028 ┆ 0.054911 │\n", "│ ffd_vnq ┆ -0.143298 ┆ -3.032813 ┆ 0.0 ┆ 0.047249 │\n", "│ ffd_iwm ┆ -0.1577 ┆ -3.099353 ┆ 0.001 ┆ 0.050882 │\n", "└────────────────────┴───────────┴───────────┴─────────┴──────────────┘\n" ] } ], "source": [ "temporal_eval = pl.DataFrame(\n", " [\n", " {\n", " \"feature\": feat,\n", " \"ic\": stats[\"ic\"],\n", " \"t_stat\": stats.get(\"t_stat\", 0.0),\n", " \"p_value\": stats.get(\"p_value\"),\n", " \"bootstrap_se\": stats.get(\"bootstrap_std\", stats.get(\"bootstrap_se\", 0.0)),\n", " }\n", " for feat, stats in temporal_ic.items()\n", " ]\n", ").sort(\"ic\", descending=True)\n", "\n", "print(\"\\nModel-Based Feature Evaluation (validation periods only):\")\n", "print(temporal_eval)" ] }, { "cell_type": "markdown", "id": "dd0d5f47", "metadata": { "papermill": { "duration": 0.00201, "end_time": "2026-03-24T02:04:29.195576+00:00", "exception": false, "start_time": "2026-03-24T02:04:29.193566+00:00", "status": "completed" } }, "source": [ "**Interpretation**: GARCH conditional volatility is evaluated via\n", "cross-sectional IC (does higher predicted vol predict lower returns?).\n", "Date-level features (HMM, FFD) are evaluated via time-series IC against\n", "the average market return. Both types add value through different channels:\n", "GARCH enables cross-sectional differentiation, while regime features\n", "enable conditional strategies (e.g., reduce exposure in stress regimes).\n", "\n", "Because all features are now computed per fold with training-only\n", "parameters, the IC values reflect genuine out-of-sample signal strength\n", "rather than in-sample fit quality." ] }, { "cell_type": "markdown", "id": "833cb810", "metadata": { "papermill": { "duration": 0.00215, "end_time": "2026-03-24T02:04:29.199725+00:00", "exception": false, "start_time": "2026-03-24T02:04:29.197575+00:00", "status": "completed" } }, "source": [ "## Key Takeaways\n", "\n", "1. **Per-fold fitting eliminates parameter look-ahead**: HMM and GARCH\n", " parameters are estimated on each fold's training window only.\n", " The `fold` column in the output lets downstream models use the\n", " correct features for each CV fold.\n", "2. **HMM on aggregate market**: 2-state Gaussian HMM on SPY captures\n", " market-wide risk-on/risk-off regimes. All ETFs inherit the same state.\n", "3. **Filtered probabilities**: Forward algorithm only -- smoothed probs\n", " use future data and would introduce observation-level look-ahead bias.\n", "4. **Label switching prevention**: States sorted by variance ensures\n", " State 0 = calm, State 1 = stressed across all folds.\n", "5. **GARCH fit-then-filter**: `model.fix(params)` applies frozen training\n", " parameters to the full fold window, producing causal conditional\n", " volatility without re-estimation.\n", "6. **Fractional differencing**: Fixed $d$ by asset class (equities=0.4,\n", " fixed income=0.5) requires no fitting, but is computed per fold\n", " window for clean warmup handling.\n", "7. **Per-ETF GARCH**: Conditional volatility provides asset-specific risk\n", " dynamics, enabling cross-sectional differentiation (unlike date-level\n", " features which are shared across all ETFs).\n", "8. **Holdout fold**: A dedicated holdout fold (fit on all pre-holdout\n", " data) provides features for the final out-of-sample evaluation.\n", "\n", "**Next**: Ch11 models will join financial features (Ch8) with\n", "model-based features (Ch9) and use the `fold` column to ensure each\n", "CV fold uses only training-fitted temporal features." ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "tags,-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.3" }, "papermill": { "default_parameters": {}, "duration": 44.951602, "end_time": "2026-03-24T02:04:30.229699+00:00", "environment_variables": {}, "exception": null, "input_path": "case_studies/etfs/04_model_based_features.ipynb", "output_path": "/tmp/ml4t_3987651_04_model_based_features.ipynb", "parameters": { "GARCH_MIN_OBS": 50, "HMM_N_RESTARTS": 1, "MAX_SAMPLE_ASSETS": 3, "MAX_SYMBOLS": 5, "MIN_OBS": 10, "MIN_TRAIN_BARS": 10, "N_RESTARTS": 1, "START_DATE": "2024-06-01" }, "start_time": "2026-03-24T02:03:45.278097+00:00", "version": "2.7.0" } }, "nbformat": 4, "nbformat_minor": 5 }