# ---
# 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]
# # Spectral and Wavelet Features
#
# **Docker image**: `ml4t`
#
# This notebook demonstrates frequency-domain feature engineering:
# wavelet decomposition for multi-resolution analysis, rolling FFT for
# production spectral features, and Welch's method for robust power
# spectral density estimation.
#
# **Learning Objectives**:
# - Decompose returns into frequency components via wavelets (research tool)
# - Compute rolling FFT features: spectral energy, dominant frequency, entropy
# - Estimate power spectral density with Welch's method
# - Build a rolling spectral heatmap for time-frequency visualization
# - Compute `low_freq_ratio` as a regime-sensitive feature
#
# **Book Reference**: Chapter 9, Section 9.2 (Signal Transforms and Spectral Features)
#
# **Prerequisites**: `01_visual_diagnostics` for stationarity concepts.

# %%
"""Spectral and Wavelet Features — frequency-domain feature engineering."""

import importlib
import warnings
from typing import Any, cast

warnings.filterwarnings("ignore")

from datetime import datetime

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
from IPython.display import display
from ml4t.engineer.features.ml import fourier_features, rolling_entropy
from scipy.signal import welch

from data import load_etfs

pywt = cast(Any, importlib.import_module("pywt"))

# %% tags=["parameters"]
# Production defaults — Papermill injects overrides for CI
START_DATE = "2015-01-01"
END_DATE = "2024-12-31"

# %% [markdown]
# ## Load Data
#
# SPY returns provide a stationary signal suitable for spectral analysis.

# %%
etfs = load_etfs(symbols=["SPY"])
spy = etfs.select(["timestamp", "close"]).sort("timestamp")

START = datetime.strptime(START_DATE, "%Y-%m-%d")
END = datetime.strptime(END_DATE, "%Y-%m-%d")
spy = spy.filter((pl.col("timestamp") >= START) & (pl.col("timestamp") <= END))

spy = spy.with_columns(returns=pl.col("close").pct_change()).drop_nulls()

spy_pd = spy.to_pandas().set_index("timestamp")
returns = spy_pd["returns"].to_numpy(dtype=float, copy=True)
prices = spy_pd["close"].to_numpy(dtype=float, copy=True)
dates = spy_pd.index

print(f"SPY: {len(spy_pd):,} obs ({dates.min().date()} to {dates.max().date()})")

# %% [markdown]
# ---
#
# # Part 1: Wavelet Analysis (Offline Research Tool)
#
# Wavelets provide **multi-resolution analysis** — they reveal what time horizons
# carry signal vs noise. This is an **analysis tool**, not a production feature
# generator: standard DWT is non-causal (uses the entire series).
#
# **Key insight**: Use wavelet analysis to understand scale relationships, then
# translate those insights into simpler **causal proxies** (rolling averages,
# filters).

# %% [markdown]
# ### Wavelet Families for Financial Data
#
# Different wavelet families offer different trade-offs:
#
# | Family | Key Property | Best For |
# |--------|-------------|----------|
# | `db6` (Daubechies 6) | Good time-frequency localization | General purpose |
# | `sym6` (Symlet 6) | Near-symmetric, reduces phase distortion | Trend analysis |
# | `coif3` (Coiflet 3) | Near-zero moments, good reconstruction | Feature design |

# %%
families = pywt.families()
print("Available wavelet families:", families)
print()
for name, wname in [("Daubechies-6", "db6"), ("Symlet-6", "sym6"), ("Coiflet-3", "coif3")]:
    w = pywt.Wavelet(wname)
    print(f"  {name:15s} ({wname}): filter length={w.dec_len}, vanishing moments={w.dec_len // 2}")

# %% [markdown]
# ### Multi-Resolution Decomposition
#
# Decompose returns into approximation (trend) and detail (fluctuation)
# components at increasing time scales.


# %%
def wavelet_decompose(
    signal: np.ndarray,
    wavelet: str = "db6",
    level: int = 5,
) -> dict[str, np.ndarray]:
    """Decompose signal into approximation and detail components.

    WARNING: Non-causal — uses entire series. For analysis only.
    """
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    result = {}

    # Approximation (lowest frequency / trend)
    approx = pywt.upcoef("a", coeffs[0], wavelet, level=level, take=len(signal))
    result[f"A{level}"] = approx

    # Details (higher frequencies)
    for i, detail in enumerate(coeffs[1:], 1):
        d = pywt.upcoef("d", detail, wavelet, level=level - i + 1, take=len(signal))
        result[f"D{level - i + 1}"] = d

    return result


# %%
decomp = wavelet_decompose(returns, wavelet="db6", level=5)

fig, axes = plt.subplots(len(decomp) + 1, 1, figsize=(14, 12), sharex=True)

axes[0].plot(dates, returns, linewidth=0.3, alpha=0.7)
axes[0].set_title("Original Returns")
axes[0].set_ylabel("Return")

scale_labels = {
    "D1": "D1 (~2-4 days)",
    "D2": "D2 (~4-8 days)",
    "D3": "D3 (~8-16 days)",
    "D4": "D4 (~16-32 days)",
    "D5": "D5 (~32-64 days)",
    "A5": "A5 (>64 days, trend)",
}

for i, name in enumerate(sorted(decomp.keys()), 1):
    axes[i].plot(dates, decomp[name], linewidth=0.5)
    label = scale_labels.get(name, name)
    axes[i].set_title(label)
    axes[i].set_ylabel("Amplitude")

plt.suptitle("Wavelet Decomposition of SPY Returns (db6, 5 levels)", fontsize=14, fontweight="bold")
plt.show()

# %% [markdown]
# ### Scale Variance Analysis
#
# Which time scales contain the most variance (and potentially signal)?

# %%
total_var = np.var(returns)

scale_names = []
scale_pcts = []
for name in sorted(decomp.keys()):
    var = np.var(decomp[name])
    pct = 100 * var / total_var
    scale_names.append(scale_labels.get(name, name))
    scale_pcts.append(pct)

fig, ax = plt.subplots(figsize=(10, 4))
ax.barh(scale_names, scale_pcts)
ax.set_xlabel("Variance Explained (%)")
ax.set_title("Variance Decomposition by Wavelet Scale")
ax.invert_yaxis()
plt.show()

# %% [markdown]
# ### Wavelet Family Comparison
#
# Compare how db6, sym6, and coif3 decompose the same signal.

# %%
fig, axes = plt.subplots(3, 2, figsize=(14, 9), sharex=True)

for row, (wname, label) in enumerate(
    [("db6", "Daubechies-6"), ("sym6", "Symlet-6"), ("coif3", "Coiflet-3")]
):
    dec = wavelet_decompose(returns, wavelet=wname, level=5)

    # Show A5 (trend) and D1 (noise)
    axes[row, 0].plot(dates, dec["A5"], linewidth=0.8)
    axes[row, 0].set_title(f"{label}: A5 (Trend)")

    axes[row, 1].plot(dates, dec["D1"], linewidth=0.3, alpha=0.7)
    axes[row, 1].set_title(f"{label}: D1 (High-Frequency)")

plt.suptitle("Wavelet Family Comparison", fontsize=14, fontweight="bold")
plt.show()

# %% [markdown]
# **Finding**: The three wavelet families produce nearly identical trend components
# (A5) but differ in high-frequency detail (D1). Daubechies-6 (`db6`) is the
# recommended default — it provides the best balance of time-frequency localization
# for financial data. The choice of wavelet family matters far less than the choice
# of decomposition level.

# %% [markdown]
# ### Scale Insights → Causal Proxies
#
# Production features must be causal. Translate wavelet scale insights into
# rolling window features:
#
# | Wavelet Scale | Time Horizon | Causal Proxy |
# |---------------|-------------|--------------|
# | D1 | 2-4 days | 3-day rolling features |
# | D2 | 4-8 days | 5-day rolling features |
# | D3 | 8-16 days | 10-day rolling features |
# | D4 | 16-32 days | 21-day rolling features |
# | D5 | 32-64 days | 63-day rolling features |
# | A5 | >64 days | 126-day rolling features |

# %%
scale_to_window = {"D1": 3, "D2": 5, "D3": 10, "D4": 21, "D5": 63, "A5": 126}

proxy_results = {}
for scale, window in scale_to_window.items():
    # Rolling volatility as causal proxy for scale energy
    vol = np.full(len(returns), np.nan)
    for t in range(window, len(returns)):
        vol[t] = np.std(returns[t - window : t])
    proxy_results[f"vol_{window}d"] = vol

    # Compare with wavelet scale energy
    wavelet_energy = np.abs(decomp[scale])
    # Rolling mean of wavelet energy
    wav_roll = np.full(len(returns), np.nan)
    for t in range(window, len(returns)):
        wav_roll[t] = np.mean(wavelet_energy[t - window : t])
    proxy_results[f"wav_energy_{scale}"] = wav_roll

# Show correlation between causal proxy and wavelet energy
proxy_corr_rows = []
for scale, window in scale_to_window.items():
    vol = proxy_results[f"vol_{window}d"]
    wav = proxy_results[f"wav_energy_{scale}"]
    mask = ~np.isnan(vol) & ~np.isnan(wav)
    if mask.sum() > 50:
        corr = np.corrcoef(vol[mask], wav[mask])[0, 1]
        proxy_corr_rows.append({"scale": scale, "rolling_window_days": window, "correlation": corr})

display(pd.DataFrame(proxy_corr_rows))

# %% [markdown]
# ---
#
# # Part 2: Rolling FFT Features (Production)
#
# Fourier analysis decomposes time series into frequency components, revealing
# periodic structure invisible in the time domain. **Rolling FFT** provides
# production-ready causal features.

# %% [markdown]
# ### Rolling FFT Implementation
#
# For each rolling window, compute the FFT and extract features from the
# power spectrum.


# %%
def rolling_fft_features(
    signal: np.ndarray,
    window: int = 63,
    target_periods: list[int] | None = None,
) -> dict[str, np.ndarray]:
    """Compute spectral features using rolling FFT.

    Features extracted:
    - spectral_energy: Total power (excluding DC)
    - dominant_period: Period of peak frequency (in samples)
    - spectral_entropy: Entropy of normalized power spectrum
    - low_freq_ratio: Fraction of energy below 1/21 cycles/sample
    - energy at target periods: Power near specific calendar frequencies
    """
    if target_periods is None:
        target_periods = [5, 21, 63]

    n = len(signal)
    spectral_energy = np.full(n, np.nan)
    dominant_period = np.full(n, np.nan)
    spectral_entropy = np.full(n, np.nan)
    low_freq_ratio = np.full(n, np.nan)
    valid = np.zeros(n, dtype=bool)
    freq_energies = {p: np.full(n, np.nan) for p in target_periods}

    for t in range(window, n):
        segment = signal[t - window : t]

        # Detrend (remove mean)
        segment = segment - segment.mean()

        # FFT
        fft_vals = np.fft.rfft(segment)
        power = np.abs(fft_vals) ** 2
        freqs = np.fft.rfftfreq(window)

        # Total spectral energy (excluding DC)
        total_power = np.sum(power[1:])
        spectral_energy[t] = total_power

        # Dominant period
        if len(power) > 1 and total_power > 0:
            dom_idx = np.argmax(power[1:]) + 1
            if freqs[dom_idx] > 0:
                dominant_period[t] = 1.0 / freqs[dom_idx]

        # Spectral entropy
        if total_power > 0:
            p_norm = power[1:] / total_power
            p_norm = p_norm[p_norm > 0]
            spectral_entropy[t] = -np.sum(p_norm * np.log(p_norm))

        # Low-frequency ratio: energy below 1/21 cycles/sample (monthly+)
        low_mask = (freqs > 0) & (freqs < 1.0 / 21)
        if low_mask.any() and total_power > 0:
            low_freq_ratio[t] = np.sum(power[low_mask]) / total_power

        # Energy at target periods
        for tp in target_periods:
            target_freq = 1.0 / tp
            freq_idx = np.argmin(np.abs(freqs - target_freq))
            low_idx = max(1, freq_idx - 1)
            high_idx = min(len(power), freq_idx + 2)
            freq_energies[tp][t] = np.sum(power[low_idx:high_idx])

        valid[t] = True

    result = {
        "spectral_energy": spectral_energy,
        "dominant_period": dominant_period,
        "spectral_entropy": spectral_entropy,
        "low_freq_ratio": low_freq_ratio,
        "valid": valid,
    }
    for tp, energy in freq_energies.items():
        result[f"energy_period_{tp}"] = energy

    return result


# %%
fft_out = rolling_fft_features(returns, window=63, target_periods=[5, 21, 63])

valid_mask = fft_out["valid"]
valid_dates = dates[valid_mask]

fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)

ax = axes[0]
ax.plot(dates, returns, linewidth=0.3, alpha=0.5, color="gray")
energy_scaled = fft_out["spectral_energy"][valid_mask]
energy_scaled = energy_scaled / energy_scaled.max() * np.abs(returns).max()
ax.plot(valid_dates, energy_scaled, linewidth=1, color="blue", label="Spectral Energy (scaled)")
ax.set_title("Returns vs Spectral Energy")
ax.legend()

ax = axes[1]
ax.plot(valid_dates, fft_out["dominant_period"][valid_mask], linewidth=0.8, color="green")
ax.axhline(5, color="red", linestyle="--", linewidth=0.5, label="Weekly (5d)")
ax.axhline(21, color="orange", linestyle="--", linewidth=0.5, label="Monthly (21d)")
ax.set_title("Dominant Period (Days)")
ax.set_ylabel("Period")
ax.legend()

ax = axes[2]
ax.plot(valid_dates, fft_out["spectral_entropy"][valid_mask], linewidth=0.8, color="purple")
ax.set_title("Spectral Entropy (Complexity)")
ax.set_ylabel("Entropy")

ax = axes[3]
ax.plot(valid_dates, fft_out["low_freq_ratio"][valid_mask], linewidth=0.8, color="teal")
ax.axhline(0.5, color="red", linestyle="--", linewidth=0.5, label="50% threshold")
ax.set_title("Low-Frequency Ratio (Regime Indicator)")
ax.set_ylabel("Ratio")
ax.legend()

plt.suptitle("Rolling FFT Spectral Features (63-Day Window)", fontsize=14, fontweight="bold")
plt.show()

# %% [markdown]
# ### Calendar-Cycle Detection
#
# Compare energy at weekly, monthly, and quarterly frequencies to detect
# periodic patterns.

# %%
fig, ax = plt.subplots(figsize=(14, 4))

for tp, color, label in [(5, "blue", "Weekly"), (21, "green", "Monthly"), (63, "red", "Quarterly")]:
    energy = fft_out[f"energy_period_{tp}"][valid_mask]
    energy_norm = energy / np.nanmax(energy)
    ax.plot(valid_dates, energy_norm, linewidth=0.8, color=color, label=label, alpha=0.8)

ax.set_title("Normalized Energy at Calendar-Cycle Frequencies")
ax.set_ylabel("Normalized Energy")
ax.legend()
plt.show()

# %% [markdown]
# ---
#
# # Part 3: Welch's Method (Robust PSD Estimation)
#
# Welch's method averages periodograms over overlapping segments, reducing
# variance at the cost of frequency resolution. It is more robust than a
# single FFT for estimating the power spectral density (PSD).

# %% [markdown]
# ### Full-Sample PSD
#
# Compute the PSD of the entire return series to identify persistent
# periodicities.

# %%
# Welch PSD of full return series
freqs_welch, psd_welch = welch(returns, fs=1.0, nperseg=128, noverlap=64)

# Convert frequency to period for interpretability
periods_welch = np.where(freqs_welch > 0, 1.0 / freqs_welch, np.inf)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

ax = axes[0]
ax.semilogy(freqs_welch[1:], psd_welch[1:], linewidth=1)
for freq, label, color in [
    (1 / 5, "Weekly", "red"),
    (1 / 21, "Monthly", "orange"),
    (1 / 63, "Quarterly", "green"),
]:
    ax.axvline(freq, color=color, linestyle="--", linewidth=0.8, label=label)
ax.set_xlabel("Frequency (cycles/day)")
ax.set_ylabel("PSD (log scale)")
ax.set_title("Power Spectral Density (Welch)")
ax.legend()

ax = axes[1]
ax.plot(periods_welch[1:50], psd_welch[1:50], linewidth=1, marker="o", markersize=3)
for period, label in [(5, "Weekly"), (21, "Monthly")]:
    ax.axvline(period, color="red", linestyle="--", linewidth=0.5, label=label)
ax.set_xlabel("Period (trading days)")
ax.set_ylabel("PSD")
ax.set_title("PSD vs Period (Low-Frequency Detail)")
ax.legend()

plt.suptitle("Welch PSD Estimation — SPY Returns", fontsize=14, fontweight="bold")
plt.show()

# %% [markdown]
# ### Rolling Welch PSD
#
# Compute Welch PSD in a rolling window to track how the spectral profile
# changes over time.


# %%
def rolling_welch_psd(
    signal: np.ndarray,
    window: int = 252,
    step: int = 5,
    nperseg: int = 64,
) -> tuple[np.ndarray, Any, np.ndarray]:
    """Compute rolling Welch PSD.

    Returns (time_indices, frequencies, psd_matrix) where psd_matrix
    has shape (n_windows, n_freqs).
    """
    n = len(signal)
    time_indices = []
    psd_list = []
    freqs_out = None

    for end in range(window, n, step):
        segment = signal[end - window : end]
        f, pxx = welch(segment, fs=1.0, nperseg=min(nperseg, window // 2), noverlap=nperseg // 2)
        if freqs_out is None:
            freqs_out = f
        psd_list.append(pxx)
        time_indices.append(end)

    return np.array(time_indices), freqs_out, np.array(psd_list)


# %%
time_idx, freqs_w, psd_matrix = rolling_welch_psd(returns, window=252, step=5, nperseg=64)

# %% [markdown]
# ### Time-Frequency Heatmap
#
# A spectrogram-like visualization showing how spectral energy shifts over
# time. Bright regions indicate periods with concentrated energy at a
# particular frequency.

# %%
# Convert to periods for y-axis, cap at 126 days
period_mask = freqs_w > 1.0 / 126  # Only show periods up to 126 days
freqs_plot = freqs_w[period_mask]
periods_plot = 1.0 / freqs_plot
psd_plot = psd_matrix[:, period_mask]

fig, ax = plt.subplots(figsize=(14, 6))

# Log-scale PSD for better visibility
log_psd = np.log10(psd_plot.T + 1e-15)

im = ax.pcolormesh(
    dates[time_idx],
    periods_plot,
    log_psd,
    cmap="YlOrRd",
    shading="auto",
)

ax.set_ylabel("Period (trading days)")
ax.set_title("Rolling Spectral Heatmap (Welch PSD)")
ax.set_ylim(2, 126)

# Mark calendar frequencies
for period, label in [(5, "Weekly"), (21, "Monthly"), (63, "Quarterly")]:
    ax.axhline(period, color="cyan", linestyle="--", linewidth=0.8, alpha=0.7)
    ax.text(dates[time_idx[-1]], period * 1.05, label, fontsize=8, color="cyan", va="bottom")

plt.colorbar(im, ax=ax, label="log10(PSD)")
fig.subplots_adjust(left=0.08, right=0.92, top=0.92, bottom=0.12)
plt.show()

# %% [markdown]
# ---
#
# # Part 4: Consolidated Feature Extraction
#
# Combine spectral features from FFT and Welch into a single feature table
# suitable for downstream ML.

# %% [markdown]
# ### Multi-Window Spectral Features
#
# Use multiple FFT window sizes to capture different frequency resolutions.

# %%
windows = [21, 63, 126]
all_features = {"timestamp": dates}

for w in windows:
    fft = rolling_fft_features(returns, window=w, target_periods=[5, 21])
    all_features[f"spectral_energy_{w}d"] = fft["spectral_energy"]
    all_features[f"spectral_entropy_{w}d"] = fft["spectral_entropy"]
    all_features[f"low_freq_ratio_{w}d"] = fft["low_freq_ratio"]
    all_features[f"dominant_period_{w}d"] = fft["dominant_period"]

feature_df = pl.DataFrame(all_features).drop_nulls()

print(f"Feature table: {len(feature_df):,} rows x {len(feature_df.columns)} columns")
feature_df.tail(10)

# %% [markdown]
# ### Predictive Value of Spectral Features
#
# Compute information coefficient (rank correlation with forward returns).

# %%
# Add forward returns
fwd_5d = np.full(len(returns), np.nan)
fwd_5d[:-5] = (prices[5:] - prices[:-5]) / prices[:-5]

fwd_df = pl.DataFrame({"timestamp": dates, "fwd_ret_5d": fwd_5d})
eval_df = feature_df.join(fwd_df, on="timestamp", how="inner").drop_nulls()

feature_cols = [c for c in eval_df.columns if c not in ("timestamp", "fwd_ret_5d")]
fwd = eval_df["fwd_ret_5d"].to_numpy()

ic_rows = []
for col in feature_cols:
    vals = eval_df[col].to_numpy()
    mask = np.isfinite(vals) & np.isfinite(fwd)
    if mask.sum() > 100:
        ic = np.corrcoef(vals[mask], fwd[mask])[0, 1]
        ic_rows.append({"feature": col, "ic_5d_fwd": ic, "n_obs": int(mask.sum())})

ic_df = pd.DataFrame(ic_rows).sort_values("ic_5d_fwd", key=lambda s: s.abs(), ascending=False)
display(ic_df)

# %% [markdown]
# ## Feature Catalog: Spectral Features
#
# | Feature | Source | Computation | Update |
# |---------|--------|-------------|--------|
# | `spectral_energy` | Rolling FFT | Total power (excl. DC) | Daily |
# | `dominant_period` | Rolling FFT | Period of peak frequency | Daily |
# | `spectral_entropy` | Rolling FFT | Normalized spectrum entropy | Daily |
# | `low_freq_ratio` | Rolling FFT | Energy below 1/21 cycles/sample | Daily |
# | `energy_period_5` | Rolling FFT | Power near weekly frequency | Daily |
# | `energy_period_21` | Rolling FFT | Power near monthly frequency | Daily |
#
# **Wavelet features** are for research only (non-causal). Translate scale
# insights into rolling-window causal proxies.

# %% [markdown]
# ## ml4t-engineer: Spectral and Entropy Features
#
# The manual rolling FFT above requires ~30 lines of NumPy code per feature.
# `ml4t-engineer` provides `fourier_features()` for spectral decomposition
# and `rolling_entropy()` for information-theoretic features — both as
# Polars expressions.

# %%
# Load SPY for ml4t-engineer demonstration
spy_ml4t = (
    etfs.filter(pl.col("symbol") == "SPY")
    .select(["timestamp", "close"])
    .sort("timestamp")
    .with_columns(returns=pl.col("close").pct_change())
    .drop_nulls()
)

# Fourier features capture periodic patterns
fourier_exprs = fourier_features("close", n_components=5)
spy_spectral = spy_ml4t.with_columns(**fourier_exprs)

# Rolling entropy measures predictability
spy_spectral = spy_spectral.with_columns(
    entropy=rolling_entropy("returns", window=50, n_bins=10),
)

print("=== ml4t-engineer: Spectral Features ===")
print(f"Fourier components: {len(fourier_exprs)} (sin/cos pairs for 5 harmonics)")
entropy_vals = spy_spectral["entropy"].drop_nulls()
print(f"Entropy: mean={entropy_vals.mean():.4f}, std={entropy_vals.std():.4f}")
print("  High entropy = unpredictable (efficient market)")
print("  Low entropy = patterned (potential signal)")

# %% [markdown]
# `fourier_features()` generates sin/cos pairs at different frequencies,
# capturing periodic structure in the price series. `rolling_entropy()`
# quantifies the predictability of returns using Shannon entropy — a
# complement to spectral energy that works in the information domain.

# %% [markdown]
# ## Key Takeaways
#
# 1. **Wavelets are analysis tools**: Multi-resolution decomposition reveals
#    which time scales carry signal, but DWT is non-causal — translate
#    insights into rolling window proxies for production
# 2. **Rolling FFT provides production features**: `spectral_energy`,
#    `dominant_period`, `spectral_entropy`, and `low_freq_ratio` are all
#    causal and can update daily
# 3. **Welch's method is more robust**: Averaging over overlapping segments
#    reduces variance in PSD estimates; use for full-sample analysis and
#    time-frequency heatmaps
# 4. **`low_freq_ratio` is regime-sensitive**: High values indicate
#    trend-dominated regimes; low values indicate noisy, mean-reverting
#    markets — useful as a conditioning feature
# 5. **Multiple window sizes capture different scales**: 21-day (short-term),
#    63-day (medium), 126-day (long-term) windows provide complementary
#    frequency resolution
#
# **Next**: See `06_path_signatures` for path-based features and
# `07_arima_features` for time-domain forecasting features.