//! Regression-based forecasting models.
//!
//! Bridges external regression estimators (e.g., OLS from `anofox-regression`)
//! into the [`Forecaster`](crate::models::Forecaster) trait, enabling them to participate in pipelines,
//! model registries, ensembles, and cross-validation.
//!
//! # Feature engineering
//!
//! Time-series forecasting with regression requires features. The
//! [`RegressionFeatures`] builder configures which features are constructed
//! from a [`TimeSeries`](crate::core::TimeSeries) before fitting:
//!
//! | Feature          | Description |
//! |-----------------|-------------|
//! | Trend index      | Linear index `0, 1, …, n-1` |
//! | Lags             | `y[t-1], y[t-2], …, y[t-max_lag]` |
//! | Exogenous regressors | From `TimeSeries::all_regressors()` |
//!
//! # Example
//!
//! ```rust,ignore
//! use anofox_forecast::models::regression::{RegressionForecaster, RegressionFeatures};
//!
//! // OLS with trend + 3 lags + exogenous regressors
//! let mut model = RegressionForecaster::ols(
//!     RegressionFeatures::new().trend().lags(3),
//! );
//! model.fit(&ts)?;
//! let forecast = model.predict(12)?;
//! ```

#[cfg(feature = "postprocess")]
mod ols_impl {
    use std::collections::HashMap;
    use std::sync::Arc;

    use anofox_regression::core::IntervalType;
    use anofox_regression::solvers::{
        BlsRegressor, ElasticNetRegressor, FittedRegressor, InformationCriterion,
        LmDynamicRegressor, OlsRegressor, PoissonRegressor, QuantileRegressor, Regressor,
        RidgeRegressor, RlsRegressor, TweedieRegressor, WlsRegressor,
    };
    use faer::{Col, Mat};

    use crate::core::{Forecast, TimeSeries};
    use crate::error::{ForecastError, Result};
    use crate::models::{validate_series_complete, Forecaster};
    use crate::seasonality::dummy::DummySeasonality;
    use crate::seasonality::exponential_trend::ExponentialTrend;
    use crate::seasonality::fourier::fourier_terms;
    use crate::seasonality::polynomial::PolynomialTrend;
    use crate::seasonality::theilsen::TheilSenTrend;
    use crate::seasonality::traits::{Recency, SeasonalComponent, TrendComponent};

    // ── Feature safety classification ────────────────────────────────

    /// Classification of a feature by its data-leakage risk in cross-validation.
    ///
    /// | Level | Examples | CV requirement |
    /// |---|---|---|
    /// | `Deterministic` | Fourier terms, raw trend index, lags | None — always safe |
    /// | `DataDependent` | Fitted PolynomialTrend, DummySeasonal | Re-fit per fold |
    /// | `Structural` | Changepoint regime indicator | Re-detect per fold; flag if break in test |
    /// | `External` | User-provided exogenous regressors | User's responsibility |
    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
    pub enum FeatureSafety {
        /// Computed from the time index only — no data leakage possible.
        Deterministic,
        /// Fitted from training data — must be re-fit per CV fold.
        DataDependent,
        /// Derived from structural analysis — forward-filled during prediction.
        Structural,
        /// User-provided exogenous regressors — user's responsibility.
        External,
    }

    // ── Structural feature trait ─────────────────────────────────────

    /// A feature derived from structural analysis that is forward-filled during prediction.
    ///
    /// Implementations compute one or more columns from training data. During
    /// prediction, each column is frozen at a constant value — the model
    /// never sees future structural changes.
    ///
    /// # Fill strategies
    ///
    /// The [`fill_values`](Self::fill_values) method returns the constant value
    /// to repeat for each column during prediction. Different features use
    /// different strategies:
    ///
    /// - **Forward-fill** (changepoints): return the value at the last training
    ///   index — the model continues in the last known regime.
    /// - **Constant fill** (outlier indicators): return a fixed default (e.g., 0.0)
    ///   — the model assumes no outliers in the forecast period.
    pub trait StructuralFeature: std::fmt::Debug + Send + Sync {
        /// Column names this feature produces.
        fn column_names(&self) -> Vec<String>;

        /// Number of output columns.
        fn n_columns(&self) -> usize {
            self.column_names().len()
        }

        /// Compute column values for observation indices `0..n`.
        ///
        /// Returns one `Vec<f64>` of length `n` per column.
        fn compute(&self, n: usize) -> Vec<Vec<f64>>;

        /// Constant values to use for each column during prediction.
        ///
        /// Returns one value per column. Each value is repeated for every
        /// forecast step.
        fn fill_values(&self, n_train: usize) -> Vec<f64>;

        /// Human-readable name for reporting.
        fn name(&self) -> &str;
    }

    // ── Changepoint encoding ─────────────────────────────────────────

    /// How changepoint locations are encoded as regression features.
    ///
    /// With `k` detected changepoints (e.g., at indices 50 and 120):
    ///
    /// | Encoding | Columns | Description |
    /// |---|---|---|
    /// | `StepFunctions` | `k` binary columns | Each column flips 0→1 at its CP |
    /// | `RegimeIndex` | 1 column | Values 0, 1, …, k per segment |
    /// | `CumulativeCount` | 1 column | Count of CPs at or before index |
    #[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
    pub enum ChangepointEncoding {
        /// One binary column per changepoint (0 before, 1 at-or-after).
        #[default]
        StepFunctions,
        /// Single column with ordinal regime index (0, 1, …, k).
        RegimeIndex,
        /// Single column counting changepoints at or before each index.
        CumulativeCount,
    }

    /// Changepoint-based structural feature with configurable encoding.
    ///
    /// Produces regression columns that encode the detected changepoint
    /// locations. During prediction, columns are forward-filled with
    /// their last training values (the model stays in the last known regime).
    #[derive(Debug, Clone)]
    pub struct ChangepointFeature {
        /// Sorted changepoint indices.
        indices: Vec<usize>,
        /// How to encode the changepoints.
        encoding: ChangepointEncoding,
    }

    impl ChangepointFeature {
        /// Create a new changepoint feature with the given encoding.
        pub fn new(mut indices: Vec<usize>, encoding: ChangepointEncoding) -> Self {
            indices.sort_unstable();
            Self { indices, encoding }
        }

        /// Create a changepoint feature with the default `StepFunctions` encoding.
        pub fn step_functions(indices: Vec<usize>) -> Self {
            Self::new(indices, ChangepointEncoding::StepFunctions)
        }
    }

    impl StructuralFeature for ChangepointFeature {
        fn column_names(&self) -> Vec<String> {
            match self.encoding {
                ChangepointEncoding::StepFunctions => self
                    .indices
                    .iter()
                    .enumerate()
                    .map(|(i, _)| format!("__cp_step_{}", i + 1))
                    .collect(),
                ChangepointEncoding::RegimeIndex => vec!["__cp_regime".into()],
                ChangepointEncoding::CumulativeCount => vec!["__cp_count".into()],
            }
        }

        fn compute(&self, n: usize) -> Vec<Vec<f64>> {
            match self.encoding {
                ChangepointEncoding::StepFunctions => self
                    .indices
                    .iter()
                    .map(|&cp| (0..n).map(|i| if i >= cp { 1.0 } else { 0.0 }).collect())
                    .collect(),
                ChangepointEncoding::RegimeIndex | ChangepointEncoding::CumulativeCount => {
                    let col: Vec<f64> = (0..n)
                        .map(|i| self.indices.iter().filter(|&&cp| cp <= i).count() as f64)
                        .collect();
                    vec![col]
                }
            }
        }

        fn fill_values(&self, n_train: usize) -> Vec<f64> {
            match self.encoding {
                ChangepointEncoding::StepFunctions => self
                    .indices
                    .iter()
                    .map(|&cp| {
                        if n_train > 0 && cp < n_train {
                            1.0
                        } else {
                            0.0
                        }
                    })
                    .collect(),
                ChangepointEncoding::RegimeIndex | ChangepointEncoding::CumulativeCount => {
                    let count = self.indices.iter().filter(|&&cp| cp < n_train).count() as f64;
                    vec![count]
                }
            }
        }

        fn name(&self) -> &str {
            "ChangepointFeature"
        }
    }

    // ── Recursive feature trait ──────────────────────────────────────

    /// A feature that depends on past observations and must be **recomputed
    /// at every horizon step** during recursive multi-step prediction.
    ///
    /// Unlike [`StructuralFeature`], which is forward-filled with a constant
    /// value, a `RecursiveFeature` consumes the rolling history buffer
    /// (training tail + predictions emitted so far) and returns a fresh
    /// value for the next prediction step. This is the right abstraction
    /// for rolling statistics, EWMs, and any other target-derived feature
    /// whose "correct" value at horizon `h` depends on forecasts at
    /// horizons `0..h`.
    ///
    /// # Contract
    ///
    /// - [`compute_fit`](Self::compute_fit) is called once per training row
    ///   with the full series and the target index. Implementations must
    ///   read only from positions strictly before `target_idx` (the feature
    ///   must not see its own target).
    /// - [`compute_predict`](Self::compute_predict) is called once per
    ///   horizon step with the rolling `recent` buffer. `recent.last()` is
    ///   the most recently known (or predicted) value **before** the
    ///   observation being forecast.
    /// - [`warmup`](Self::warmup) is the minimum history required before
    ///   the feature is valid. The enclosing [`RegressionFeatures`] uses
    ///   this to grow the lag offset and drops unusable warmup rows.
    pub trait RecursiveFeature: std::fmt::Debug + Send + Sync {
        /// Column names this feature produces.
        fn column_names(&self) -> Vec<String>;

        /// Number of output columns.
        fn n_columns(&self) -> usize {
            self.column_names().len()
        }

        /// Minimum history length required for the feature to be valid.
        ///
        /// Rows where `target_idx < warmup()` are unusable at fit time
        /// and must be dropped.
        fn warmup(&self) -> usize;

        /// Populate the row corresponding to predicting `values[target_idx]`.
        ///
        /// Implementations must read only from `values[..target_idx]`
        /// (strictly-before, no leakage). `out.len() == n_columns()`.
        fn compute_fit(&self, values: &[f64], target_idx: usize, out: &mut [f64]);

        /// Populate the next prediction row from the rolling history buffer.
        ///
        /// `recent` is the training tail extended with predictions emitted
        /// so far; `recent.last()` is the most recent known or predicted
        /// value before the observation being forecast. `out.len() == n_columns()`.
        fn compute_predict(&self, recent: &[f64], out: &mut [f64]);

        /// Populate the next prediction row with knowledge of the **absolute**
        /// horizon step index.
        ///
        /// Default impl delegates to [`compute_predict`](Self::compute_predict)
        /// — most features (e.g. [`RollingFeature`]) only need the rolling
        /// buffer. Features that depend on the absolute timestep (e.g.
        /// [`EventDistanceFeature`]) override this to use `step_h + n_train`.
        ///
        /// - `step_h`: 0-based offset within the forecast horizon.
        /// - `n_train`: total length of the training data (so the absolute
        ///   index of the row being predicted is `n_train + step_h`).
        fn compute_predict_at(
            &self,
            recent: &[f64],
            _step_h: usize,
            _n_train: usize,
            out: &mut [f64],
        ) {
            self.compute_predict(recent, out);
        }

        /// Human-readable name for reporting.
        fn name(&self) -> &str;

        /// Clone into a boxed trait object (for `Vec<Arc<dyn RecursiveFeature>>`).
        fn clone_box(&self) -> Box<dyn RecursiveFeature>;
    }

    // ── Rolling statistic kinds ──────────────────────────────────────

    /// Which rolling statistic to compute in a [`RollingFeature`].
    #[derive(Debug, Clone, Copy, PartialEq)]
    pub enum RollingStatKind {
        /// Arithmetic mean of the window.
        Mean,
        /// Sample standard deviation with Bessel's correction (ddof=1).
        Std,
        /// Sample variance with Bessel's correction (ddof=1).
        Var,
        /// Minimum of the window.
        Min,
        /// Maximum of the window.
        Max,
        /// Median of the window (linear interpolation of order statistics).
        Median,
        /// Sum of the window.
        Sum,
        /// Exponentially weighted mean with smoothing factor `alpha` ∈ (0, 1].
        ///
        /// Updated iteratively over the window as
        /// `s_k = alpha · x_k + (1 - alpha) · s_{k-1}`, seeded with the first
        /// value. The `window` parameter determines how many past observations
        /// contribute to the EWM — larger windows produce smoother values.
        EwmMean { alpha: f64 },
        /// Exponentially weighted standard deviation over the window.
        ///
        /// Equivalent to `EwmVar` followed by `sqrt`.
        EwmStd { alpha: f64 },
        /// Empirical τ-quantile (linear interpolation between order statistics).
        /// `tau ∈ [0, 1]`. NumPy-default behaviour.
        Quantile { tau: f64 },
        /// Range = max − min.
        Range,
        /// Interquartile range = Q₀.₇₅ − Q₀.₂₅.
        Iqr,
        /// Sample skewness (Fisher's, biased moment estimator).
        ///
        /// `(1/n) Σ ((x − μ) / σ)³` where σ is the population std (ddof = 0).
        /// Returns 0 when the window is constant or has fewer than 2 points.
        Skew,
        /// Sample excess kurtosis (Fisher's, biased moment estimator).
        ///
        /// `(1/n) Σ ((x − μ) / σ)⁴ − 3`. Returns 0 when the window is constant
        /// or has fewer than 2 points.
        Kurt,
        /// Slope of an OLS regression of the window values against
        /// `t = 0..n` (i.e. how fast the series is moving locally).
        Slope,
        /// Fractional rank of the **most recent** value in the window, in
        /// `[0, 1]`. 0 → smallest in window, 1 → largest. Useful for detecting
        /// regime extremes without depending on absolute scale.
        Rank,
        /// Z-score of the most recent value relative to the rest of the
        /// window: `(window.last() − mean(window)) / std(window)`. Returns 0
        /// when the window is constant.
        ZScore,
        /// Count of values strictly greater than `threshold`.
        CountAbove { threshold: f64 },
        /// Count of values strictly less than `threshold`.
        CountBelow { threshold: f64 },
    }

    impl RollingStatKind {
        fn short_name(&self) -> String {
            match self {
                Self::Mean => "mean".into(),
                Self::Std => "std".into(),
                Self::Var => "var".into(),
                Self::Min => "min".into(),
                Self::Max => "max".into(),
                Self::Median => "median".into(),
                Self::Sum => "sum".into(),
                Self::EwmMean { .. } => "ewm_mean".into(),
                Self::EwmStd { .. } => "ewm_std".into(),
                Self::Quantile { tau } => format!("q{}", tau),
                Self::Range => "range".into(),
                Self::Iqr => "iqr".into(),
                Self::Skew => "skew".into(),
                Self::Kurt => "kurt".into(),
                Self::Slope => "slope".into(),
                Self::Rank => "rank".into(),
                Self::ZScore => "zscore".into(),
                Self::CountAbove { threshold } => format!("countabove_{}", threshold),
                Self::CountBelow { threshold } => format!("countbelow_{}", threshold),
            }
        }

        fn compute(&self, window: &[f64]) -> f64 {
            let n = window.len();
            if n == 0 {
                return 0.0;
            }
            match *self {
                Self::Mean => window.iter().sum::<f64>() / n as f64,
                Self::Sum => window.iter().sum::<f64>(),
                Self::Min => window.iter().cloned().fold(f64::INFINITY, f64::min),
                Self::Max => window.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
                Self::Var => sample_variance_of(window),
                Self::Std => sample_variance_of(window).sqrt(),
                Self::Median => median_of(window),
                Self::EwmMean { alpha } => ewm_mean_of(window, alpha),
                Self::EwmStd { alpha } => ewm_var_of(window, alpha).sqrt(),
                Self::Quantile { tau } => quantile_of(window, tau),
                Self::Range => {
                    let mn = window.iter().cloned().fold(f64::INFINITY, f64::min);
                    let mx = window.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
                    mx - mn
                }
                Self::Iqr => quantile_of(window, 0.75) - quantile_of(window, 0.25),
                Self::Skew => skewness_of(window),
                Self::Kurt => excess_kurtosis_of(window),
                Self::Slope => ols_slope_of(window),
                Self::Rank => fractional_rank_of_last(window),
                Self::ZScore => z_score_of_last(window),
                Self::CountAbove { threshold } => {
                    window.iter().filter(|&&x| x > threshold).count() as f64
                }
                Self::CountBelow { threshold } => {
                    window.iter().filter(|&&x| x < threshold).count() as f64
                }
            }
        }
    }

    #[inline]
    fn sample_variance_of(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n < 2 {
            return 0.0;
        }
        let m = xs.iter().sum::<f64>() / n as f64;
        let sum_sq: f64 = xs.iter().map(|x| (x - m) * (x - m)).sum();
        sum_sq / (n - 1) as f64
    }

    #[inline]
    fn median_of(xs: &[f64]) -> f64 {
        let mut v: Vec<f64> = xs.to_vec();
        v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let n = v.len();
        if n == 0 {
            return 0.0;
        }
        if n % 2 == 1 {
            v[n / 2]
        } else {
            0.5 * (v[n / 2 - 1] + v[n / 2])
        }
    }

    #[inline]
    fn ewm_mean_of(xs: &[f64], alpha: f64) -> f64 {
        if xs.is_empty() {
            return 0.0;
        }
        let mut s = xs[0];
        for &x in &xs[1..] {
            s = alpha * x + (1.0 - alpha) * s;
        }
        s
    }

    #[inline]
    fn quantile_of(xs: &[f64], tau: f64) -> f64 {
        let n = xs.len();
        if n == 0 {
            return 0.0;
        }
        if n == 1 {
            return xs[0];
        }
        let mut v: Vec<f64> = xs.to_vec();
        v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let pos = tau.clamp(0.0, 1.0) * (n - 1) as f64;
        let lo = pos.floor() as usize;
        let hi = pos.ceil() as usize;
        if lo == hi {
            v[lo]
        } else {
            let frac = pos - lo as f64;
            v[lo] * (1.0 - frac) + v[hi] * frac
        }
    }

    #[inline]
    fn skewness_of(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n < 2 {
            return 0.0;
        }
        let nf = n as f64;
        let m = xs.iter().sum::<f64>() / nf;
        let mut m2 = 0.0;
        let mut m3 = 0.0;
        for &x in xs {
            let d = x - m;
            m2 += d * d;
            m3 += d * d * d;
        }
        m2 /= nf;
        m3 /= nf;
        if m2 <= 0.0 {
            return 0.0;
        }
        m3 / m2.powf(1.5)
    }

    #[inline]
    fn excess_kurtosis_of(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n < 2 {
            return 0.0;
        }
        let nf = n as f64;
        let m = xs.iter().sum::<f64>() / nf;
        let mut m2 = 0.0;
        let mut m4 = 0.0;
        for &x in xs {
            let d = x - m;
            let d2 = d * d;
            m2 += d2;
            m4 += d2 * d2;
        }
        m2 /= nf;
        m4 /= nf;
        if m2 <= 0.0 {
            return 0.0;
        }
        m4 / (m2 * m2) - 3.0
    }

    #[inline]
    fn ols_slope_of(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n < 2 {
            return 0.0;
        }
        let nf = n as f64;
        let t_mean = (nf - 1.0) * 0.5;
        let y_mean = xs.iter().sum::<f64>() / nf;
        let mut num = 0.0;
        let mut den = 0.0;
        for (i, &y) in xs.iter().enumerate() {
            let dt = i as f64 - t_mean;
            num += dt * (y - y_mean);
            den += dt * dt;
        }
        if den == 0.0 {
            0.0
        } else {
            num / den
        }
    }

    #[inline]
    fn fractional_rank_of_last(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n == 0 {
            return 0.0;
        }
        if n == 1 {
            return 0.5;
        }
        let last = xs[n - 1];
        let mut less = 0usize;
        let mut equal = 0usize;
        for &x in xs {
            if x < last {
                less += 1;
            } else if x == last {
                equal += 1;
            }
        }
        // Average rank when ties are present, normalised to [0, 1].
        ((less as f64) + 0.5 * (equal as f64 - 1.0)) / (n - 1) as f64
    }

    #[inline]
    fn z_score_of_last(xs: &[f64]) -> f64 {
        let n = xs.len();
        if n < 2 {
            return 0.0;
        }
        let m = xs.iter().sum::<f64>() / n as f64;
        let var = sample_variance_of(xs);
        if var <= 0.0 {
            return 0.0;
        }
        (xs[n - 1] - m) / var.sqrt()
    }

    #[inline]
    fn ewm_var_of(xs: &[f64], alpha: f64) -> f64 {
        if xs.len() < 2 {
            return 0.0;
        }
        // Pandas-style EWM variance with adjust=False, bias=False is complex;
        // we use the simpler biased recursion:
        // s_k = α x_k + (1-α) s_{k-1}
        // v_k = (1-α) (v_{k-1} + α (x_k - s_{k-1})²)
        let mut s = xs[0];
        let mut v = 0.0;
        for &x in &xs[1..] {
            let diff = x - s;
            v = (1.0 - alpha) * (v + alpha * diff * diff);
            s = alpha * x + (1.0 - alpha) * s;
        }
        v
    }

    // ── Rolling feature ──────────────────────────────────────────────

    /// A rolling-window statistic of the target used as a regression feature.
    ///
    /// Computes `kind(values[target - lag - window + 1 ..= target - lag])` for
    /// each training row, and is recomputed at every horizon step during
    /// multi-step prediction using the rolling history buffer.
    ///
    /// # Leakage guard
    ///
    /// `lag == 0` would include the target value in its own feature window,
    /// which is data leakage for any rolling statistic of the target series.
    /// Both [`RollingFeature::new`] and [`RollingFeature::with_lag`] reject
    /// `lag == 0` at construction time. The default `lag = 1` is always safe.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// use anofox_forecast::models::regression::{RollingFeature, RollingStatKind};
    ///
    /// // Rolling mean of the last 7 values (lag=1 by default)
    /// let f = RollingFeature::new(7, RollingStatKind::Mean).unwrap();
    ///
    /// // Rolling std of values from 2..=8 steps ago
    /// let f = RollingFeature::with_lag(7, 2, RollingStatKind::Std).unwrap();
    /// ```
    #[derive(Debug, Clone)]
    pub struct RollingFeature {
        window: usize,
        lag: usize,
        kind: RollingStatKind,
    }

    impl RollingFeature {
        /// Create a rolling feature with the given `window` and `kind`, using
        /// the default `lag = 1` (safe — window excludes the current target).
        pub fn new(window: usize, kind: RollingStatKind) -> Result<Self> {
            Self::with_lag(window, 1, kind)
        }

        /// Create a rolling feature with an explicit `lag`.
        ///
        /// `lag >= 1` is required to prevent target leakage.
        pub fn with_lag(window: usize, lag: usize, kind: RollingStatKind) -> Result<Self> {
            if window == 0 {
                return Err(ForecastError::InvalidParameter(
                    "RollingFeature: window must be >= 1".to_string(),
                ));
            }
            if lag == 0 {
                return Err(ForecastError::InvalidParameter(
                    "RollingFeature: lag must be >= 1 to avoid target leakage. \
                     The default lag = 1 (via RollingFeature::new) is the \
                     standard choice."
                        .to_string(),
                ));
            }
            if let RollingStatKind::EwmMean { alpha } | RollingStatKind::EwmStd { alpha } = kind {
                if !(0.0 < alpha && alpha <= 1.0) {
                    return Err(ForecastError::InvalidParameter(format!(
                        "RollingFeature: EWM alpha must satisfy 0 < α ≤ 1, got {}",
                        alpha
                    )));
                }
            }
            if let RollingStatKind::Quantile { tau } = kind {
                if !(0.0..=1.0).contains(&tau) || tau.is_nan() {
                    return Err(ForecastError::InvalidParameter(format!(
                        "RollingFeature: Quantile tau must satisfy 0 ≤ τ ≤ 1, got {}",
                        tau
                    )));
                }
            }
            Ok(Self { window, lag, kind })
        }

        /// Return the window size.
        pub fn window(&self) -> usize {
            self.window
        }

        /// Return the lag.
        pub fn lag(&self) -> usize {
            self.lag
        }

        /// Return the statistic kind.
        pub fn kind(&self) -> RollingStatKind {
            self.kind
        }
    }

    impl RecursiveFeature for RollingFeature {
        fn column_names(&self) -> Vec<String> {
            vec![format!(
                "__rolling_{}_w{}_l{}",
                self.kind.short_name(),
                self.window,
                self.lag
            )]
        }

        fn warmup(&self) -> usize {
            // Block is values[t - lag - window + 1 ..= t - lag].
            // We need t - lag - window + 1 >= 0, i.e. t >= lag + window - 1.
            // warmup returns the minimum valid target_idx for which
            // compute_fit can be called.
            self.lag + self.window - 1
        }

        fn compute_fit(&self, values: &[f64], target_idx: usize, out: &mut [f64]) {
            debug_assert_eq!(out.len(), 1);
            // Block: values[start..end_excl]
            //   end_excl = target_idx - lag + 1   (so last included index is target - lag)
            //   start    = end_excl - window
            // Safe under the warmup contract: target_idx >= lag + window - 1 ≥ lag.
            let end_excl = (target_idx + 1).saturating_sub(self.lag);
            let start = end_excl.saturating_sub(self.window);
            let end = end_excl.min(values.len());
            out[0] = self.kind.compute(&values[start..end]);
        }

        fn compute_predict(&self, recent: &[f64], out: &mut [f64]) {
            debug_assert_eq!(out.len(), 1);
            // At predict step we are about to forecast position P, with
            // recent[len - 1] == values[P - 1]. The block ends at
            // values[P - lag] = recent[len - lag] (inclusive), which is
            // exclusive-end len - lag + 1.
            let len = recent.len();
            let end_excl = (len + 1).saturating_sub(self.lag).min(len);
            let start = end_excl.saturating_sub(self.window);
            out[0] = self.kind.compute(&recent[start..end_excl]);
        }

        fn name(&self) -> &str {
            "RollingFeature"
        }

        fn clone_box(&self) -> Box<dyn RecursiveFeature> {
            Box::new(self.clone())
        }
    }

    // ── Event distance feature ───────────────────────────────────────

    /// Which event-distance columns to emit.
    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
    pub enum EventDistanceMode {
        /// Steps since the most recent event at index ≤ t.
        Since,
        /// Steps until the next event at index > t.
        Until,
        /// Both columns: `[since, until]`.
        Both,
    }

    /// Distance-to-event feature. Given a sorted list of **absolute** event
    /// indices (e.g. holidays, promotions, regime changes), emits
    /// steps-since-last-event and/or steps-until-next-event at each row.
    ///
    /// Pass the full event list spanning training **and** forecast horizon —
    /// at predict step `h`, the distance is computed against `n_train + h`,
    /// so future events must be declared up front.
    ///
    /// When no preceding event exists (or no following event exists), the
    /// `cap` value is emitted (default `usize::MAX as f64`). Use a smaller
    /// `cap` to keep features bounded.
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// use std::sync::Arc;
    /// use anofox_forecast::models::regression::{
    ///     EventDistanceFeature, EventDistanceMode, RegressionFeatures,
    /// };
    ///
    /// // Holidays at indices 30, 60, 120 (covering both train and horizon).
    /// let events = EventDistanceFeature::new(
    ///     vec![30, 60, 120],
    ///     EventDistanceMode::Both,
    /// )
    /// .with_cap(365);
    ///
    /// let features = RegressionFeatures::new()
    ///     .with_recursive(Arc::new(events));
    /// ```
    #[derive(Debug, Clone)]
    pub struct EventDistanceFeature {
        events: Vec<usize>,
        mode: EventDistanceMode,
        cap: f64,
    }

    impl EventDistanceFeature {
        /// Create a new event-distance feature. `events` is sorted ascending
        /// internally; duplicates are de-duplicated.
        pub fn new(mut events: Vec<usize>, mode: EventDistanceMode) -> Self {
            events.sort_unstable();
            events.dedup();
            Self {
                events,
                mode,
                cap: f64::from(u32::MAX),
            }
        }

        /// Override the cap used when no preceding/following event exists.
        pub fn with_cap(mut self, cap: usize) -> Self {
            self.cap = cap as f64;
            self
        }

        /// Sorted unique event indices.
        pub fn events(&self) -> &[usize] {
            &self.events
        }

        /// Output mode.
        pub fn mode(&self) -> EventDistanceMode {
            self.mode
        }

        fn since_at(&self, t: usize) -> f64 {
            // largest event index e such that e ≤ t
            match self.events.binary_search(&t) {
                Ok(_) => 0.0,
                Err(pos) => {
                    if pos == 0 {
                        self.cap
                    } else {
                        (t - self.events[pos - 1]) as f64
                    }
                }
            }
        }

        fn until_at(&self, t: usize) -> f64 {
            // smallest event index e such that e > t
            let pos = self.events.partition_point(|&e| e <= t);
            if pos == self.events.len() {
                self.cap
            } else {
                (self.events[pos] - t) as f64
            }
        }

        fn fill_row(&self, t: usize, out: &mut [f64]) {
            match self.mode {
                EventDistanceMode::Since => {
                    out[0] = self.since_at(t);
                }
                EventDistanceMode::Until => {
                    out[0] = self.until_at(t);
                }
                EventDistanceMode::Both => {
                    out[0] = self.since_at(t);
                    out[1] = self.until_at(t);
                }
            }
        }
    }

    impl RecursiveFeature for EventDistanceFeature {
        fn column_names(&self) -> Vec<String> {
            match self.mode {
                EventDistanceMode::Since => vec!["__event_since".into()],
                EventDistanceMode::Until => vec!["__event_until".into()],
                EventDistanceMode::Both => {
                    vec!["__event_since".into(), "__event_until".into()]
                }
            }
        }

        fn n_columns(&self) -> usize {
            match self.mode {
                EventDistanceMode::Since | EventDistanceMode::Until => 1,
                EventDistanceMode::Both => 2,
            }
        }

        fn warmup(&self) -> usize {
            0
        }

        fn compute_fit(&self, _values: &[f64], target_idx: usize, out: &mut [f64]) {
            self.fill_row(target_idx, out);
        }

        fn compute_predict(&self, recent: &[f64], out: &mut [f64]) {
            // Without n_train context, fall back to using recent.len() as the
            // absolute index. This is only correct when the caller passes the
            // full series as `recent`; the predict-time path uses
            // `compute_predict_at` instead.
            self.fill_row(recent.len(), out);
        }

        fn compute_predict_at(
            &self,
            _recent: &[f64],
            step_h: usize,
            n_train: usize,
            out: &mut [f64],
        ) {
            self.fill_row(n_train + step_h, out);
        }

        fn name(&self) -> &str {
            "EventDistanceFeature"
        }

        fn clone_box(&self) -> Box<dyn RecursiveFeature> {
            Box::new(self.clone())
        }
    }

    // ── Regression backend ────────────────────────────────────────────

    /// Strategy for generating observation weights in WLS.
    #[derive(Debug, Clone)]
    pub enum WeightStrategy {
        /// Exponential decay: `w_i = decay^(n-1-i)`. Most recent observation
        /// gets weight 1.0, oldest gets `decay^(n-1)`.
        ExponentialDecay(f64),
        /// Custom weight vector (must match training row count after lag offset).
        Custom(Vec<f64>),
        /// Logistic (sigmoid) recency weights centred `offset` steps from
        /// the end of the series:
        ///
        /// `w_i = 1 / (1 + exp(-((i as f64) - (n as f64 - offset))))`
        ///
        /// Produces a "use the recent window" effect with smooth boundaries —
        /// the last ~`offset` observations contribute near 1.0, older ones
        /// decay through a single inflection point and plateau near 0.
        /// `offset` is the distance (in observations) from `n` at which the
        /// sigmoid is centred; typical values are 6–24.
        Logistic {
            /// Distance from the end of the series at which the sigmoid is
            /// centred. Must be finite; negative values invert the sigmoid.
            offset: f64,
        },
    }

    /// Specifies which regression estimator backs the forecaster.
    ///
    /// All backends share the same feature engineering ([`RegressionFeatures`]),
    /// design matrix construction, and recursive prediction logic. They differ
    /// only in the loss function / regularization applied during coefficient
    /// estimation.
    ///
    /// # Examples
    ///
    /// ```rust,ignore
    /// use anofox_forecast::models::regression::*;
    ///
    /// // Ridge with λ = 0.1
    /// let model = RegressionForecaster::ridge(0.1, RegressionFeatures::new().trend().fourier(7, 3));
    ///
    /// // Quantile regression at the median
    /// let model = RegressionForecaster::quantile(0.5, RegressionFeatures::new().trend());
    /// ```
    #[derive(Debug, Clone, Default)]
    pub enum RegressionBackend {
        /// Ordinary Least Squares (default).
        #[default]
        Ols,
        /// Ridge regression (L2 regularization).
        Ridge {
            /// Regularization strength (λ ≥ 0).
            lambda: f64,
        },
        /// Elastic Net (L1 + L2 regularization).
        ElasticNet {
            /// Combined regularization strength.
            lambda: f64,
            /// Mixing parameter: 0 = pure Ridge, 1 = pure Lasso.
            alpha: f64,
        },
        /// Quantile regression — estimate a specific conditional quantile.
        Quantile {
            /// Quantile to estimate (0 < τ < 1). 0.5 = median.
            tau: f64,
        },
        /// Weighted Least Squares — down-weight older or less reliable observations.
        Wls {
            /// How observation weights are generated.
            strategy: WeightStrategy,
        },
        /// Combined logistic-recency weighting + L2 regularization. See
        /// [`RegressionForecaster::wls_logistic_ridge`]. Solves the
        /// weighted Ridge normal equations
        /// `β = (Xᵀ W X + λI)⁻¹ Xᵀ W y` in a single fit.
        WlsLogisticRidge {
            /// Logistic sigmoid centre offset (same meaning as
            /// [`WeightStrategy::Logistic`]).
            offset: f64,
            /// L2 penalty strength (λ ≥ 0).
            lambda: f64,
        },
        /// Recursive Least Squares — adaptive coefficients via forgetting factor.
        Rls {
            /// Exponential forgetting (0 < λ ≤ 1). 1.0 = equal weights, <1 = recent emphasis.
            forgetting_factor: f64,
        },
        /// Tweedie GLM — handles count, continuous, and zero-inflated data.
        Tweedie {
            /// Variance power: 0 = Gaussian, 1 = Poisson, 2 = Gamma, 3 = Inv-Gaussian.
            var_power: f64,
            /// Link function power: None = canonical. 0 = log, 1 = identity.
            link_power: Option<f64>,
        },
        /// Poisson GLM — for count data (non-negative integers).
        Poisson,
        /// Bounded Least Squares — box constraints on coefficients.
        Bls {
            /// Lower bound for all coefficients (None = unconstrained).
            lower: Option<f64>,
            /// Upper bound for all coefficients (None = unconstrained).
            upper: Option<f64>,
        },
        /// Dynamic linear model — time-varying parameters via IC-weighted model averaging.
        ///
        /// Automatically generates candidate models from variable subsets, fits each,
        /// and computes observation-level IC weights. Coefficients vary over time,
        /// giving the model the ability to adapt to structural changes.
        Dynamic {
            /// Information criterion for model weighting (default: AICc).
            ic: InformationCriterion,
            /// LOWESS smoothing span for weights (None = no smoothing).
            lowess_span: Option<f64>,
        },
    }

    impl RegressionBackend {
        /// Human-readable name for this backend.
        fn name(&self) -> &str {
            match self {
                Self::Ols => "OLS",
                Self::Ridge { .. } => "Ridge",
                Self::ElasticNet { .. } => "ElasticNet",
                Self::Quantile { .. } => "Quantile",
                Self::Wls { .. } => "WLS",
                Self::WlsLogisticRidge { .. } => "WlsLogisticRidge",
                Self::Rls { .. } => "RLS",
                Self::Tweedie { .. } => "Tweedie",
                Self::Poisson => "Poisson",
                Self::Bls { .. } => "BLS",
                Self::Dynamic { .. } => "Dynamic",
            }
        }

        /// Fit the backend to data, returning a boxed fitted regressor.
        fn fit_to(
            &self,
            x: &Mat<f64>,
            y: &Col<f64>,
        ) -> std::result::Result<Box<dyn FittedRegressor + Send>, String> {
            match self {
                Self::Ols => {
                    let model = OlsRegressor::builder().with_intercept(true).build();
                    let fitted = model.fit(x, y).map_err(|e| format!("OLS: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Ridge { lambda } => {
                    let model = RidgeRegressor::builder()
                        .with_intercept(true)
                        .lambda(*lambda)
                        .build();
                    let fitted = model.fit(x, y).map_err(|e| format!("Ridge: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::ElasticNet { lambda, alpha } => {
                    let model = ElasticNetRegressor::builder()
                        .with_intercept(true)
                        .lambda(*lambda)
                        .alpha(*alpha)
                        .build();
                    let fitted = model.fit(x, y).map_err(|e| format!("ElasticNet: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Quantile { tau } => {
                    let model = QuantileRegressor::builder()
                        .with_intercept(true)
                        .tau(*tau)
                        .build();
                    let fitted = model.fit(x, y).map_err(|e| format!("Quantile: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Wls { strategy } => {
                    let n = y.nrows();
                    let weights = match strategy {
                        WeightStrategy::ExponentialDecay(decay) => {
                            let mut w = Col::zeros(n);
                            for i in 0..n {
                                w[i] = decay.powi((n - 1 - i) as i32);
                            }
                            w
                        }
                        WeightStrategy::Custom(v) => {
                            if v.len() != n {
                                return Err(format!(
                                    "WLS: weight vector length {} != training rows {}",
                                    v.len(),
                                    n
                                ));
                            }
                            let mut w = Col::zeros(n);
                            for (i, &val) in v.iter().enumerate() {
                                w[i] = val;
                            }
                            w
                        }
                        WeightStrategy::Logistic { offset } => {
                            if !offset.is_finite() {
                                return Err(format!(
                                    "WLS Logistic: offset must be finite, got {}",
                                    offset
                                ));
                            }
                            let mut w = Col::zeros(n);
                            let centre = n as f64 - *offset;
                            for i in 0..n {
                                let z = (i as f64) - centre;
                                w[i] = 1.0 / (1.0 + (-z).exp());
                            }
                            w
                        }
                    };
                    let model = WlsRegressor::builder()
                        .with_intercept(true)
                        .weights(weights)
                        .build();
                    let fitted = model.fit(x, y).map_err(|e| format!("WLS: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::WlsLogisticRidge { offset, lambda } => {
                    if !offset.is_finite() {
                        return Err(format!(
                            "WlsLogisticRidge: offset must be finite, got {}",
                            offset
                        ));
                    }
                    if *lambda < 0.0 {
                        return Err(format!(
                            "WlsLogisticRidge: lambda must be non-negative, got {}",
                            lambda
                        ));
                    }
                    let n = y.nrows();
                    let p = x.ncols();
                    // Build logistic weights and accumulators for weighted means.
                    let centre = n as f64 - *offset;
                    let mut w = vec![0.0_f64; n];
                    let mut sum_w = 0.0;
                    for i in 0..n {
                        let z = (i as f64) - centre;
                        let wi = 1.0 / (1.0 + (-z).exp());
                        w[i] = wi;
                        sum_w += wi;
                    }
                    if sum_w == 0.0 {
                        return Err("WlsLogisticRidge: weights summed to zero".into());
                    }
                    let mut x_bar_w = vec![0.0_f64; p];
                    let mut y_bar_w = 0.0;
                    for i in 0..n {
                        y_bar_w += w[i] * y[i];
                        for j in 0..p {
                            x_bar_w[j] += w[i] * x[(i, j)];
                        }
                    }
                    y_bar_w /= sum_w;
                    for v in x_bar_w.iter_mut() {
                        *v /= sum_w;
                    }
                    // Augmented design: top n rows are sqrt(W) * (X − X̄_w),
                    // bottom p rows are √λ · I_p (Tikhonov penalty, intercept
                    // not augmented since we re-add it analytically below).
                    let sqrt_lambda = lambda.sqrt();
                    let n_aug = n + p;
                    let mut x_aug = Mat::zeros(n_aug, p);
                    let mut y_aug = Col::zeros(n_aug);
                    for i in 0..n {
                        let s = w[i].sqrt();
                        for j in 0..p {
                            x_aug[(i, j)] = s * (x[(i, j)] - x_bar_w[j]);
                        }
                        y_aug[i] = s * (y[i] - y_bar_w);
                    }
                    for j in 0..p {
                        x_aug[(n + j, j)] = sqrt_lambda;
                    }
                    // OLS on augmented system → β solves
                    // (Xᵀ W X + λI) β = Xᵀ W (y − ȳ_w · 1), which after
                    // re-adding ȳ_w to the intercept is the weighted Ridge
                    // solution.
                    let ols = OlsRegressor::builder().with_intercept(false).build();
                    let ols_fitted = ols
                        .fit(&x_aug, &y_aug)
                        .map_err(|e| format!("WlsLogisticRidge: {}", e))?;
                    let beta = &ols_fitted.result().coefficients;
                    // Reconstruct the intercept and predictions on the
                    // original (unweighted, uncentred) scale.
                    let mut intercept = y_bar_w;
                    for j in 0..p {
                        intercept -= x_bar_w[j] * beta[j];
                    }
                    let mut fitted_values = Col::zeros(n);
                    let mut residuals = Col::zeros(n);
                    let mut rss = 0.0;
                    for i in 0..n {
                        let mut pred = intercept;
                        for j in 0..p {
                            pred += beta[j] * x[(i, j)];
                        }
                        fitted_values[i] = pred;
                        let r = y[i] - pred;
                        residuals[i] = r;
                        rss += r * r;
                    }
                    let y_mean: f64 = (0..n).map(|i| y[i]).sum::<f64>() / n as f64;
                    let tss: f64 = (0..n).map(|i| (y[i] - y_mean).powi(2)).sum();
                    let r_squared = if tss > 0.0 { 1.0 - rss / tss } else { 0.0 };

                    // Compute standard errors on the original-scale design.
                    // The augmented OLS that produced β was fitted on a
                    // centred + √W-scaled + Tikhonov-augmented design, so
                    // its result.std_errors are in the wrong units and
                    // must NOT be propagated. Instead, populate from the
                    // weighted covariance
                    //   Cov(β̂) = σ̂² · (Xfullᵀ W Xfull + λI_β)⁻¹
                    //   σ̂² = Σ wᵢ rᵢ² / (Σ wᵢ − (p + 1))
                    // where Xfull prepends the intercept column and
                    // λI_β has 0 on the intercept diagonal (intercept is
                    // unpenalised). This is the "plain weighted Ridge"
                    // form from issue #115 — gives correct SEs at λ = 0
                    // and a numerically-stable approximation for λ > 0.
                    let p_full = p + 1;
                    let mut xtwx: Mat<f64> = Mat::zeros(p_full, p_full);
                    let mut sum_w_rsq = 0.0;
                    for i in 0..n {
                        let wi = w[i];
                        sum_w_rsq += wi * residuals[i] * residuals[i];
                        xtwx[(0, 0)] += wi;
                        for j in 0..p {
                            let wij = wi * x[(i, j)];
                            xtwx[(0, j + 1)] += wij;
                            xtwx[(j + 1, 0)] += wij;
                            for k in 0..p {
                                xtwx[(j + 1, k + 1)] += wij * x[(i, k)];
                            }
                        }
                    }
                    for j in 0..p {
                        xtwx[(j + 1, j + 1)] += *lambda;
                    }
                    let df_eff = sum_w - p_full as f64;
                    let sigma_sq = if df_eff > 0.0 {
                        sum_w_rsq / df_eff
                    } else {
                        f64::NAN
                    };
                    let (std_errors_opt, intercept_se_opt) = match xtwx.llt(faer::Side::Lower) {
                        Ok(llt) => {
                            use faer::linalg::solvers::DenseSolveCore;
                            let inv: Mat<f64> = llt.inverse();
                            let mut se_coef = Col::<f64>::zeros(p);
                            for j in 0..p {
                                let var = sigma_sq * inv[(j + 1, j + 1)];
                                se_coef[j] = if var.is_finite() && var >= 0.0 {
                                    var.sqrt()
                                } else {
                                    f64::NAN
                                };
                            }
                            let var_int = sigma_sq * inv[(0, 0)];
                            let se_int = if var_int.is_finite() && var_int >= 0.0 {
                                var_int.sqrt()
                            } else {
                                f64::NAN
                            };
                            (Some(se_coef), Some(se_int))
                        }
                        Err(_) => (None, None),
                    };

                    let mut result = ols_fitted.result().clone();
                    result.intercept = Some(intercept);
                    result.fitted_values = fitted_values;
                    result.residuals = residuals;
                    result.n_observations = n;
                    result.r_squared = r_squared;
                    result.std_errors = std_errors_opt;
                    result.intercept_std_error = intercept_se_opt;
                    Ok(Box::new(WlsLogisticRidgeFitted { result }))
                }
                Self::Rls { forgetting_factor } => {
                    let model = RlsRegressor::builder()
                        .with_intercept(true)
                        .forgetting_factor(*forgetting_factor)
                        .build();
                    let fitted = model.fit(x, y).map_err(|e| format!("RLS: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Tweedie {
                    var_power,
                    link_power,
                } => {
                    let mut builder = TweedieRegressor::builder()
                        .with_intercept(true)
                        .var_power(*var_power);
                    if let Some(lp) = link_power {
                        builder = builder.link_power(*lp);
                    }
                    let model = builder.build();
                    let fitted = model.fit(x, y).map_err(|e| format!("Tweedie: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Poisson => {
                    let model = PoissonRegressor::builder().with_intercept(true).build();
                    let fitted = model.fit(x, y).map_err(|e| format!("Poisson: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Bls { lower, upper } => {
                    let mut builder = BlsRegressor::builder().with_intercept(true);
                    if let Some(lb) = lower {
                        builder = builder.lower_bound_all(*lb);
                    }
                    if let Some(ub) = upper {
                        builder = builder.upper_bound_all(*ub);
                    }
                    let model = builder.build();
                    let fitted = model.fit(x, y).map_err(|e| format!("BLS: {}", e))?;
                    Ok(Box::new(fitted))
                }
                Self::Dynamic { ic, lowess_span } => {
                    let mut builder = LmDynamicRegressor::builder().with_intercept(true).ic(*ic);
                    if let Some(span) = lowess_span {
                        builder = builder.lowess_span(*span);
                    } else {
                        builder = builder.no_smoothing();
                    }
                    let model = builder.build();
                    let fitted = model.fit(x, y).map_err(|e| format!("Dynamic: {}", e))?;
                    Ok(Box::new(fitted))
                }
            }
        }
    }

    // ── Component specifications ────────────────────────────────────

    /// Specifies a trend model to include as a regression feature column.
    ///
    /// Each trend type produces **one column** in the design matrix containing
    /// the fitted trend values. During prediction, `predict_trend(horizon)` is
    /// used to generate future values.
    ///
    /// All trend components are fitted with `Recency::Full` so the regression
    /// model sees the component's view of the entire training window.
    ///
    /// # Cross-validation caveat
    ///
    /// Trend features are functions of the training data only. If you use
    /// [`Recency::Auto`] (changepoint-based) on a manually created component,
    /// ensure the changepoint detection only sees the training fold — otherwise
    /// the changepoint location leaks future information into the features.
    #[derive(Debug, Clone)]
    pub enum TrendType {
        /// Linear trend via [`PolynomialTrend`](crate::seasonality::PolynomialTrend) degree 1.
        Linear,
        /// Quadratic trend via [`PolynomialTrend`](crate::seasonality::PolynomialTrend) degree 2.
        Quadratic,
        /// Cubic trend via [`PolynomialTrend`](crate::seasonality::PolynomialTrend) degree 3.
        Cubic,
        /// Exponential trend via [`ExponentialTrend`](crate::seasonality::ExponentialTrend).
        /// Requires positive values.
        Exponential,
        /// Theil-Sen robust linear trend via [`TheilSenTrend`](crate::seasonality::TheilSenTrend).
        TheilSen,
    }

    impl TrendType {
        /// Feature safety classification.
        pub fn safety(&self) -> FeatureSafety {
            FeatureSafety::DataDependent
        }
    }

    /// Specifies a seasonal component to include as regression feature column(s).
    ///
    /// - [`Fourier`](SeasonalSpec::Fourier): `2 * order` sin/cos columns — deterministic
    ///   functions of the time index, no fitting required (Prophet-style).
    /// - [`Dummy`](SeasonalSpec::Dummy): 1 column of per-period seasonal means,
    ///   fitted from training data.
    #[derive(Debug, Clone)]
    pub enum SeasonalSpec {
        /// Fourier seasonality: `2 * order` sin/cos columns at the given period.
        ///
        /// Period is in observation units (e.g., 7 for weekly with daily data,
        /// 12 for yearly with monthly data).
        Fourier {
            /// Seasonal period in observation units.
            period: usize,
            /// Number of Fourier pairs (total columns = 2 * order).
            order: usize,
        },
        /// Dummy seasonal encoding with the given period.
        ///
        /// One column of per-position means (averaged over all full cycles
        /// in the training data).
        Dummy(usize),
    }

    impl SeasonalSpec {
        /// Feature safety classification.
        pub fn safety(&self) -> FeatureSafety {
            match self {
                SeasonalSpec::Fourier { .. } => FeatureSafety::Deterministic,
                SeasonalSpec::Dummy(_) => FeatureSafety::DataDependent,
            }
        }
    }

    // ── Fitted component storage ────────────────────────────────────

    /// Internal storage for a fitted trend/seasonal component.
    #[derive(Debug)]
    enum FittedComponentState {
        Polynomial(PolynomialTrend),
        Exponential(ExponentialTrend),
        TheilSen(TheilSenTrend),
        Dummy(DummySeasonality),
        Fourier {
            period: usize,
            order: usize,
        },
        /// Structural feature — forward-filled during prediction.
        Structural {
            /// Prediction fill values — one per column, repeated for every forecast step.
            fill_values: Vec<f64>,
        },
    }

    impl FittedComponentState {
        /// Number of columns this component contributes to the design matrix.
        fn n_columns(&self) -> usize {
            match self {
                Self::Polynomial(_) | Self::Exponential(_) | Self::TheilSen(_) | Self::Dummy(_) => {
                    1
                }
                Self::Fourier { order, .. } => 2 * order,
                Self::Structural { fill_values } => fill_values.len(),
            }
        }

        /// Generate future feature columns for this component.
        ///
        /// Returns one `Vec<f64>` per column (most components produce 1 column,
        /// Fourier produces `2 * order`).
        fn predict(&self, horizon: usize, n_train: usize) -> Vec<Vec<f64>> {
            match self {
                Self::Polynomial(p) => vec![p.predict_trend(horizon)],
                Self::Exponential(e) => vec![e.predict_trend(horizon)],
                Self::TheilSen(t) => vec![t.predict_trend(horizon)],
                Self::Dummy(d) => vec![d.predict_seasonal(horizon)],
                Self::Fourier { period, order } => {
                    let timestamps: Vec<f64> = (0..horizon).map(|h| (n_train + h) as f64).collect();
                    fourier_terms(&timestamps, *period as f64, *order).unwrap_or_default()
                }
                Self::Structural { fill_values } => {
                    fill_values.iter().map(|&v| vec![v; horizon]).collect()
                }
            }
        }
    }

    // ── Feature specification ───────────────────────────────────────

    /// Criterion for automatic lag selection.
    #[derive(Debug, Clone, Copy, PartialEq, Eq)]
    pub enum LagSelectionCriterion {
        /// Bayesian Information Criterion (default). Penalizes complexity more than AIC.
        Bic,
        /// Akaike Information Criterion. Penalizes complexity less than BIC.
        Aic,
    }

    /// Spec for derived exogenous-regressor features.
    ///
    /// Allows declaring lagged or rolling-window transforms of exog columns
    /// without preprocessing the input series. The model materializes the
    /// columns at fit-time and resolves them at predict-time using the
    /// stored exog tail plus the user-supplied future regressors.
    #[derive(Debug, Clone)]
    pub enum ExogFeatureSpec {
        /// Lagged values of the named exog column at the given lag offsets.
        ///
        /// Each entry in `lags` produces one column `{col}_lag{k}`.
        Lag {
            /// Exog column name to lag.
            col: String,
            /// Lag offsets (must be ≥ 1; `0` would be the contemporaneous value).
            lags: Vec<usize>,
        },
        /// Rolling-window statistic of the named exog column.
        ///
        /// Window ends at `t - lag` (inclusive) and spans `window` observations.
        /// Produces one column `{col}_roll_{kind}_w{window}_l{lag}`.
        Rolling {
            /// Exog column name to roll.
            col: String,
            /// Window size (≥ 1).
            window: usize,
            /// Lag offset of the window end (≥ 0; `0` includes the
            /// contemporaneous value, which is fine for *exogenous* signals
            /// known at time t).
            lag: usize,
            /// Statistic kind.
            kind: RollingStatKind,
        },
        /// Polynomial expansion of the named exog column: emits `col^k` for
        /// `k = 2..=degree`. Use `degree = 1` for no-op (no columns).
        ///
        /// Each emitted column is named `{col}_pow{k}`.
        Polynomial {
            /// Exog column name to expand.
            col: String,
            /// Maximum power. Must be ≥ 2 to produce any columns.
            degree: usize,
        },
        /// Pairwise interaction `col_a × col_b` between two exog columns.
        ///
        /// Produces one column named `{col_a}_x_{col_b}`.
        Interaction {
            /// First exog column name.
            col_a: String,
            /// Second exog column name.
            col_b: String,
        },
        /// Categorical encoding of an integer-coded exog column.
        ///
        /// The exog column carries integer category codes (stored as `f64`,
        /// rounded at fit/predict to nearest `i64`). Categories not in the
        /// declared `categories` list are treated as **unknown** at both fit
        /// and predict time; the unknown handling depends on the strategy:
        ///
        /// - `OneHot`: all dummy columns are zero (interpreted as the
        ///   `drop_first` baseline if dropped).
        /// - `Ordinal`: encoded as `categories.len()` (the next unused index).
        /// - `Count`: encoded as `0` (no observations seen).
        /// - `Target`: encoded as the grand mean of the target.
        ///
        /// Predeclaring categories keeps the design-matrix shape deterministic
        /// across CV folds and prevents leakage from train→test category
        /// drift.
        Categorical {
            /// Exog column name to encode.
            col: String,
            /// Declared category codes. Order matters (used by `OneHot` and
            /// `Ordinal`).
            categories: Vec<i64>,
            /// Encoding strategy.
            strategy: CategoricalStrategy,
        },
    }

    /// Fitted categorical encoder for `Count` / `Target` strategies.
    ///
    /// Stateless strategies (`OneHot`, `Ordinal`) do not produce one of these
    /// — their values are reconstructible from the spec alone.
    #[derive(Debug, Clone)]
    pub(crate) struct FittedCategoricalEncoder {
        /// Per-category encoded value (count for Count, mean-with-smoothing for Target).
        map: HashMap<i64, f64>,
        /// Value used for unseen categories.
        fallback: f64,
    }

    /// Encoding strategy for [`ExogFeatureSpec::Categorical`].
    #[derive(Debug, Clone, Copy, PartialEq)]
    pub enum CategoricalStrategy {
        /// Dummy / one-hot encoding. With `drop_first = true`, the first
        /// declared category is omitted (baseline reference).
        OneHot {
            /// Drop the first declared category to avoid the dummy-variable
            /// trap with an intercept term.
            drop_first: bool,
        },
        /// Ordinal — emit the position of the value in `categories`. Unknown
        /// codes map to `categories.len()`.
        Ordinal,
        /// Count encoding — emit the number of fit-time observations of each
        /// category. Unknown codes map to 0.
        Count,
        /// Target (mean) encoding with additive smoothing.
        ///
        /// `encoded(c) = (n_c · mean_c + smoothing · grand_mean) / (n_c + smoothing)`
        ///
        /// Higher `smoothing` shrinks rare-category encodings toward the
        /// grand mean. `smoothing = 0` is plain group-mean encoding.
        Target {
            /// Additive prior weight (≥ 0). 0 = no smoothing.
            smoothing: f64,
        },
    }

    impl ExogFeatureSpec {
        /// All exog column names this spec references (1 or 2).
        fn referenced_cols(&self) -> Vec<&str> {
            match self {
                Self::Lag { col, .. }
                | Self::Rolling { col, .. }
                | Self::Polynomial { col, .. }
                | Self::Categorical { col, .. } => vec![col],
                Self::Interaction { col_a, col_b } => vec![col_a, col_b],
            }
        }

        fn warmup(&self) -> usize {
            match self {
                Self::Lag { lags, .. } => lags.iter().copied().max().unwrap_or(0),
                Self::Rolling { window, lag, .. } => lag + window.saturating_sub(1),
                Self::Polynomial { .. } | Self::Interaction { .. } | Self::Categorical { .. } => 0,
            }
        }

        #[allow(dead_code)]
        fn n_columns(&self) -> usize {
            match self {
                Self::Lag { lags, .. } => lags.len(),
                Self::Rolling { .. } | Self::Interaction { .. } => 1,
                Self::Polynomial { degree, .. } => degree.saturating_sub(1),
                Self::Categorical {
                    categories,
                    strategy,
                    ..
                } => match strategy {
                    CategoricalStrategy::OneHot { drop_first } => {
                        if *drop_first {
                            categories.len().saturating_sub(1)
                        } else {
                            categories.len()
                        }
                    }
                    CategoricalStrategy::Ordinal
                    | CategoricalStrategy::Count
                    | CategoricalStrategy::Target { .. } => 1,
                },
            }
        }

        fn column_names(&self) -> Vec<String> {
            match self {
                Self::Lag { col, lags } => {
                    lags.iter().map(|k| format!("{}_lag{}", col, k)).collect()
                }
                Self::Rolling {
                    col,
                    window,
                    lag,
                    kind,
                } => vec![format!(
                    "{}_roll_{}_w{}_l{}",
                    col,
                    kind.short_name(),
                    window,
                    lag
                )],
                Self::Polynomial { col, degree } => {
                    (2..=*degree).map(|k| format!("{}_pow{}", col, k)).collect()
                }
                Self::Interaction { col_a, col_b } => {
                    vec![format!("{}_x_{}", col_a, col_b)]
                }
                Self::Categorical {
                    col,
                    categories,
                    strategy,
                } => match strategy {
                    CategoricalStrategy::OneHot { drop_first } => {
                        let start = if *drop_first { 1 } else { 0 };
                        categories[start..]
                            .iter()
                            .map(|cat| format!("{}_eq_{}", col, cat))
                            .collect()
                    }
                    CategoricalStrategy::Ordinal => vec![format!("{}_ord", col)],
                    CategoricalStrategy::Count => vec![format!("{}_count", col)],
                    CategoricalStrategy::Target { .. } => vec![format!("{}_target", col)],
                },
            }
        }
    }

    /// Configures which features are built from a [`TimeSeries`] for
    /// the regression model.
    ///
    /// Features are added to the design matrix in this order:
    /// 1. Linear trend index (if `use_trend`)
    /// 2. Autoregressive lags
    /// 3. Trend component columns ([`TrendType`])
    /// 4. Seasonal component columns ([`SeasonalSpec`])
    /// 5. Structural feature columns ([`StructuralFeature`])
    /// 6. Recursive feature columns ([`RecursiveFeature`], e.g. [`RollingFeature`])
    /// 7. Exogenous regressors (if `use_exog`)
    /// 8. Derived exog features ([`ExogFeatureSpec`], lag/rolling)
    #[derive(Debug, Clone)]
    pub struct RegressionFeatures {
        /// Include a linear trend index (0, 1, …, n-1).
        pub use_trend: bool,
        /// Number of autoregressive lags to include (legacy; prefer `lag_indices`).
        pub max_lag: usize,
        /// Specific lag indices to include (e.g. `vec![12]` for only lag-12).
        /// When non-empty, this takes precedence over `max_lag`.
        pub lag_indices: Vec<usize>,
        /// Automatic lag selection: (max_lag, criterion). When set, overrides max_lag/lag_indices.
        pub auto_lag_config: Option<(usize, LagSelectionCriterion)>,
        /// Include exogenous regressors from the TimeSeries (if present).
        pub use_exog: bool,
        /// Trend components to include as feature columns.
        pub trend_components: Vec<TrendType>,
        /// Seasonal components to include as feature columns.
        pub seasonal_components: Vec<SeasonalSpec>,
        /// Structural features (forward-filled during prediction).
        pub structural_features: Vec<Arc<dyn StructuralFeature>>,
        /// Recursive features (recomputed from the rolling history buffer at
        /// every horizon step — e.g. [`RollingFeature`]).
        pub recursive_features: Vec<Arc<dyn RecursiveFeature>>,
        /// Derived exog features (lag / rolling transforms of named exog
        /// columns). Materialized after the contemporaneous exog block.
        pub exog_features: Vec<ExogFeatureSpec>,
        /// Regular differencing order (d). Applied before fitting, integrated after predict.
        pub diff_order: usize,
        /// Seasonal differencing specs: Vec of (order D, period s). Applied in order before fitting, integrated in reverse after predict.
        pub seasonal_diffs: Vec<(usize, usize)>,
        /// Fractional differencing order (0 < d < 1). Applied before fitting.
        /// Note: fractional differencing is NOT invertible for prediction —
        /// forecasts are on the differenced scale. Use for feature engineering
        /// in regression, not for direct forecasting.
        pub frac_diff_order: Option<f64>,
    }

    impl Default for RegressionFeatures {
        fn default() -> Self {
            Self {
                use_trend: true,
                max_lag: 0,
                lag_indices: Vec::new(),
                auto_lag_config: None,
                use_exog: true,
                trend_components: Vec::new(),
                seasonal_components: Vec::new(),
                structural_features: Vec::new(),
                recursive_features: Vec::new(),
                exog_features: Vec::new(),
                diff_order: 0,
                seasonal_diffs: Vec::new(),
                frac_diff_order: None,
            }
        }
    }

    impl RegressionFeatures {
        /// Create a new feature configuration (trend only by default).
        pub fn new() -> Self {
            Self::default()
        }

        /// Include a linear trend index.
        pub fn trend(mut self) -> Self {
            self.use_trend = true;
            self
        }

        /// Do not include a trend index.
        pub fn no_trend(mut self) -> Self {
            self.use_trend = false;
            self
        }

        /// Include autoregressive lags `y[t-1] … y[t-max_lag]`.
        pub fn lags(mut self, max_lag: usize) -> Self {
            self.max_lag = max_lag;
            self.lag_indices = Vec::new(); // clear specific lags
            self
        }

        /// Automatically select the best lag order up to `max_lag` using BIC.
        ///
        /// Tries each lag order from 0 to `max_lag`, fits OLS, computes BIC, and
        /// selects the order with the lowest BIC. This is resolved during `fit()`.
        pub fn auto_lags(mut self, max_lag: usize) -> Self {
            self.auto_lag_config = Some((max_lag, LagSelectionCriterion::Bic));
            self.max_lag = 0;
            self.lag_indices = Vec::new();
            self
        }

        /// Automatically select the best lag order using the specified criterion.
        pub fn auto_lags_with(mut self, max_lag: usize, criterion: LagSelectionCriterion) -> Self {
            self.auto_lag_config = Some((max_lag, criterion));
            self.max_lag = 0;
            self.lag_indices = Vec::new();
            self
        }

        /// Include only the specified lag indices (e.g. `&[12]` for lag-12 only).
        pub fn specific_lags(mut self, lags: &[usize]) -> Self {
            let mut sorted = lags.to_vec();
            sorted.sort_unstable();
            sorted.dedup();
            self.lag_indices = sorted;
            self.max_lag = 0; // specific_lags takes precedence
            self
        }

        /// Return the effective list of lag indices.
        fn effective_lags(&self) -> Vec<usize> {
            if !self.lag_indices.is_empty() {
                self.lag_indices.clone()
            } else if self.max_lag > 0 {
                (1..=self.max_lag).collect()
            } else {
                Vec::new()
            }
        }

        /// Maximum lag value (for offset calculation and tail storage).
        fn max_effective_lag(&self) -> usize {
            if !self.lag_indices.is_empty() {
                self.lag_indices.iter().copied().max().unwrap_or(0)
            } else {
                self.max_lag
            }
        }

        /// Resolve auto-lag selection if configured.
        ///
        /// Tries each lag order 0..=max_lag, builds the feature matrix with that
        /// lag order, fits OLS, computes the information criterion, and selects
        /// the order with the lowest IC. Mutates `self.max_lag` with the result.
        fn resolve_auto_lags(&mut self, series: &TimeSeries) -> Result<()> {
            let (max_lag, criterion) = match self.auto_lag_config {
                Some(cfg) => cfg,
                None => return Ok(()),
            };

            let values = series.primary_values();
            let n = values.len();

            if max_lag == 0 || n < 4 {
                self.max_lag = 0;
                self.auto_lag_config = None;
                return Ok(());
            }

            let mut best_ic = f64::INFINITY;
            let mut best_order = 0_usize;

            for p in 0..=max_lag {
                // Temporarily set lag order
                self.max_lag = p;
                self.lag_indices.clear();

                let offset = self.lag_offset();
                if n <= offset + 1 {
                    continue; // Not enough data for this lag order
                }

                // Build matrices and fit OLS
                let result = self.build_matrices(series);
                let (x, y, n_train, _, _, _) = match result {
                    Ok(v) => v,
                    Err(_) => continue,
                };

                if n_train < 3 {
                    continue;
                }

                let fitted = match RegressionBackend::Ols.fit_to(&x, &y) {
                    Ok(f) => f,
                    Err(_) => continue,
                };

                // Compute RSS from in-sample predictions
                let preds = fitted.predict(&x);
                let mut rss = 0.0_f64;
                for i in 0..n_train {
                    let r = y[i] - preds[i];
                    rss += r * r;
                }

                if !rss.is_finite() || rss <= 0.0 {
                    continue;
                }

                let n_f = n_train as f64;
                let k = x.ncols() as f64 + 1.0; // +1 for intercept

                let ic = match criterion {
                    LagSelectionCriterion::Bic => n_f * (rss / n_f).ln() + k * n_f.ln(),
                    LagSelectionCriterion::Aic => n_f * (rss / n_f).ln() + 2.0 * k,
                };

                if ic.is_finite() && ic < best_ic {
                    best_ic = ic;
                    best_order = p;
                }
            }

            // Set the winning lag order (keep auto_lag_config so re-fit works)
            self.max_lag = best_order;
            self.lag_indices.clear();

            Ok(())
        }

        /// Include exogenous regressors from the TimeSeries.
        pub fn exog(mut self) -> Self {
            self.use_exog = true;
            self
        }

        /// Do not include exogenous regressors.
        pub fn no_exog(mut self) -> Self {
            self.use_exog = false;
            self
        }

        /// Add a trend component as a regression feature column.
        pub fn with_trend_component(mut self, trend: TrendType) -> Self {
            self.trend_components.push(trend);
            self
        }

        /// Add a seasonal component as regression feature column(s).
        pub fn with_seasonal(mut self, seasonal: SeasonalSpec) -> Self {
            self.seasonal_components.push(seasonal);
            self
        }

        /// Add Fourier seasonality features (shorthand for
        /// `with_seasonal(SeasonalSpec::Fourier { period, order })`).
        pub fn fourier(self, period: usize, order: usize) -> Self {
            self.with_seasonal(SeasonalSpec::Fourier { period, order })
        }

        /// Add dummy seasonal encoding (shorthand for
        /// `with_seasonal(SeasonalSpec::Dummy(period))`).
        pub fn dummy_seasonal(self, period: usize) -> Self {
            self.with_seasonal(SeasonalSpec::Dummy(period))
        }

        /// Add a structural feature (forward-filled during prediction).
        pub fn with_structural(mut self, feature: Arc<dyn StructuralFeature>) -> Self {
            self.structural_features.push(feature);
            self
        }

        /// Add changepoint step-function features (convenience).
        pub fn with_changepoint_steps(self, indices: Vec<usize>) -> Self {
            self.with_structural(Arc::new(ChangepointFeature::step_functions(indices)))
        }

        /// Add changepoint features with specified encoding (convenience).
        pub fn with_changepoints(self, indices: Vec<usize>, encoding: ChangepointEncoding) -> Self {
            self.with_structural(Arc::new(ChangepointFeature::new(indices, encoding)))
        }

        /// Add a recursive feature (recomputed at every horizon step during
        /// multi-step prediction — see [`RecursiveFeature`]).
        pub fn with_recursive(mut self, feature: Arc<dyn RecursiveFeature>) -> Self {
            self.recursive_features.push(feature);
            self
        }

        /// Add a rolling statistic as a feature.
        ///
        /// Uses the default `lag = 1` (the window excludes the current target).
        ///
        /// Returns an error if the rolling feature parameters are invalid.
        pub fn with_rolling(self, window: usize, kind: RollingStatKind) -> Result<Self> {
            let feat = RollingFeature::new(window, kind)?;
            Ok(self.with_recursive(Arc::new(feat)))
        }

        /// Add a rolling statistic with an explicit lag.
        ///
        /// `lag >= 1` is required (target-leakage guard).
        pub fn with_rolling_lagged(
            self,
            window: usize,
            lag: usize,
            kind: RollingStatKind,
        ) -> Result<Self> {
            let feat = RollingFeature::with_lag(window, lag, kind)?;
            Ok(self.with_recursive(Arc::new(feat)))
        }

        /// Add a rolling mean feature (shorthand for
        /// `with_rolling(window, RollingStatKind::Mean)`).
        pub fn with_rolling_mean(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Mean)
        }

        /// Add a rolling sample standard deviation feature.
        pub fn with_rolling_std(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Std)
        }

        /// Add a rolling sample variance feature.
        pub fn with_rolling_var(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Var)
        }

        /// Add a rolling minimum feature.
        pub fn with_rolling_min(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Min)
        }

        /// Add a rolling maximum feature.
        pub fn with_rolling_max(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Max)
        }

        /// Add a rolling median feature.
        pub fn with_rolling_median(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Median)
        }

        /// Add a rolling sum feature.
        pub fn with_rolling_sum(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Sum)
        }

        /// Add an exponentially weighted mean feature with smoothing factor
        /// `alpha ∈ (0, 1]` over a window of `window` observations.
        pub fn with_ewm_mean(self, window: usize, alpha: f64) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::EwmMean { alpha })
        }

        /// Add an exponentially weighted standard deviation feature.
        pub fn with_ewm_std(self, window: usize, alpha: f64) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::EwmStd { alpha })
        }

        /// Add a rolling τ-quantile feature. `tau ∈ [0, 1]`.
        pub fn with_rolling_quantile(self, window: usize, tau: f64) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Quantile { tau })
        }

        /// Add a rolling range (max − min) feature.
        pub fn with_rolling_range(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Range)
        }

        /// Add a rolling interquartile range (Q₀.₇₅ − Q₀.₂₅) feature.
        pub fn with_rolling_iqr(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Iqr)
        }

        /// Add a rolling skewness feature.
        pub fn with_rolling_skew(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Skew)
        }

        /// Add a rolling excess-kurtosis feature.
        pub fn with_rolling_kurt(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Kurt)
        }

        /// Add a rolling OLS slope feature (local trend coefficient).
        pub fn with_rolling_slope(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Slope)
        }

        /// Add a rolling fractional-rank feature for the most recent value.
        pub fn with_rolling_rank(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::Rank)
        }

        /// Add a rolling z-score feature for the most recent value relative
        /// to its window predecessors.
        pub fn with_rolling_zscore(self, window: usize) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::ZScore)
        }

        /// Add a rolling count of values strictly greater than `threshold`.
        pub fn with_rolling_count_above(self, window: usize, threshold: f64) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::CountAbove { threshold })
        }

        /// Add a rolling count of values strictly less than `threshold`.
        pub fn with_rolling_count_below(self, window: usize, threshold: f64) -> Result<Self> {
            self.with_rolling(window, RollingStatKind::CountBelow { threshold })
        }

        /// Add an event-distance feature (steps-since-last and/or
        /// steps-until-next event) — see [`EventDistanceFeature`].
        ///
        /// Pass the full event index list spanning training **and** the
        /// forecast horizon: at predict step `h`, the distance is computed
        /// against absolute index `n_train + h`.
        pub fn with_event_distance(self, events: Vec<usize>, mode: EventDistanceMode) -> Self {
            self.with_recursive(Arc::new(EventDistanceFeature::new(events, mode)))
        }

        /// Add lagged copies of an exogenous regressor column.
        ///
        /// Each entry in `lags` produces one column `{col}_lag{k}`.
        /// Requires `use_exog = true` for predict-time column lookup.
        pub fn with_exog_lags(mut self, col: impl Into<String>, lags: &[usize]) -> Self {
            self.exog_features.push(ExogFeatureSpec::Lag {
                col: col.into(),
                lags: lags.to_vec(),
            });
            self
        }

        /// Add a rolling-window statistic over an exogenous regressor column.
        ///
        /// Window ends at `t - lag` (inclusive). `lag = 0` is allowed for
        /// exog signals (unlike target rolling), since exog values at time t
        /// are observed by definition.
        pub fn with_exog_rolling(
            mut self,
            col: impl Into<String>,
            window: usize,
            lag: usize,
            kind: RollingStatKind,
        ) -> Result<Self> {
            if window == 0 {
                return Err(ForecastError::InvalidParameter(
                    "exog_rolling: window must be >= 1".into(),
                ));
            }
            if let RollingStatKind::Quantile { tau } = kind {
                if !(0.0..=1.0).contains(&tau) || tau.is_nan() {
                    return Err(ForecastError::InvalidParameter(format!(
                        "exog_rolling: Quantile tau must satisfy 0 ≤ τ ≤ 1, got {}",
                        tau
                    )));
                }
            }
            if let RollingStatKind::EwmMean { alpha } | RollingStatKind::EwmStd { alpha } = kind {
                if !(0.0 < alpha && alpha <= 1.0) {
                    return Err(ForecastError::InvalidParameter(format!(
                        "exog_rolling: EWM alpha must satisfy 0 < α ≤ 1, got {}",
                        alpha
                    )));
                }
            }
            self.exog_features.push(ExogFeatureSpec::Rolling {
                col: col.into(),
                window,
                lag,
                kind,
            });
            Ok(self)
        }

        /// Add a polynomial expansion of an exogenous regressor column.
        ///
        /// Emits one column per power `k = 2..=degree`, named
        /// `{col}_pow{k}`. `degree = 1` is a no-op (no columns added);
        /// `degree = 0` is rejected.
        ///
        /// Use to capture nonlinear effects of an exog signal in linear
        /// regression backends (OLS, Ridge, ElasticNet).
        pub fn with_exog_polynomial(
            mut self,
            col: impl Into<String>,
            degree: usize,
        ) -> Result<Self> {
            if degree == 0 {
                return Err(ForecastError::InvalidParameter(
                    "exog_polynomial: degree must be >= 1".into(),
                ));
            }
            self.exog_features.push(ExogFeatureSpec::Polynomial {
                col: col.into(),
                degree,
            });
            Ok(self)
        }

        /// Add a pairwise interaction between two exogenous regressor
        /// columns. Emits one column `{col_a}_x_{col_b}` containing the
        /// elementwise product.
        pub fn with_exog_interaction(
            mut self,
            col_a: impl Into<String>,
            col_b: impl Into<String>,
        ) -> Self {
            self.exog_features.push(ExogFeatureSpec::Interaction {
                col_a: col_a.into(),
                col_b: col_b.into(),
            });
            self
        }

        /// Add a categorical encoding of an integer-coded exog column.
        ///
        /// Categories must be declared up front so the design-matrix shape
        /// is deterministic across CV folds. Values not in `categories` are
        /// treated as **unknown** (see [`ExogFeatureSpec::Categorical`]).
        ///
        /// # Errors
        ///
        /// Returns an error if `categories` is empty, contains duplicates,
        /// or `Target { smoothing }` has `smoothing < 0`.
        pub fn with_categorical(
            mut self,
            col: impl Into<String>,
            categories: Vec<i64>,
            strategy: CategoricalStrategy,
        ) -> Result<Self> {
            if categories.is_empty() {
                return Err(ForecastError::InvalidParameter(
                    "with_categorical: categories must be non-empty".into(),
                ));
            }
            let mut sorted = categories.clone();
            sorted.sort_unstable();
            sorted.dedup();
            if sorted.len() != categories.len() {
                return Err(ForecastError::InvalidParameter(
                    "with_categorical: categories must be unique".into(),
                ));
            }
            if let CategoricalStrategy::Target { smoothing } = strategy {
                if !(smoothing.is_finite() && smoothing >= 0.0) {
                    return Err(ForecastError::InvalidParameter(format!(
                        "with_categorical: Target smoothing must be ≥ 0 and finite, got {}",
                        smoothing
                    )));
                }
            }
            self.exog_features.push(ExogFeatureSpec::Categorical {
                col: col.into(),
                categories,
                strategy,
            });
            Ok(self)
        }

        /// Apply regular differencing of order `d` before fitting.
        ///
        /// The model automatically integrates (undoes differencing) during predict.
        pub fn differencing(mut self, d: usize) -> Self {
            self.diff_order = d;
            self
        }

        /// Apply seasonal differencing of order `D` with the given period before fitting.
        ///
        /// Can be called multiple times for multi-seasonal series (e.g., weekly + yearly).
        /// Each call adds a differencing step; they are applied in order during fit
        /// and integrated in reverse order during predict.
        pub fn seasonal_differencing(mut self, d: usize, period: usize) -> Self {
            self.seasonal_diffs.push((d, period));
            self
        }

        /// Apply fractional differencing before fitting.
        ///
        /// Differences the series by order `d` (typically 0 < d < 1) using the
        /// binomial series expansion `(1-B)^d`. This removes just enough memory
        /// to achieve stationarity while preserving predictive signal.
        ///
        /// Note: fractional differencing is NOT exactly invertible, so forecasts
        /// remain on the differenced scale. Use this for feature engineering in
        /// regression models, not for direct point forecasting that needs
        /// original-scale output.
        ///
        /// Reference: Lopez de Prado, *Advances in Financial Machine Learning* (2018).
        pub fn fractional_differencing(mut self, d: f64) -> Self {
            self.frac_diff_order = Some(d);
            self
        }

        /// Maximum warmup required by any recursive feature (0 if none).
        fn max_recursive_warmup(&self) -> usize {
            self.recursive_features
                .iter()
                .map(|f| f.warmup())
                .max()
                .unwrap_or(0)
        }

        /// Maximum warmup required by any derived exog feature.
        fn max_exog_warmup(&self) -> usize {
            self.exog_features
                .iter()
                .map(|s| s.warmup())
                .max()
                .unwrap_or(0)
        }

        /// Number of observations lost to lagging and recursive-feature warmup.
        ///
        /// The design matrix drops the first `lag_offset()` rows.
        fn lag_offset(&self) -> usize {
            self.max_effective_lag()
                .max(self.max_recursive_warmup())
                .max(self.max_exog_warmup())
        }

        /// Build feature column names for a given TimeSeries.
        pub fn feature_names(&self, exog_names: &[String]) -> Vec<String> {
            let mut names = Vec::new();
            if self.use_trend {
                names.push("__trend".to_string());
            }
            for lag in self.effective_lags() {
                names.push(format!("__lag_{}", lag));
            }
            // Trend component columns
            for trend in &self.trend_components {
                match trend {
                    TrendType::Linear => names.push("__linear_trend".to_string()),
                    TrendType::Quadratic => names.push("__quadratic_trend".to_string()),
                    TrendType::Cubic => names.push("__cubic_trend".to_string()),
                    TrendType::Exponential => names.push("__exp_trend".to_string()),
                    TrendType::TheilSen => names.push("__theilsen_trend".to_string()),
                }
            }
            // Seasonal component columns
            for seasonal in &self.seasonal_components {
                match seasonal {
                    SeasonalSpec::Fourier { period, order } => {
                        for k in 1..=*order {
                            names.push(format!("__fourier_p{}_sin_{}", period, k));
                            names.push(format!("__fourier_p{}_cos_{}", period, k));
                        }
                    }
                    SeasonalSpec::Dummy(period) => {
                        names.push(format!("__seasonal_{}", period));
                    }
                }
            }
            // Structural feature columns
            for sf in &self.structural_features {
                names.extend(sf.column_names());
            }
            // Recursive feature columns (rolling stats, etc.)
            for rf in &self.recursive_features {
                names.extend(rf.column_names());
            }
            if self.use_exog {
                for name in exog_names {
                    names.push(name.clone());
                }
            }
            for spec in &self.exog_features {
                names.extend(spec.column_names());
            }
            names
        }

        /// Classify all feature columns by their data-leakage risk.
        ///
        /// Returns `(name, safety)` pairs in design-matrix column order.
        pub fn classify_features(&self, exog_names: &[String]) -> Vec<(String, FeatureSafety)> {
            let mut result = Vec::new();
            if self.use_trend {
                result.push(("__trend".to_string(), FeatureSafety::Deterministic));
            }
            for lag in self.effective_lags() {
                result.push((format!("__lag_{}", lag), FeatureSafety::Deterministic));
            }
            for trend in &self.trend_components {
                let name = match trend {
                    TrendType::Linear => "__linear_trend",
                    TrendType::Quadratic => "__quadratic_trend",
                    TrendType::Cubic => "__cubic_trend",
                    TrendType::Exponential => "__exp_trend",
                    TrendType::TheilSen => "__theilsen_trend",
                };
                result.push((name.to_string(), trend.safety()));
            }
            for seasonal in &self.seasonal_components {
                match seasonal {
                    SeasonalSpec::Fourier { period, order } => {
                        for k in 1..=*order {
                            result.push((
                                format!("__fourier_p{}_sin_{}", period, k),
                                FeatureSafety::Deterministic,
                            ));
                            result.push((
                                format!("__fourier_p{}_cos_{}", period, k),
                                FeatureSafety::Deterministic,
                            ));
                        }
                    }
                    SeasonalSpec::Dummy(period) => {
                        result.push((
                            format!("__seasonal_{}", period),
                            FeatureSafety::DataDependent,
                        ));
                    }
                }
            }
            for sf in &self.structural_features {
                for col_name in sf.column_names() {
                    result.push((col_name, FeatureSafety::Structural));
                }
            }
            for rf in &self.recursive_features {
                for col_name in rf.column_names() {
                    // Recursive features (rolling stats, EWMs) are derived
                    // deterministically from past target values at predict
                    // time — no fitted parameters, no external data.
                    result.push((col_name, FeatureSafety::Deterministic));
                }
            }
            if self.use_exog {
                for name in exog_names {
                    result.push((name.clone(), FeatureSafety::External));
                }
            }
            for spec in &self.exog_features {
                let safety = match spec {
                    ExogFeatureSpec::Categorical { strategy, .. } => match strategy {
                        // Target / Count encodings depend on training data
                        // (y or category frequencies) — must refit on each CV fold.
                        CategoricalStrategy::Target { .. } | CategoricalStrategy::Count => {
                            FeatureSafety::DataDependent
                        }
                        CategoricalStrategy::OneHot { .. } | CategoricalStrategy::Ordinal => {
                            FeatureSafety::External
                        }
                    },
                    _ => FeatureSafety::External,
                };
                for col_name in spec.column_names() {
                    result.push((col_name, safety));
                }
            }
            result
        }

        /// Build the design matrix and target vector from a TimeSeries.
        ///
        /// Returns `(X, y, n_train, exog_names, fitted_components)` where `n_train`
        /// is the number of usable rows (= n - max_lag) and `fitted_components`
        /// holds any trend/seasonal components fitted during matrix construction.
        fn build_matrices(
            &self,
            series: &TimeSeries,
        ) -> Result<(
            Mat<f64>,
            Col<f64>,
            usize,
            Vec<String>,
            Vec<FittedComponentState>,
            HashMap<usize, FittedCategoricalEncoder>,
        )> {
            let values = series.primary_values();
            let n = values.len();
            let offset = self.lag_offset();

            if n <= offset {
                return Err(ForecastError::InsufficientData {
                    needed: offset + 2,
                    got: n,
                    hint: Some(format!(
                        "need > {} observations for lags {:?}",
                        offset,
                        self.effective_lags()
                    )),
                });
            }

            let n_train = n - offset;

            // Collect exogenous regressor names (sorted for determinism).
            // Columns already claimed by an exog_features spec
            // (Categorical / Lag / Rolling / Polynomial / Interaction)
            // are excluded from the raw exog block — they're already
            // represented downstream by the spec's derived columns, so
            // including them as raw would cause collinearity (clearest
            // case: Categorical OneHot + raw == singular). Closes #113.
            let exog_names = if self.use_exog && series.has_regressors() {
                let claimed: std::collections::HashSet<&str> = self
                    .exog_features
                    .iter()
                    .flat_map(|s| s.referenced_cols())
                    .collect();
                let mut names: Vec<String> = series
                    .all_regressors()
                    .keys()
                    .filter(|k| !claimed.contains(k.as_str()))
                    .cloned()
                    .collect();
                names.sort();
                names
            } else {
                Vec::new()
            };

            let feature_names = self.feature_names(&exog_names);
            let n_features = feature_names.len();

            if n_features == 0 {
                return Err(ForecastError::InvalidParameter(
                    "No features configured — enable at least one of: trend, lags, components, or exog"
                        .to_string(),
                ));
            }

            // ── Fit trend/seasonal components ───────────────────────────
            let mut fitted_components = Vec::new();

            for trend_type in &self.trend_components {
                let comp = match trend_type {
                    TrendType::Linear => {
                        let mut p = PolynomialTrend::new(1).with_recency(Recency::Full);
                        p.fit_trend(values)?;
                        FittedComponentState::Polynomial(p)
                    }
                    TrendType::Quadratic => {
                        let mut p = PolynomialTrend::new(2).with_recency(Recency::Full);
                        p.fit_trend(values)?;
                        FittedComponentState::Polynomial(p)
                    }
                    TrendType::Cubic => {
                        let mut p = PolynomialTrend::new(3).with_recency(Recency::Full);
                        p.fit_trend(values)?;
                        FittedComponentState::Polynomial(p)
                    }
                    TrendType::Exponential => {
                        let mut e = ExponentialTrend::new().with_recency(Recency::Full);
                        e.fit_trend(values)?;
                        FittedComponentState::Exponential(e)
                    }
                    TrendType::TheilSen => {
                        let mut t = TheilSenTrend::new().with_recency(Recency::Full);
                        t.fit_trend(values)?;
                        FittedComponentState::TheilSen(t)
                    }
                };
                fitted_components.push(comp);
            }

            for seasonal in &self.seasonal_components {
                match seasonal {
                    SeasonalSpec::Fourier { period, order } => {
                        fitted_components.push(FittedComponentState::Fourier {
                            period: *period,
                            order: *order,
                        });
                    }
                    SeasonalSpec::Dummy(period) => {
                        let mut d = DummySeasonality::new();
                        d.fit_seasonal(values, *period)?;
                        fitted_components.push(FittedComponentState::Dummy(d));
                    }
                }
            }

            // ── Build design matrix ─────────────────────────────────────
            let mut x = Mat::zeros(n_train, n_features);
            let mut y = Col::zeros(n_train);

            // Populate target
            for i in 0..n_train {
                y[i] = values[offset + i];
            }

            // Populate features
            let mut col_idx = 0;

            // Trend: index of the observation (relative to full series)
            if self.use_trend {
                for i in 0..n_train {
                    x[(i, col_idx)] = (offset + i) as f64;
                }
                col_idx += 1;
            }

            // Lags: y[t-k] for each specified lag k
            for lag in self.effective_lags() {
                for i in 0..n_train {
                    x[(i, col_idx)] = values[offset + i - lag];
                }
                col_idx += 1;
            }

            // Trend/seasonal component columns
            for comp in &fitted_components {
                match comp {
                    FittedComponentState::Polynomial(p) => {
                        let fitted = p.fitted_trend();
                        for i in 0..n_train {
                            x[(i, col_idx)] = fitted[offset + i];
                        }
                        col_idx += 1;
                    }
                    FittedComponentState::Exponential(e) => {
                        let fitted = e.fitted_trend();
                        for i in 0..n_train {
                            x[(i, col_idx)] = fitted[offset + i];
                        }
                        col_idx += 1;
                    }
                    FittedComponentState::TheilSen(t) => {
                        let fitted = t.fitted_trend();
                        for i in 0..n_train {
                            x[(i, col_idx)] = fitted[offset + i];
                        }
                        col_idx += 1;
                    }
                    FittedComponentState::Dummy(d) => {
                        let fitted = d.fitted_seasonal();
                        for i in 0..n_train {
                            x[(i, col_idx)] = fitted[offset + i];
                        }
                        col_idx += 1;
                    }
                    FittedComponentState::Fourier { period, order } => {
                        let timestamps: Vec<f64> = (0..n).map(|i| i as f64).collect();
                        let basis = fourier_terms(&timestamps, *period as f64, *order)?;
                        for basis_vec in &basis {
                            for i in 0..n_train {
                                x[(i, col_idx)] = basis_vec[offset + i];
                            }
                            col_idx += 1;
                        }
                    }
                    FittedComponentState::Structural { .. } => {
                        // Handled below via structural_features iteration
                    }
                }
            }

            // Structural feature columns
            for sf in &self.structural_features {
                let columns = sf.compute(n);
                let fill = sf.fill_values(n);
                for col_vals in &columns {
                    for i in 0..n_train {
                        x[(i, col_idx)] = col_vals[offset + i];
                    }
                    col_idx += 1;
                }
                fitted_components.push(FittedComponentState::Structural { fill_values: fill });
            }

            // Recursive feature columns (rolling statistics, EWMs, …).
            // For row i, the "target" is values[offset + i] and the feature
            // must be computed from strictly-prior values only.
            for rf in &self.recursive_features {
                let n_cols = rf.n_columns();
                let mut scratch = vec![0.0_f64; n_cols];
                for i in 0..n_train {
                    rf.compute_fit(values, offset + i, &mut scratch);
                    for (k, &v) in scratch.iter().enumerate() {
                        x[(i, col_idx + k)] = v;
                    }
                }
                col_idx += n_cols;
            }

            // Exogenous regressors (sliced to match after lag offset)
            if self.use_exog {
                let regressors = series.all_regressors();
                for name in &exog_names {
                    if let Some(reg_values) = regressors.get(name) {
                        for i in 0..n_train {
                            let idx = offset + i;
                            if idx < reg_values.len() {
                                x[(i, col_idx)] = reg_values[idx];
                            }
                        }
                    }
                    col_idx += 1;
                }
            }

            // Derived exog features (lags, rolling, polynomial, interaction, categorical)
            let mut categorical_encoders: HashMap<usize, FittedCategoricalEncoder> = HashMap::new();
            if !self.exog_features.is_empty() {
                let regressors = series.all_regressors();
                let lookup_col = |name: &str| -> Result<&Vec<f64>> {
                    regressors.get(name).ok_or_else(|| {
                        ForecastError::InvalidParameter(format!(
                            "exog_features: column '{}' not found in series regressors",
                            name
                        ))
                    })
                };
                for (spec_idx, spec) in self.exog_features.iter().enumerate() {
                    match spec {
                        ExogFeatureSpec::Lag { col, lags } => {
                            let reg_values = lookup_col(col)?;
                            for &k in lags {
                                for i in 0..n_train {
                                    let src_idx = offset + i - k;
                                    x[(i, col_idx)] = reg_values[src_idx];
                                }
                                col_idx += 1;
                            }
                        }
                        ExogFeatureSpec::Rolling {
                            col,
                            window,
                            lag,
                            kind,
                        } => {
                            let reg_values = lookup_col(col)?;
                            for i in 0..n_train {
                                let end_excl = offset + i + 1 - lag;
                                let start = end_excl - window;
                                x[(i, col_idx)] = kind.compute(&reg_values[start..end_excl]);
                            }
                            col_idx += 1;
                        }
                        ExogFeatureSpec::Polynomial { col, degree } => {
                            let reg_values = lookup_col(col)?;
                            for k in 2..=*degree {
                                for i in 0..n_train {
                                    x[(i, col_idx)] = reg_values[offset + i].powi(k as i32);
                                }
                                col_idx += 1;
                            }
                        }
                        ExogFeatureSpec::Interaction { col_a, col_b } => {
                            let va = lookup_col(col_a)?;
                            let vb = lookup_col(col_b)?;
                            for i in 0..n_train {
                                x[(i, col_idx)] = va[offset + i] * vb[offset + i];
                            }
                            col_idx += 1;
                        }
                        ExogFeatureSpec::Categorical {
                            col,
                            categories,
                            strategy,
                        } => {
                            let reg_values = lookup_col(col)?;
                            // Pre-fit encoder for Count / Target strategies.
                            let encoder = match strategy {
                                CategoricalStrategy::OneHot { .. }
                                | CategoricalStrategy::Ordinal => None,
                                CategoricalStrategy::Count => {
                                    let mut counts: HashMap<i64, f64> = HashMap::new();
                                    for i in 0..n_train {
                                        let cat = reg_values[offset + i].round() as i64;
                                        *counts.entry(cat).or_insert(0.0) += 1.0;
                                    }
                                    Some(FittedCategoricalEncoder {
                                        map: counts,
                                        fallback: 0.0,
                                    })
                                }
                                CategoricalStrategy::Target { smoothing } => {
                                    // Per-category sums and counts, then smooth toward grand mean.
                                    let mut sums: HashMap<i64, f64> = HashMap::new();
                                    let mut counts: HashMap<i64, f64> = HashMap::new();
                                    let mut total = 0.0_f64;
                                    for i in 0..n_train {
                                        let cat = reg_values[offset + i].round() as i64;
                                        *sums.entry(cat).or_insert(0.0) += y[i];
                                        *counts.entry(cat).or_insert(0.0) += 1.0;
                                        total += y[i];
                                    }
                                    let grand_mean = if n_train > 0 {
                                        total / n_train as f64
                                    } else {
                                        0.0
                                    };
                                    let map: HashMap<i64, f64> = counts
                                        .iter()
                                        .map(|(cat, &n_c)| {
                                            let mean_c = sums[cat] / n_c;
                                            (
                                                *cat,
                                                (n_c * mean_c + smoothing * grand_mean)
                                                    / (n_c + smoothing),
                                            )
                                        })
                                        .collect();
                                    Some(FittedCategoricalEncoder {
                                        map,
                                        fallback: grand_mean,
                                    })
                                }
                            };
                            // Materialize columns for training rows.
                            match strategy {
                                CategoricalStrategy::OneHot { drop_first } => {
                                    let start = if *drop_first { 1 } else { 0 };
                                    for cat in &categories[start..] {
                                        for i in 0..n_train {
                                            let v = reg_values[offset + i].round() as i64;
                                            x[(i, col_idx)] = if v == *cat { 1.0 } else { 0.0 };
                                        }
                                        col_idx += 1;
                                    }
                                }
                                CategoricalStrategy::Ordinal => {
                                    let unknown = categories.len() as f64;
                                    for i in 0..n_train {
                                        let v = reg_values[offset + i].round() as i64;
                                        x[(i, col_idx)] = categories
                                            .iter()
                                            .position(|c| *c == v)
                                            .map(|p| p as f64)
                                            .unwrap_or(unknown);
                                    }
                                    col_idx += 1;
                                }
                                CategoricalStrategy::Count | CategoricalStrategy::Target { .. } => {
                                    let enc = encoder.as_ref().unwrap();
                                    for i in 0..n_train {
                                        let v = reg_values[offset + i].round() as i64;
                                        x[(i, col_idx)] = *enc.map.get(&v).unwrap_or(&enc.fallback);
                                    }
                                    col_idx += 1;
                                }
                            }
                            if let Some(enc) = encoder {
                                categorical_encoders.insert(spec_idx, enc);
                            }
                        }
                    }
                }
            }

            Ok((
                x,
                y,
                n_train,
                exog_names,
                fitted_components,
                categorical_encoders,
            ))
        }

        /// Build a design matrix for the forecast horizon.
        ///
        /// For lags: uses the last values from training + predicted values
        /// for multi-step recursive forecasting. Component columns are
        /// populated from the fitted components' `predict_trend` / `predict_seasonal`.
        fn build_future_matrix(
            &self,
            horizon: usize,
            n_total: usize,
            tail_values: &[f64],
            future_regressors: Option<&HashMap<String, Vec<f64>>>,
            exog_names: &[String],
            components: &[FittedComponentState],
            exog_tails: Option<&HashMap<String, Vec<f64>>>,
            categorical_encoders: Option<&HashMap<usize, FittedCategoricalEncoder>>,
        ) -> Result<Mat<f64>> {
            let feature_names = self.feature_names(exog_names);
            let n_features = feature_names.len();
            let mut x = Mat::zeros(horizon, n_features);

            let mut col_idx = 0;

            // Trend: continue the index
            if self.use_trend {
                for h in 0..horizon {
                    x[(h, col_idx)] = (n_total + h) as f64;
                }
                col_idx += 1;
            }

            // Lags: filled during recursive prediction (column indices stored)
            // Pre-fill from tail_values where possible
            for lag in self.effective_lags() {
                for h in 0..horizon {
                    if h >= lag {
                        // Will be filled recursively during prediction
                        x[(h, col_idx)] = f64::NAN; // placeholder
                    } else {
                        // Use known historical values
                        let idx = tail_values.len() as isize - lag as isize + h as isize;
                        if idx >= 0 {
                            x[(h, col_idx)] = tail_values[idx as usize];
                        }
                    }
                }
                col_idx += 1;
            }

            // Trend/seasonal component columns
            for comp in components {
                let future_cols = comp.predict(horizon, n_total);
                for col_vals in &future_cols {
                    for h in 0..horizon.min(col_vals.len()) {
                        x[(h, col_idx)] = col_vals[h];
                    }
                    col_idx += 1;
                }
            }

            // Recursive feature columns — filled per step during predict_recursive.
            // Leave placeholders for predict_recursive to overwrite.
            for rf in &self.recursive_features {
                col_idx += rf.n_columns();
            }

            // Exogenous regressors
            if self.use_exog {
                for name in exog_names {
                    if let Some(regs) = future_regressors {
                        if let Some(vals) = regs.get(name) {
                            for h in 0..horizon.min(vals.len()) {
                                x[(h, col_idx)] = vals[h];
                            }
                        }
                    }
                    col_idx += 1;
                }
            }

            // Derived exog features: resolve against exog_tails ++
            // future_regressors as a virtual full sequence (Lag, Rolling) or
            // pure future_regressors lookups (Polynomial, Interaction).
            if !self.exog_features.is_empty() {
                let warmup = self.max_exog_warmup();
                let virt_at = |col: &str, v_idx: isize| -> Option<f64> {
                    let tail = exog_tails
                        .and_then(|t| t.get(col))
                        .map(|v| v.as_slice())
                        .unwrap_or(&[]);
                    let future = future_regressors
                        .and_then(|m| m.get(col))
                        .map(|v| v.as_slice())
                        .unwrap_or(&[]);
                    let tlen = tail.len() as isize;
                    if v_idx < 0 {
                        None
                    } else if v_idx < tlen {
                        Some(tail[v_idx as usize])
                    } else {
                        let f_idx = (v_idx - tlen) as usize;
                        future.get(f_idx).copied()
                    }
                };
                let future_at = |col: &str, h: usize| -> Option<f64> {
                    future_regressors
                        .and_then(|m| m.get(col))
                        .and_then(|v| v.get(h).copied())
                };
                for (spec_idx, spec) in self.exog_features.iter().enumerate() {
                    match spec {
                        ExogFeatureSpec::Lag { col, lags } => {
                            for &k in lags {
                                for h in 0..horizon {
                                    // tail covers [n_train - warmup .. n_train),
                                    // so virtual_idx = warmup + h - k.
                                    let v_idx = warmup as isize + h as isize - k as isize;
                                    if let Some(v) = virt_at(col, v_idx) {
                                        x[(h, col_idx)] = v;
                                    }
                                }
                                col_idx += 1;
                            }
                        }
                        ExogFeatureSpec::Rolling {
                            col,
                            window,
                            lag,
                            kind,
                        } => {
                            for h in 0..horizon {
                                let end_v = warmup as isize + h as isize - *lag as isize;
                                let start_v = end_v - *window as isize + 1;
                                let mut buf = Vec::with_capacity(*window);
                                for v_idx in start_v..=end_v {
                                    if let Some(v) = virt_at(col, v_idx) {
                                        buf.push(v);
                                    }
                                }
                                if buf.len() == *window {
                                    x[(h, col_idx)] = kind.compute(&buf);
                                }
                            }
                            col_idx += 1;
                        }
                        ExogFeatureSpec::Polynomial { col, degree } => {
                            for k in 2..=*degree {
                                for h in 0..horizon {
                                    if let Some(v) = future_at(col, h) {
                                        x[(h, col_idx)] = v.powi(k as i32);
                                    }
                                }
                                col_idx += 1;
                            }
                        }
                        ExogFeatureSpec::Interaction { col_a, col_b } => {
                            for h in 0..horizon {
                                if let (Some(va), Some(vb)) =
                                    (future_at(col_a, h), future_at(col_b, h))
                                {
                                    x[(h, col_idx)] = va * vb;
                                }
                            }
                            col_idx += 1;
                        }
                        ExogFeatureSpec::Categorical {
                            col,
                            categories,
                            strategy,
                        } => match strategy {
                            CategoricalStrategy::OneHot { drop_first } => {
                                let start = if *drop_first { 1 } else { 0 };
                                for cat in &categories[start..] {
                                    for h in 0..horizon {
                                        if let Some(v) = future_at(col, h) {
                                            let cv = v.round() as i64;
                                            x[(h, col_idx)] = if cv == *cat { 1.0 } else { 0.0 };
                                        }
                                    }
                                    col_idx += 1;
                                }
                            }
                            CategoricalStrategy::Ordinal => {
                                let unknown = categories.len() as f64;
                                for h in 0..horizon {
                                    if let Some(v) = future_at(col, h) {
                                        let cv = v.round() as i64;
                                        x[(h, col_idx)] = categories
                                            .iter()
                                            .position(|c| *c == cv)
                                            .map(|p| p as f64)
                                            .unwrap_or(unknown);
                                    }
                                }
                                col_idx += 1;
                            }
                            CategoricalStrategy::Count | CategoricalStrategy::Target { .. } => {
                                let enc = categorical_encoders.and_then(|m| m.get(&spec_idx));
                                for h in 0..horizon {
                                    if let Some(v) = future_at(col, h) {
                                        let cv = v.round() as i64;
                                        x[(h, col_idx)] = enc
                                            .map(|e| *e.map.get(&cv).unwrap_or(&e.fallback))
                                            .unwrap_or(0.0);
                                    }
                                }
                                col_idx += 1;
                            }
                        },
                    }
                }
            }

            Ok(x)
        }
    }

    // ── Custom WLS-Logistic-Ridge fitted regressor ──────────────────

    /// FittedRegressor wrapper for the WLS-Logistic-Ridge backend.
    ///
    /// Holds a manually-constructed [`RegressionResult`] with the
    /// weighted-Ridge slopes (from OLS on the centred + sqrt-weighted +
    /// Tikhonov-augmented design) and the weighted-mean intercept
    /// (`ȳ_w − β · X̄_w`). Implements [`FittedRegressor::predict`] as the
    /// standard `intercept + Xβ` on the original (unweighted) scale.
    #[derive(Debug)]
    struct WlsLogisticRidgeFitted {
        result: anofox_regression::core::RegressionResult,
    }

    impl FittedRegressor for WlsLogisticRidgeFitted {
        fn predict(&self, x: &Mat<f64>) -> Col<f64> {
            let n = x.nrows();
            let p = x.ncols();
            let mut out = Col::zeros(n);
            let intercept = self.result.intercept.unwrap_or(0.0);
            for i in 0..n {
                let mut pred = intercept;
                for j in 0..p {
                    pred += self.result.coefficients[j] * x[(i, j)];
                }
                out[i] = pred;
            }
            out
        }

        fn result(&self) -> &anofox_regression::core::RegressionResult {
            &self.result
        }

        fn predict_with_interval(
            &self,
            x: &Mat<f64>,
            _interval: Option<IntervalType>,
            _level: f64,
        ) -> anofox_regression::core::PredictionResult {
            // Weighted-Ridge inference intervals are not standard
            // (sandwich / nominal forms vary by literature) so we return
            // point predictions with zero-width bounds — downstream uses
            // residual-based conformal intervals when needed.
            let fit = self.predict(x);
            let n = fit.nrows();
            anofox_regression::core::PredictionResult {
                fit: fit.clone(),
                lower: fit.clone(),
                upper: fit,
                se: Col::zeros(n),
            }
        }
    }

    // ── Fitted state ────────────────────────────────────────────────

    /// Internal state stored after fitting.
    struct FittedState {
        /// The fitted regression model (any backend).
        model: Box<dyn FittedRegressor + Send>,
        /// Feature configuration used.
        features: RegressionFeatures,
        /// Number of observations in the full series (before differencing).
        n_total: usize,
        /// Last `max_lag` values for recursive prediction (from differenced series).
        tail_values: Vec<f64>,
        /// In-sample fitted values (full length, NaN-padded for lags/differencing).
        fitted_values: Vec<f64>,
        /// In-sample residuals (full length, NaN-padded for lags/differencing).
        residuals: Vec<f64>,
        /// Exogenous regressor names (sorted).
        exog_names: Vec<String>,
        /// Tail of exogenous regressor values (for resolving lag/rolling
        /// derived-exog columns at predict time). Keyed by column name; each
        /// value has length `max_exog_warmup`. Empty when no derived-exog
        /// features are configured.
        exog_tails: HashMap<String, Vec<f64>>,
        /// Fitted categorical encoders (Count/Target strategies). Keyed by
        /// spec index in `features.exog_features`. Empty for OneHot/Ordinal.
        categorical_encoders: HashMap<usize, FittedCategoricalEncoder>,
        /// Fitted trend/seasonal components for generating future feature columns.
        components: Vec<FittedComponentState>,
        /// Original series values (stored when differencing is used, for integration).
        original_values: Option<Vec<f64>>,
    }

    impl std::fmt::Debug for FittedState {
        fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
            f.debug_struct("FittedState")
                .field("n_total", &self.n_total)
                .field(
                    "n_features",
                    &self.features.feature_names(&self.exog_names).len(),
                )
                .field("r_squared", &self.model.r_squared())
                .finish()
        }
    }

    // ── RegressionForecaster ────────────────────────────────────────

    /// A forecasting model backed by an external regression estimator.
    ///
    /// Wraps regression estimators from `anofox-regression` behind the
    /// [`Forecaster`] trait, enabling them to participate in pipelines,
    /// registries, ensembles, and cross-validation.
    ///
    /// # Backends
    ///
    /// | Backend | Use case |
    /// |---|---|
    /// | [`Ols`](RegressionBackend::Ols) | Default — unregularized |
    /// | [`Ridge`](RegressionBackend::Ridge) | L2 regularization — many features |
    /// | [`ElasticNet`](RegressionBackend::ElasticNet) | L1+L2 — feature selection |
    /// | [`Quantile`](RegressionBackend::Quantile) | Conditional quantile estimation |
    /// | [`Wls`](RegressionBackend::Wls) | Observation weighting / recency |
    /// | [`Rls`](RegressionBackend::Rls) | Adaptive / online coefficients |
    /// | [`Tweedie`](RegressionBackend::Tweedie) | GLM for count / continuous data |
    /// | [`Poisson`](RegressionBackend::Poisson) | GLM for count data |
    /// | [`Bls`](RegressionBackend::Bls) | Box-constrained coefficients |
    ///
    /// # Example
    ///
    /// ```rust,ignore
    /// use anofox_forecast::models::regression::{RegressionForecaster, RegressionFeatures};
    ///
    /// // OLS with trend + 3 lags
    /// let mut model = RegressionForecaster::ols(
    ///     RegressionFeatures::new().trend().lags(3),
    /// );
    ///
    /// // Ridge with Fourier seasonality
    /// let mut model = RegressionForecaster::ridge(
    ///     0.1,
    ///     RegressionFeatures::new().trend().fourier(7, 3),
    /// );
    /// ```
    #[derive(Debug)]
    pub struct RegressionForecaster {
        features: RegressionFeatures,
        backend: RegressionBackend,
        state: Option<FittedState>,
    }

    impl RegressionForecaster {
        /// Create a regression forecaster with the given backend and features.
        pub fn new(backend: RegressionBackend, features: RegressionFeatures) -> Self {
            Self {
                features,
                backend,
                state: None,
            }
        }

        /// Create a regression forecaster using OLS with the given features.
        pub fn ols(features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Ols, features)
        }

        /// Create a Ridge regression forecaster (L2 regularization).
        pub fn ridge(lambda: f64, features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Ridge { lambda }, features)
        }

        /// Create an Elastic Net forecaster (L1 + L2 regularization).
        ///
        /// `alpha` controls the L1/L2 mix: 0 = pure Ridge, 1 = pure Lasso.
        pub fn elastic_net(lambda: f64, alpha: f64, features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::ElasticNet { lambda, alpha }, features)
        }

        /// Create a quantile regression forecaster.
        ///
        /// `tau` is the quantile to estimate (0.5 = median).
        pub fn quantile(tau: f64, features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Quantile { tau }, features)
        }

        /// Create a WLS forecaster with exponential decay weighting.
        ///
        /// `decay` controls recency emphasis (e.g., 0.95 = moderate, 0.99 = mild).
        pub fn wls_decay(decay: f64, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Wls {
                    strategy: WeightStrategy::ExponentialDecay(decay),
                },
                features,
            )
        }

        /// Create a WLS forecaster with custom weights.
        pub fn wls(weights: Vec<f64>, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Wls {
                    strategy: WeightStrategy::Custom(weights),
                },
                features,
            )
        }

        /// Create a WLS forecaster with logistic (sigmoid) recency weights
        /// centred `offset` steps from the end of the training series.
        ///
        /// See [`WeightStrategy::Logistic`] for the formula. Typical
        /// `offset` values are 6–24 — e.g. `wls_logistic(12.0, …)` makes
        /// the last ~12 observations dominate the fit on a long series.
        pub fn wls_logistic(offset: f64, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Wls {
                    strategy: WeightStrategy::Logistic { offset },
                },
                features,
            )
        }

        /// Create a Recursive Least Squares forecaster (adaptive coefficients).
        ///
        /// `forgetting_factor` controls adaptation speed (0 < λ ≤ 1).
        /// 1.0 = equal weights, 0.95 = moderate adaptation.
        pub fn rls(forgetting_factor: f64, features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Rls { forgetting_factor }, features)
        }

        /// Create a forecaster that combines logistic recency weighting with
        /// L2 (Ridge) regularization in a single fit.
        ///
        /// Solves `β = (Xᵀ W X + λI)⁻¹ Xᵀ W y` where `W` is the diagonal
        /// of logistic recency weights (see [`WeightStrategy::Logistic`])
        /// and `λI` is the L2 penalty (intercept *is* penalised, matching
        /// standard Ridge convention).
        ///
        /// Useful when recency weighting reduces the effective `N` enough
        /// that an unregularised WLS fit becomes near-singular under a
        /// heavy feature design — see issue #103.
        ///
        /// Typical usage: `offset` in 6–24, `lambda` in 0.1–10.
        pub fn wls_logistic_ridge(offset: f64, lambda: f64, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::WlsLogisticRidge { offset, lambda },
                features,
            )
        }

        /// Create a Tweedie GLM forecaster.
        ///
        /// `var_power`: 0 = Gaussian, 1 = Poisson, 2 = Gamma, 3 = Inv-Gaussian.
        pub fn tweedie(var_power: f64, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Tweedie {
                    var_power,
                    link_power: None,
                },
                features,
            )
        }

        /// Create a Poisson GLM forecaster (for count data).
        pub fn poisson(features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Poisson, features)
        }

        /// Create a Bounded Least Squares forecaster (box-constrained coefficients).
        pub fn bls(lower: Option<f64>, upper: Option<f64>, features: RegressionFeatures) -> Self {
            Self::new(RegressionBackend::Bls { lower, upper }, features)
        }

        /// Create a Non-Negative Least Squares forecaster (all coefficients ≥ 0).
        pub fn nnls(features: RegressionFeatures) -> Self {
            Self::bls(Some(0.0), None, features)
        }

        /// Create a Dynamic Linear Model forecaster (time-varying parameters).
        ///
        /// Automatically generates candidate models from variable subsets and
        /// computes observation-level IC weights. Use `lowess_span` to smooth
        /// the weights over time (e.g., 0.3), or `None` for no smoothing.
        pub fn dynamic(features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Dynamic {
                    ic: InformationCriterion::AICc,
                    lowess_span: None,
                },
                features,
            )
        }

        /// Create a Dynamic Linear Model with LOWESS-smoothed weights.
        pub fn dynamic_smoothed(lowess_span: f64, features: RegressionFeatures) -> Self {
            Self::new(
                RegressionBackend::Dynamic {
                    ic: InformationCriterion::AICc,
                    lowess_span: Some(lowess_span),
                },
                features,
            )
        }

        // ── OLS convenience constructors (backward-compatible) ──────

        /// Create a trend-only OLS forecaster (linear regression on time index).
        pub fn linear_trend() -> Self {
            Self::ols(RegressionFeatures::new().trend().no_exog())
        }

        /// Create an autoregressive OLS forecaster with the given number of lags.
        pub fn ar(lags: usize) -> Self {
            Self::ols(RegressionFeatures::new().no_trend().lags(lags).no_exog())
        }

        /// Create a trend + autoregressive OLS forecaster.
        pub fn trend_ar(lags: usize) -> Self {
            Self::ols(RegressionFeatures::new().trend().lags(lags))
        }

        /// Create a trend + Fourier seasonality OLS forecaster.
        pub fn trend_fourier(period: usize, order: usize) -> Self {
            Self::ols(
                RegressionFeatures::new()
                    .trend()
                    .fourier(period, order)
                    .no_exog(),
            )
        }

        /// Get the feature configuration.
        pub fn features(&self) -> &RegressionFeatures {
            &self.features
        }

        /// Get the backend configuration.
        pub fn backend(&self) -> &RegressionBackend {
            &self.backend
        }

        /// Get the regression result (coefficients, R², etc.) if fitted.
        pub fn fitted_result(&self) -> Option<&anofox_regression::core::RegressionResult> {
            self.state.as_ref().map(|s| s.model.result())
        }

        /// Get R² of the fitted model.
        pub fn r_squared(&self) -> Option<f64> {
            self.state.as_ref().map(|s| s.model.r_squared())
        }

        /// Apply configured differencing to a series.
        fn apply_differencing(&self, values: &[f64]) -> Vec<f64> {
            use crate::models::arima::{difference, fractional_difference, seasonal_difference};

            let mut result = values.to_vec();

            // Seasonal differencing first (standard ARIMA convention), in order
            for &(d, period) in &self.features.seasonal_diffs {
                result = seasonal_difference(&result, d, period);
            }

            // Then regular differencing
            if self.features.diff_order > 0 {
                result = difference(&result, self.features.diff_order);
            }

            // Then fractional differencing (applied last)
            if let Some(d) = self.features.frac_diff_order {
                if d > 0.0 {
                    result = fractional_difference(&result, d, 1e-4);
                }
            }

            result
        }

        /// Integrate (undo differencing) forecast values back to original scale.
        fn apply_integration(&self, state: &FittedState, predictions: &[f64]) -> Vec<f64> {
            let original = match &state.original_values {
                Some(v) => v,
                None => return predictions.to_vec(),
            };

            use crate::models::arima::{integrate, seasonal_integrate};

            let mut result = predictions.to_vec();

            // Undo regular differencing first (reverse of application order)
            if self.features.diff_order > 0 {
                // The reference for regular integration is the original after all seasonal diffs
                let mut reference = original.clone();
                for &(d, period) in &self.features.seasonal_diffs {
                    reference = crate::models::arima::seasonal_difference(&reference, d, period);
                }
                result = integrate(&result, &reference, self.features.diff_order);
            }

            // Then undo seasonal differencing in reverse order
            for &(d, period) in self.features.seasonal_diffs.iter().rev() {
                // Reference for each level is original with all prior seasonal diffs applied
                // For the last applied (first undone), reference is original with all-but-last
                // For simplicity, use original — seasonal_integrate handles the seed correctly
                result = seasonal_integrate(&result, original, d, period);
            }

            result
        }

        /// Recursive multi-step prediction for models with lag features.
        fn predict_recursive(
            &self,
            state: &FittedState,
            horizon: usize,
            future_regressors: Option<&HashMap<String, Vec<f64>>>,
        ) -> Result<Vec<f64>> {
            let mut x_future = state.features.build_future_matrix(
                horizon,
                state.n_total,
                &state.tail_values,
                future_regressors,
                &state.exog_names,
                &state.components,
                Some(&state.exog_tails),
                Some(&state.categorical_encoders),
            )?;

            let eff_lags = state.features.effective_lags();
            let has_recursive = !state.features.recursive_features.is_empty();

            if eff_lags.is_empty() && !has_recursive {
                // No lags and no recursive features — direct prediction.
                let preds = state.model.predict(&x_future);
                return Ok(preds.iter().copied().collect());
            }

            // Recursive: predict one step at a time, feeding predictions back.
            let trend_offset = if state.features.use_trend { 1 } else { 0 };
            // Column layout: [trend?] [lags] [components…] [recursive features] [exog]
            // — matches build_future_matrix/build_matrices.
            let component_cols: usize = state.components.iter().map(|c| c.n_columns()).sum();
            let recursive_start_col = trend_offset + eff_lags.len() + component_cols;

            let mut predictions = Vec::with_capacity(horizon);
            let mut recent: Vec<f64> = state.tail_values.clone();

            for h in 0..horizon {
                // Update lag columns with most recent known/predicted values.
                for (col_offset, &lag) in eff_lags.iter().enumerate() {
                    let col = trend_offset + col_offset;
                    let idx = recent.len() as isize - lag as isize;
                    if idx >= 0 {
                        x_future[(h, col)] = recent[idx as usize];
                    }
                }

                // Update recursive feature columns from the rolling buffer.
                let mut rcol = recursive_start_col;
                for rf in &state.features.recursive_features {
                    let n_cols = rf.n_columns();
                    let mut scratch = vec![0.0_f64; n_cols];
                    rf.compute_predict_at(&recent, h, state.n_total, &mut scratch);
                    for (k, &v) in scratch.iter().enumerate() {
                        x_future[(h, rcol + k)] = v;
                    }
                    rcol += n_cols;
                }

                // Predict this single step.
                let row = x_future.submatrix(h, 0, 1, x_future.ncols());
                let row_mat = Mat::from_fn(1, row.ncols(), |r, c| row[(r, c)]);
                let pred = state.model.predict(&row_mat);
                let y_hat = pred[0];
                predictions.push(y_hat);
                recent.push(y_hat);
            }

            Ok(predictions)
        }
    }

    impl Clone for RegressionForecaster {
        fn clone(&self) -> Self {
            // State is not Clone (Box<dyn FittedRegressor>), so we only clone config
            Self {
                features: self.features.clone(),
                backend: self.backend.clone(),
                state: None,
            }
        }
    }

    impl Forecaster for RegressionForecaster {
        fn fit(&mut self, series: &TimeSeries) -> Result<()> {
            validate_series_complete(series)?;

            // Resolve auto-lag selection before anything else
            self.features.resolve_auto_lags(series)?;

            let values = series.primary_values();
            let n_original = values.len();

            // Apply differencing if configured
            let uses_diff = self.features.diff_order > 0
                || !self.features.seasonal_diffs.is_empty()
                || self.features.frac_diff_order.is_some();
            let original_values = if uses_diff {
                Some(values.to_vec())
            } else {
                None
            };

            let working_values = self.apply_differencing(values);

            // Build a temporary TimeSeries from differenced values for matrix construction.
            // Slice the original series to preserve regressors (trimmed to differenced length).
            let fit_series = if uses_diff {
                let diff_offset = n_original - working_values.len();
                let trimmed_ts = series.slice(diff_offset, n_original)?;
                // Replace primary values with differenced values, keeping regressors
                TimeSeries::new(
                    trimmed_ts.timestamps().to_vec(),
                    vec![working_values.clone()],
                    crate::core::ValueLayout::Column,
                    trimmed_ts.labels().to_vec(),
                    trimmed_ts.metadata().clone(),
                    Vec::new(),
                    None,
                    None,
                    trimmed_ts.calendar().cloned(),
                )?
            } else {
                series.clone()
            };

            let n = fit_series.primary_values().len();
            let (x, y, n_train, exog_names, components, categorical_encoders) =
                self.features.build_matrices(&fit_series)?;

            // Fit via the configured backend
            let fitted = self.backend.fit_to(&x, &y).map_err(|e| {
                ForecastError::ComputationError(format!(
                    "{} fit failed: {}",
                    self.backend.name(),
                    e
                ))
            })?;

            // In-sample predictions (on differenced scale)
            let in_sample_preds = fitted.predict(&x);

            // Build full-length fitted values (NaN-padded for lag offset + differencing)
            let diff_offset = n_original - n;
            let lag_offset = self.features.lag_offset();
            let total_offset = diff_offset + lag_offset;
            let mut fitted_values = vec![f64::NAN; n_original];
            let mut residuals = vec![f64::NAN; n_original];
            let diff_values = fit_series.primary_values();
            for i in 0..n_train {
                fitted_values[total_offset + i] = in_sample_preds[i];
                residuals[total_offset + i] = diff_values[lag_offset + i] - in_sample_preds[i];
            }

            // Store tail values for recursive prediction (from differenced series).
            // Must cover the largest history reach: max lag AND largest
            // rolling window + lag (i.e. recursive warmup).
            let tail_len = self
                .features
                .max_effective_lag()
                .max(self.features.max_recursive_warmup())
                .max(1);
            let tail_values =
                working_values[working_values.len().saturating_sub(tail_len)..].to_vec();

            // Capture exog tails for derived-exog features (lag / rolling)
            // resolution at predict time. Each tail is the last
            // `max_exog_warmup` values of the named regressor.
            let mut exog_tails: HashMap<String, Vec<f64>> = HashMap::new();
            let exog_warmup = self.features.max_exog_warmup();
            if exog_warmup > 0 {
                let regressors = series.all_regressors();
                let mut seen: std::collections::HashSet<String> = Default::default();
                for spec in &self.features.exog_features {
                    for col in spec.referenced_cols() {
                        if !seen.insert(col.to_string()) {
                            continue;
                        }
                        if let Some(vals) = regressors.get(col) {
                            let start = vals.len().saturating_sub(exog_warmup);
                            exog_tails.insert(col.to_string(), vals[start..].to_vec());
                        }
                    }
                }
            }

            self.state = Some(FittedState {
                model: fitted,
                features: self.features.clone(),
                n_total: n,
                tail_values,
                fitted_values,
                residuals,
                exog_names,
                exog_tails,
                categorical_encoders,
                components,
                original_values,
            });

            Ok(())
        }

        fn predict(&self, horizon: usize) -> Result<Forecast> {
            let state = self
                .state
                .as_ref()
                .ok_or(ForecastError::FitRequired { model: None })?;

            if horizon == 0 {
                return Ok(Forecast::new());
            }

            // If model has exog and was fit with exog, require predict_with_exog
            if !state.exog_names.is_empty() {
                return Err(ForecastError::InvalidParameter(
                    "Model was fit with exogenous regressors; use predict_with_exog() \
                     to provide future regressor values"
                        .to_string(),
                ));
            }

            let predictions = self.predict_recursive(state, horizon, None)?;
            let predictions = self.apply_integration(state, &predictions);
            Ok(Forecast::from_values(predictions))
        }

        fn predict_with_intervals(&self, horizon: usize, level: f64) -> Result<Forecast> {
            let state = self
                .state
                .as_ref()
                .ok_or(ForecastError::FitRequired { model: None })?;

            if horizon == 0 {
                return Ok(Forecast::new());
            }

            // OLS prediction intervals are only valid for direct (non-recursive)
            // forecasts.  Recursive models (max_lag > 0) feed predicted values
            // back as features, so the standard interval formula does not apply;
            // fall back to point-only predictions in that case.
            if !state.features.effective_lags().is_empty() {
                return self.predict(horizon);
            }

            if !state.exog_names.is_empty() {
                return Err(ForecastError::InvalidParameter(
                    "Model was fit with exogenous regressors; use predict_with_exog() \
                     to provide future regressor values"
                        .to_string(),
                ));
            }

            let x_future = state.features.build_future_matrix(
                horizon,
                state.n_total,
                &state.tail_values,
                None,
                &state.exog_names,
                &state.components,
                Some(&state.exog_tails),
                Some(&state.categorical_encoders),
            )?;

            let pred_result =
                state
                    .model
                    .predict_with_interval(&x_future, Some(IntervalType::Prediction), level);

            let values: Vec<f64> = pred_result.fit.iter().copied().collect();
            let lower: Vec<f64> = pred_result.lower.iter().copied().collect();
            let upper: Vec<f64> = pred_result.upper.iter().copied().collect();

            // If intervals contain NaN (e.g. xtx_inverse unavailable), return
            // point-only to avoid misleading results.
            if lower.iter().any(|v| v.is_nan()) || upper.iter().any(|v| v.is_nan()) {
                return Ok(Forecast::from_values(values));
            }

            Ok(Forecast::from_values_with_intervals(values, lower, upper))
        }

        fn supports_exog(&self) -> bool {
            self.features.use_exog
        }

        fn has_exog(&self) -> bool {
            self.state
                .as_ref()
                .map(|s| !s.exog_names.is_empty())
                .unwrap_or(false)
        }

        fn exog_names(&self) -> Option<&[String]> {
            self.state
                .as_ref()
                .filter(|s| !s.exog_names.is_empty())
                .map(|s| s.exog_names.as_slice())
        }

        fn predict_with_exog(
            &self,
            horizon: usize,
            future_regressors: &HashMap<String, Vec<f64>>,
        ) -> Result<Forecast> {
            let state = self
                .state
                .as_ref()
                .ok_or(ForecastError::FitRequired { model: None })?;

            if horizon == 0 {
                return Ok(Forecast::new());
            }

            // Validate that all required regressors are provided —
            // both raw exog columns (state.exog_names) AND columns
            // claimed by exog_features specs (Categorical / Lag /
            // Rolling / Polynomial / Interaction). The two are
            // disjoint after the #113 fix; we union them here so a
            // claimed column being absent still surfaces an error
            // instead of silently producing zeros.
            let mut required: std::collections::BTreeSet<&str> =
                state.exog_names.iter().map(|s| s.as_str()).collect();
            for spec in &state.features.exog_features {
                for col in spec.referenced_cols() {
                    required.insert(col);
                }
            }
            for name in &required {
                match future_regressors.get(*name) {
                    None => {
                        return Err(ForecastError::InvalidParameter(format!(
                            "Missing future regressor '{}'. Required: {:?}",
                            name, required
                        )));
                    }
                    Some(vals) if vals.len() < horizon => {
                        return Err(ForecastError::DimensionMismatch {
                            expected: horizon,
                            got: vals.len(),
                        });
                    }
                    _ => {}
                }
            }

            let predictions = self.predict_recursive(state, horizon, Some(future_regressors))?;
            let predictions = self.apply_integration(state, &predictions);
            Ok(Forecast::from_values(predictions))
        }

        fn fitted_values(&self) -> Option<&[f64]> {
            self.state.as_ref().map(|s| s.fitted_values.as_slice())
        }

        fn residuals(&self) -> Option<&[f64]> {
            self.state.as_ref().map(|s| s.residuals.as_slice())
        }

        fn name(&self) -> &str {
            self.backend.name()
        }
    }

    impl crate::models::inspect::Inspectable for RegressionForecaster {
        fn explanation(&self) -> Result<crate::models::inspect::Explanation> {
            use crate::models::inspect::{Explanation, RegressionExplanation};

            let state = self.state.as_ref().ok_or(ForecastError::FitRequired {
                model: Some("RegressionForecaster".to_string()),
            })?;

            let result = state.model.result();
            let feature_names = state.features.feature_names(&state.exog_names);

            let n_coef = result.coefficients.nrows();
            let mut coefficients = Vec::with_capacity(n_coef);
            for i in 0..n_coef {
                coefficients.push(result.coefficients[i]);
            }

            // Standard errors are populated when the backend ran inference
            // (OLS / Ridge / ElasticNet / WLS by default). For non-linear
            // backends (Poisson / Tweedie / BLS / RLS) std_errors stays
            // None and we emit empty vec + NaN, per RegressionExplanation
            // doc.
            let coef_std_errors: Vec<f64> = match &result.std_errors {
                Some(se) => {
                    let n = se.nrows();
                    let mut out = Vec::with_capacity(n);
                    for i in 0..n {
                        out.push(se[i]);
                    }
                    out
                }
                None => Vec::new(),
            };
            let intercept_std_error = result.intercept_std_error.unwrap_or(f64::NAN);

            Ok(Explanation::Regression(RegressionExplanation {
                feature_names,
                coefficients,
                intercept: result.intercept.unwrap_or(0.0),
                r_squared: result.r_squared,
                backend: self.backend.name().to_string(),
                coef_std_errors,
                intercept_std_error,
                fitted_values: state.fitted_values.clone(),
                residuals: state.residuals.clone(),
            }))
        }
    }

    #[cfg(test)]
    mod tests {
        use super::*;
        use crate::core::{CalendarAnnotations, TimeSeriesBuilder};
        use approx::assert_relative_eq;
        use chrono::{Duration, TimeZone, Utc};

        fn make_timestamps(n: usize) -> Vec<chrono::DateTime<Utc>> {
            let start = Utc.with_ymd_and_hms(2024, 1, 1, 0, 0, 0).unwrap();
            (0..n).map(|i| start + Duration::days(i as i64)).collect()
        }

        fn make_linear_ts(n: usize) -> TimeSeries {
            // y = 2*t + 10 + small noise
            let values: Vec<f64> = (0..n)
                .map(|i| 2.0 * i as f64 + 10.0 + 0.01 * (i as f64 * 0.7).sin())
                .collect();
            TimeSeries::univariate(make_timestamps(n), values).unwrap()
        }

        #[test]
        fn inspectable_linear_trend_returns_regression_explanation() {
            use crate::models::inspect::{Explanation, Inspectable};

            let ts = make_linear_ts(40);
            let mut model = RegressionForecaster::linear_trend();
            assert!(
                matches!(model.explanation(), Err(ForecastError::FitRequired { .. })),
                "explanation() before fit must return FitRequired"
            );

            model.fit(&ts).unwrap();
            let Explanation::Regression(e) = model.explanation().unwrap() else {
                panic!("expected Explanation::Regression");
            };

            assert_eq!(e.backend, "OLS");
            assert_eq!(e.coefficients.len(), e.feature_names.len());
            assert!(e.feature_names.iter().any(|n| n == "__trend"));
            assert_eq!(e.fitted_values.len(), 40);
            assert_eq!(e.residuals.len(), 40);
            // Linear-trend OLS on a nearly-linear series should fit R² ≈ 1.
            assert!(e.r_squared > 0.99, "R² unexpectedly low: {}", e.r_squared);
            // The trend slope coefficient should recover ≈ 2.0.
            let trend_idx = e.feature_names.iter().position(|n| n == "__trend").unwrap();
            assert_relative_eq!(e.coefficients[trend_idx], 2.0, epsilon = 0.05);
        }

        #[test]
        fn inspectable_ols_populates_standard_errors() {
            use crate::models::inspect::{Explanation, Inspectable};

            let ts = make_linear_ts(40);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let Explanation::Regression(e) = model.explanation().unwrap() else {
                panic!("expected Explanation::Regression");
            };

            // OLS computes inference by default → SEs aligned with coefficients.
            assert_eq!(
                e.coef_std_errors.len(),
                e.coefficients.len(),
                "coef_std_errors must align with coefficients"
            );
            for (i, &se) in e.coef_std_errors.iter().enumerate() {
                assert!(se.is_finite(), "SE[{}] not finite: {}", i, se);
                assert!(se > 0.0, "SE[{}] not positive: {}", i, se);
            }
            assert!(
                e.intercept_std_error.is_finite(),
                "intercept SE not finite: {}",
                e.intercept_std_error
            );
            assert!(
                e.intercept_std_error > 0.0,
                "intercept SE not positive: {}",
                e.intercept_std_error
            );
        }

        #[test]
        fn ols_linear_trend_fit_predict() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Should continue the linear trend: y ≈ 2*t + 10
            for (h, &pred) in forecast.primary().iter().enumerate() {
                let expected = 2.0 * (50 + h) as f64 + 10.0;
                assert_relative_eq!(pred, expected, epsilon = 0.5);
            }
        }

        #[test]
        fn ols_linear_trend_fitted_values() {
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let fitted = model.fitted_values().unwrap();
            assert_eq!(fitted.len(), 30);

            // All should be finite (no lags = no NaN padding)
            for &v in fitted {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn ols_ar_model() {
            // AR(1) process: y[t] = 0.8 * y[t-1] + 1.0
            let n = 100;
            let mut values = vec![10.0];
            for i in 1..n {
                values.push(0.8 * values[i - 1] + 1.0 + 0.01 * (i as f64).sin());
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ar(1);
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Predictions should converge toward the stationary mean ≈ 5.0
            for &pred in forecast.primary() {
                assert!(pred.is_finite());
                assert!(pred > 0.0 && pred < 20.0);
            }
        }

        #[test]
        fn ols_ar_fitted_has_nan_padding() {
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::ar(3);
            model.fit(&ts).unwrap();

            let fitted = model.fitted_values().unwrap();
            assert_eq!(fitted.len(), 30);

            // First 3 values should be NaN (lag offset)
            assert!(fitted[0].is_nan());
            assert!(fitted[1].is_nan());
            assert!(fitted[2].is_nan());
            // Rest should be finite
            assert!(fitted[3].is_finite());
        }

        #[test]
        fn ols_trend_ar_combined() {
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::trend_ar(2);
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Trend+AR model should produce finite, increasing forecasts
            for &pred in forecast.primary() {
                assert!(pred.is_finite());
            }
            // Should be increasing (upward trend)
            assert!(forecast.primary()[4] > forecast.primary()[0]);
        }

        #[test]
        fn ols_with_exogenous_regressors() {
            // y = 3*x + 5 + trend
            let n = 50;
            let x_vals: Vec<f64> = (0..n).map(|i| (i as f64 * 0.3).sin()).collect();
            let values: Vec<f64> = (0..n)
                .map(|i| 3.0 * x_vals[i] + 5.0 + 0.1 * i as f64)
                .collect();

            let cal = CalendarAnnotations::new()
                .with_regressor("temperature".to_string(), x_vals.clone());

            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(values)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(RegressionFeatures::new().trend().no_exog());
            // First verify it works without exog
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Now with exog
            let mut model_exog = RegressionForecaster::ols(RegressionFeatures::new().trend());
            model_exog.fit(&ts).unwrap();

            assert!(model_exog.supports_exog());
            assert!(model_exog.has_exog());
            assert_eq!(model_exog.exog_names().unwrap(), &["temperature"]);

            // predict() should error because exog regressors are needed
            assert!(model_exog.predict(5).is_err());

            // predict_with_exog() should work
            let future_x: Vec<f64> = (n..n + 5).map(|i| (i as f64 * 0.3).sin()).collect();
            let mut future_regs = HashMap::new();
            future_regs.insert("temperature".to_string(), future_x);
            let forecast = model_exog.predict_with_exog(5, &future_regs).unwrap();
            assert_eq!(forecast.primary().len(), 5);
        }

        #[test]
        fn ols_with_exog_lags_recovers_lagged_relationship() {
            // y_t = 2 * x_{t-3} + noise — model with exog_lags(x, [3]) should fit cleanly.
            let n = 80;
            let x_vals: Vec<f64> = (0..n).map(|i| (i as f64 * 0.2).sin()).collect();
            let mut y = vec![0.0; n];
            for i in 3..n {
                let noise = ((i * 7 + 11) % 13) as f64 * 0.001;
                y[i] = 2.0 * x_vals[i - 3] + noise;
            }
            let cal = CalendarAnnotations::new().with_regressor("x".to_string(), x_vals.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_exog_lags("x", &[3]),
            );
            model.fit(&ts).unwrap();

            // Future exog: continue the same generator.
            let future_x: Vec<f64> = (n..n + 5).map(|i| (i as f64 * 0.2).sin()).collect();
            let mut future_regs = HashMap::new();
            future_regs.insert("x".to_string(), future_x.clone());
            let forecast = model.predict_with_exog(5, &future_regs).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Step h=0 should ≈ 2 * x_{n-3}; step h=2 should ≈ 2 * x_{n-1};
            // step h=3 should ≈ 2 * future_x[0].
            let preds = forecast.primary();
            assert_relative_eq!(preds[0], 2.0 * x_vals[n - 3], epsilon = 0.05);
            assert_relative_eq!(preds[2], 2.0 * x_vals[n - 1], epsilon = 0.05);
            assert_relative_eq!(preds[3], 2.0 * future_x[0], epsilon = 0.05);
        }

        #[test]
        fn ols_with_exog_rolling_emits_correct_columns() {
            // Sanity check that the column appears with the correct name and
            // that fit/predict produces finite output.
            let n = 60;
            let x_vals: Vec<f64> = (0..n).map(|i| i as f64).collect();
            let y: Vec<f64> = (0..n).map(|i| i as f64 * 0.5 + 1.0).collect();
            let cal = CalendarAnnotations::new().with_regressor("x".to_string(), x_vals.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_exog_rolling("x", 5, 0, RollingStatKind::Mean)
                    .unwrap(),
            );
            model.fit(&ts).unwrap();

            let names = model.features().feature_names(&["x".to_string()]);
            assert!(
                names.contains(&"x_roll_mean_w5_l0".to_string()),
                "expected x_roll_mean_w5_l0 in {:?}",
                names
            );

            let future_x: Vec<f64> = (n..n + 3).map(|i| i as f64).collect();
            let mut future_regs = HashMap::new();
            future_regs.insert("x".to_string(), future_x);
            let forecast = model.predict_with_exog(3, &future_regs).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn ols_with_exog_polynomial_recovers_quadratic() {
            // y = 0.5 * x^2 — needs a quadratic feature, not just x.
            let n = 80;
            let x_vals: Vec<f64> = (0..n).map(|i| (i as f64 * 0.05) - 2.0).collect();
            let y: Vec<f64> = x_vals.iter().map(|x| 0.5 * x * x).collect();
            let cal = CalendarAnnotations::new().with_regressor("x".to_string(), x_vals.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_exog_polynomial("x", 2)
                    .unwrap(),
            );
            model.fit(&ts).unwrap();

            let names = model.features().feature_names(&["x".to_string()]);
            assert!(
                names.contains(&"x_pow2".to_string()),
                "expected x_pow2 in {:?}",
                names
            );

            let future_x: Vec<f64> = (n..n + 5).map(|i| (i as f64 * 0.05) - 2.0).collect();
            let mut future_regs = HashMap::new();
            future_regs.insert("x".to_string(), future_x.clone());
            let forecast = model.predict_with_exog(5, &future_regs).unwrap();
            // Predictions should track 0.5 * x^2 closely
            for (h, &p) in forecast.primary().iter().enumerate() {
                let expected = 0.5 * future_x[h] * future_x[h];
                assert_relative_eq!(p, expected, epsilon = 0.05);
            }
        }

        #[test]
        fn ols_with_exog_interaction_recovers_product() {
            // y = 1.5 * a * b — pure interaction term.
            let n = 80;
            let a: Vec<f64> = (0..n).map(|i| (i as f64 * 0.1).sin()).collect();
            let b: Vec<f64> = (0..n).map(|i| (i as f64 * 0.13).cos()).collect();
            let y: Vec<f64> = (0..n).map(|i| 1.5 * a[i] * b[i]).collect();

            let cal = CalendarAnnotations::new()
                .with_regressor("a".to_string(), a.clone())
                .with_regressor("b".to_string(), b.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_exog_interaction("a", "b"),
            );
            model.fit(&ts).unwrap();

            let names = model.features().feature_names(&[]);
            assert!(
                names.contains(&"a_x_b".to_string()),
                "expected a_x_b in {:?}",
                names
            );

            let future_a: Vec<f64> = (n..n + 5).map(|i| (i as f64 * 0.1).sin()).collect();
            let future_b: Vec<f64> = (n..n + 5).map(|i| (i as f64 * 0.13).cos()).collect();
            let mut future_regs = HashMap::new();
            future_regs.insert("a".to_string(), future_a.clone());
            future_regs.insert("b".to_string(), future_b.clone());
            let forecast = model.predict_with_exog(5, &future_regs).unwrap();
            for (h, &p) in forecast.primary().iter().enumerate() {
                let expected = 1.5 * future_a[h] * future_b[h];
                assert_relative_eq!(p, expected, epsilon = 0.02);
            }
        }

        #[test]
        fn exog_polynomial_validates_degree_zero() {
            let f = RegressionFeatures::new().with_exog_polynomial("x", 0);
            assert!(f.is_err());
        }

        #[test]
        fn exog_polynomial_degree_one_emits_no_columns() {
            let feats = RegressionFeatures::new()
                .no_trend()
                .no_exog()
                .with_exog_polynomial("x", 1)
                .unwrap();
            // Only the polynomial spec is added, but degree=1 produces no columns.
            let names = feats.feature_names(&[]);
            assert!(
                names.iter().all(|n| !n.starts_with("x_pow")),
                "degree=1 should emit no columns; got {:?}",
                names
            );
        }

        #[test]
        fn categorical_one_hot_recovers_per_category_offsets() {
            // y = baseline + offset[category], with category in {0, 1, 2}
            let n = 60;
            let mut cat = Vec::with_capacity(n);
            let mut y = Vec::with_capacity(n);
            for i in 0..n {
                let c = (i % 3) as i64;
                cat.push(c as f64);
                y.push(match c {
                    0 => 1.0,
                    1 => 5.0,
                    2 => 9.0,
                    _ => unreachable!(),
                });
            }
            let cal = CalendarAnnotations::new().with_regressor("region".to_string(), cat.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_categorical(
                        "region",
                        vec![0, 1, 2],
                        CategoricalStrategy::OneHot { drop_first: false },
                    )
                    .unwrap(),
            );
            model.fit(&ts).unwrap();

            let names = model.features().feature_names(&[]);
            assert!(names.contains(&"region_eq_0".to_string()));
            assert!(names.contains(&"region_eq_1".to_string()));
            assert!(names.contains(&"region_eq_2".to_string()));

            // Forecast: cat values 0, 1, 2 → expect 1, 5, 9 (within OLS noise tolerance).
            let mut fr = HashMap::new();
            fr.insert("region".to_string(), vec![0.0, 1.0, 2.0]);
            let forecast = model.predict_with_exog(3, &fr).unwrap();
            let preds = forecast.primary();
            assert_relative_eq!(preds[0], 1.0, epsilon = 1e-6);
            assert_relative_eq!(preds[1], 5.0, epsilon = 1e-6);
            assert_relative_eq!(preds[2], 9.0, epsilon = 1e-6);
        }

        #[test]
        fn issue_113_categorical_excludes_raw_exog_for_same_column() {
            // Before #113 fix: declaring a column via with_categorical(...)
            // while keeping .exog() on caused the column to appear BOTH as
            // raw numeric AND as its one-hot dummies, producing perfect
            // collinearity. After the fix, the claimed column is excluded
            // from the raw exog block so the design stays well-conditioned
            // and a *different* raw regressor can still flow through.
            let n = 60;
            let mut cat = Vec::with_capacity(n);
            let mut promo = Vec::with_capacity(n);
            let mut y = Vec::with_capacity(n);
            for i in 0..n {
                let c = (i % 3) as i64;
                cat.push(c as f64);
                // independent binary regressor — bumps y by 2 when on
                let p = ((i / 4) % 2) as i64;
                promo.push(p as f64);
                y.push(
                    match c {
                        0 => 1.0,
                        1 => 5.0,
                        2 => 9.0,
                        _ => unreachable!(),
                    } + 2.0 * p as f64,
                );
            }
            let cal = CalendarAnnotations::new()
                .with_regressor("is_chp".to_string(), cat)
                .with_regressor("promo".to_string(), promo);
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            // .exog() ON, plus .with_categorical("is_chp", …)
            // — the pre-fix path would produce a singular X.
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_categorical(
                        "is_chp",
                        vec![0, 1, 2],
                        CategoricalStrategy::OneHot { drop_first: true },
                    )
                    .unwrap(),
            );
            model.fit(&ts).unwrap();

            // Feature-name sanity: raw `promo` present, raw `is_chp` absent,
            // is_chp dummies present.
            let names = model.features().feature_names(&["promo".to_string()]);
            assert!(
                names.contains(&"promo".to_string()),
                "expected raw 'promo' in feature names, got {:?}",
                names,
            );
            assert!(
                !names.contains(&"is_chp".to_string()),
                "raw 'is_chp' must be excluded (claimed by categorical), got {:?}",
                names,
            );
            assert!(
                names
                    .iter()
                    .any(|n| n == "is_chp_eq_1" || n == "is_chp_eq_2"),
                "expected at least one is_chp dummy in feature names, got {:?}",
                names,
            );

            // Forecast: should recover the four expected combos
            //   (cat, promo) → y: (0,0)→1, (1,1)→7, (2,0)→9, (0,1)→3
            let mut fr = HashMap::new();
            fr.insert("is_chp".to_string(), vec![0.0, 1.0, 2.0, 0.0]);
            fr.insert("promo".to_string(), vec![0.0, 1.0, 0.0, 1.0]);
            let forecast = model.predict_with_exog(4, &fr).unwrap();
            let preds = forecast.primary();
            assert_relative_eq!(preds[0], 1.0, epsilon = 1e-4);
            assert_relative_eq!(preds[1], 7.0, epsilon = 1e-4);
            assert_relative_eq!(preds[2], 9.0, epsilon = 1e-4);
            assert_relative_eq!(preds[3], 3.0, epsilon = 1e-4);

            // Validation: missing a claimed column at predict time must
            // surface an error (regression test for the extended
            // predict_with_exog validation in #113).
            let mut fr_missing = HashMap::new();
            fr_missing.insert("promo".to_string(), vec![0.0, 1.0, 0.0, 1.0]);
            // 'is_chp' deliberately absent
            let err = model.predict_with_exog(4, &fr_missing);
            assert!(
                err.is_err(),
                "predict_with_exog should error when a claimed column ('is_chp') is missing from future regressors",
            );
        }

        #[test]
        fn categorical_one_hot_drop_first() {
            let feats = RegressionFeatures::new()
                .no_trend()
                .no_exog()
                .with_categorical(
                    "c",
                    vec![10, 20, 30],
                    CategoricalStrategy::OneHot { drop_first: true },
                )
                .unwrap();
            let names = feats.feature_names(&[]);
            // First category dropped → only c_eq_20 and c_eq_30
            assert!(!names.contains(&"c_eq_10".to_string()));
            assert!(names.contains(&"c_eq_20".to_string()));
            assert!(names.contains(&"c_eq_30".to_string()));
        }

        #[test]
        fn categorical_ordinal_emits_index() {
            let n = 30;
            let cat: Vec<f64> = (0..n).map(|i| (i % 3) as f64).collect();
            let y: Vec<f64> = cat.clone();
            let cal = CalendarAnnotations::new().with_regressor("c".to_string(), cat.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_categorical("c", vec![0, 1, 2], CategoricalStrategy::Ordinal)
                    .unwrap(),
            );
            model.fit(&ts).unwrap();
            let names = model.features().feature_names(&[]);
            assert!(names.contains(&"c_ord".to_string()));
        }

        #[test]
        fn categorical_target_encoding_uses_per_category_mean() {
            // Means: cat0 → mean 0.0, cat1 → mean 10.0
            let n = 40;
            let mut cat = Vec::with_capacity(n);
            let mut y = Vec::with_capacity(n);
            for i in 0..n {
                let c = (i % 2) as i64;
                cat.push(c as f64);
                y.push(if c == 0 { 0.0 } else { 10.0 });
            }
            let cal = CalendarAnnotations::new().with_regressor("c".to_string(), cat.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_categorical(
                        "c",
                        vec![0, 1],
                        CategoricalStrategy::Target { smoothing: 0.0 },
                    )
                    .unwrap(),
            );
            model.fit(&ts).unwrap();

            // Target encoding maps cat 0 → 0.0 and cat 1 → 10.0; with no
            // smoothing the OLS regression on this single feature should
            // recover y exactly via slope=1, intercept=0.
            let mut fr = HashMap::new();
            fr.insert("c".to_string(), vec![0.0, 1.0]);
            let preds = model.predict_with_exog(2, &fr).unwrap();
            assert_relative_eq!(preds.primary()[0], 0.0, epsilon = 1e-6);
            assert_relative_eq!(preds.primary()[1], 10.0, epsilon = 1e-6);
        }

        #[test]
        fn categorical_count_encoding() {
            // 30 obs, cat=0 in 20, cat=1 in 10
            let n = 30;
            let mut cat = vec![0.0; 20];
            cat.extend(vec![1.0; 10]);
            let y: Vec<f64> = cat.iter().map(|&c| c + 1.0).collect();
            let cal = CalendarAnnotations::new().with_regressor("c".to_string(), cat.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_categorical("c", vec![0, 1], CategoricalStrategy::Count)
                    .unwrap(),
            );
            model.fit(&ts).unwrap();
            let names = model.features().feature_names(&[]);
            assert!(names.contains(&"c_count".to_string()));

            // Predict on cat values 0 and 1 — encoded as counts 20 and 10.
            // Should produce finite forecasts.
            let mut fr = HashMap::new();
            fr.insert("c".to_string(), vec![0.0, 1.0]);
            let preds = model.predict_with_exog(2, &fr).unwrap();
            for &p in preds.primary() {
                assert!(p.is_finite());
            }
        }

        #[test]
        fn categorical_unknown_category_handled() {
            // Train on cat in {0, 1}, predict on unknown cat 99.
            let n = 20;
            let cat: Vec<f64> = (0..n).map(|i| (i % 2) as f64).collect();
            let y: Vec<f64> = cat.clone();
            let cal = CalendarAnnotations::new().with_regressor("c".to_string(), cat.clone());
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(y)
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .no_exog()
                    .with_categorical(
                        "c",
                        vec![0, 1],
                        CategoricalStrategy::OneHot { drop_first: false },
                    )
                    .unwrap(),
            );
            model.fit(&ts).unwrap();
            // Future has cat=99 which is unknown — all dummy columns 0 → predict baseline.
            let mut fr = HashMap::new();
            fr.insert("c".to_string(), vec![99.0]);
            let preds = model.predict_with_exog(1, &fr).unwrap();
            assert!(preds.primary()[0].is_finite());
        }

        #[test]
        fn categorical_validates_inputs() {
            // empty categories
            assert!(RegressionFeatures::new()
                .with_categorical("c", vec![], CategoricalStrategy::Ordinal)
                .is_err());
            // duplicate categories
            assert!(RegressionFeatures::new()
                .with_categorical("c", vec![1, 1, 2], CategoricalStrategy::Ordinal)
                .is_err());
            // negative target smoothing
            assert!(RegressionFeatures::new()
                .with_categorical(
                    "c",
                    vec![0, 1],
                    CategoricalStrategy::Target { smoothing: -1.0 }
                )
                .is_err());
        }

        #[test]
        fn classify_features_marks_target_encoding_data_dependent() {
            let feats = RegressionFeatures::new()
                .no_trend()
                .no_exog()
                .with_categorical(
                    "c",
                    vec![0, 1],
                    CategoricalStrategy::Target { smoothing: 1.0 },
                )
                .unwrap();
            let classified = feats.classify_features(&[]);
            let target_safety = classified
                .iter()
                .find(|(name, _)| name == "c_target")
                .map(|(_, s)| *s);
            assert_eq!(target_safety, Some(FeatureSafety::DataDependent));
        }

        #[test]
        fn ols_exog_missing_regressor_errors() {
            let n = 30;
            let cal = CalendarAnnotations::new().with_regressor("x".to_string(), vec![1.0; n]);
            let ts = TimeSeriesBuilder::new()
                .timestamps(make_timestamps(n))
                .values(vec![1.0; n])
                .calendar(cal)
                .build()
                .unwrap();

            let mut model = RegressionForecaster::ols(RegressionFeatures::new().trend());
            model.fit(&ts).unwrap();

            let future_regs = HashMap::new(); // missing "x"
            assert!(model.predict_with_exog(5, &future_regs).is_err());
        }

        #[test]
        fn ols_name() {
            let model = RegressionForecaster::linear_trend();
            assert_eq!(model.name(), "OLS");
        }

        #[test]
        fn ols_residuals_sum_near_zero() {
            let ts = make_linear_ts(40);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let residuals = model.residuals().unwrap();
            let sum: f64 = residuals.iter().filter(|r| r.is_finite()).sum();
            assert!(sum.abs() < 1.0, "residuals sum = {}", sum);
        }

        #[test]
        fn ols_predict_with_intervals_linear_trend() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let forecast = model.predict_with_intervals(5, 0.95).unwrap();
            assert_eq!(forecast.primary().len(), 5);

            // Intervals should be present for a non-recursive model
            let lower = forecast.lower().expect("lower bounds should be present");
            let upper = forecast.upper().expect("upper bounds should be present");
            assert_eq!(lower[0].len(), 5);
            assert_eq!(upper[0].len(), 5);

            // Lower < point < upper
            for h in 0..5 {
                assert!(
                    lower[0][h] < forecast.primary()[h],
                    "lower[{}] = {} should be < point = {}",
                    h,
                    lower[0][h],
                    forecast.primary()[h],
                );
                assert!(
                    forecast.primary()[h] < upper[0][h],
                    "point = {} should be < upper[{}] = {}",
                    forecast.primary()[h],
                    h,
                    upper[0][h],
                );
            }
        }

        #[test]
        fn ols_predict_with_intervals_ar_falls_back() {
            // AR models use recursive prediction; intervals should be absent
            let n = 100;
            let mut values = vec![10.0];
            for i in 1..n {
                values.push(0.8 * values[i - 1] + 1.0 + 0.01 * (i as f64).sin());
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ar(1);
            model.fit(&ts).unwrap();

            let forecast = model.predict_with_intervals(5, 0.95).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            // Recursive model falls back to point-only
            assert!(forecast.lower().is_none());
            assert!(forecast.upper().is_none());
        }

        #[test]
        fn ols_r_squared() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();

            let r2 = model.r_squared().unwrap();
            assert!(
                r2 > 0.99,
                "R² should be near 1.0 for linear data, got {}",
                r2
            );
        }

        #[test]
        fn ols_insufficient_data() {
            let ts = TimeSeries::univariate(make_timestamps(2), vec![1.0, 2.0]).unwrap();
            let mut model = RegressionForecaster::ar(3);
            assert!(model.fit(&ts).is_err());
        }

        #[test]
        fn ols_no_features_errors() {
            let ts = make_linear_ts(30);
            let mut model =
                RegressionForecaster::ols(RegressionFeatures::new().no_trend().no_exog());
            assert!(model.fit(&ts).is_err());
        }

        #[test]
        fn ols_zero_horizon() {
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::linear_trend();
            model.fit(&ts).unwrap();
            let forecast = model.predict(0).unwrap();
            assert!(forecast.primary().is_empty());
        }

        #[test]
        fn ols_model_registry_integration() {
            use crate::models::{ModelRegistry, ModelSpec};

            let mut reg = ModelRegistry::new();
            reg.register(ModelSpec::new(
                "OLS(trend)",
                || Box::new(RegressionForecaster::linear_trend()),
                false,
            ));
            reg.register(ModelSpec::new(
                "OLS(AR3)",
                || Box::new(RegressionForecaster::ar(3)),
                false,
            ));

            assert_eq!(reg.len(), 2);

            let ts = make_linear_ts(50);
            for spec in reg.iter() {
                let mut model = spec.create();
                model.fit(&ts).unwrap();
                // AR model won't have exog, so predict should work
                if !model.has_exog() {
                    let fc = model.predict(5).unwrap();
                    assert_eq!(fc.primary().len(), 5);
                }
            }
        }

        // ── Component integration tests ─────────────────────────────

        fn make_seasonal_ts(n: usize, period: usize) -> TimeSeries {
            // y = 2*t + 10 * sin(2*pi*t/period) + 5
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let t = i as f64;
                    2.0 * t + 10.0 * (2.0 * std::f64::consts::PI * t / period as f64).sin() + 5.0
                })
                .collect();
            TimeSeries::univariate(make_timestamps(n), values).unwrap()
        }

        #[test]
        fn ols_fourier_seasonality() {
            let ts = make_seasonal_ts(100, 7);
            let mut model = RegressionForecaster::trend_fourier(7, 3);
            model.fit(&ts).unwrap();

            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Fourier prediction should be finite");
            }

            // R² should be high since y = trend + sin is well modeled by Fourier + trend
            let r2 = model.r_squared().unwrap();
            assert!(
                r2 > 0.95,
                "R² should be > 0.95 for sinusoidal data, got {}",
                r2
            );
        }

        #[test]
        fn ols_dummy_seasonal() {
            let ts = make_seasonal_ts(56, 7); // 8 full weeks
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .dummy_seasonal(7)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn ols_theilsen_trend_component() {
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_trend_component(TrendType::TheilSen)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            // Should produce increasing predictions for linear data
            assert!(forecast.primary()[4] > forecast.primary()[0]);
        }

        #[test]
        fn ols_exponential_trend_component() {
            // y = exp(0.05 * t)
            let n = 60;
            let values: Vec<f64> = (0..n).map(|i| (0.05 * i as f64).exp()).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_trend_component(TrendType::Exponential)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
                assert!(v > 0.0, "exponential trend should predict positive values");
            }
        }

        #[test]
        fn ols_trend_plus_fourier_plus_lags() {
            let ts = make_seasonal_ts(100, 7);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .lags(2)
                    .fourier(7, 3)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(10).unwrap();
            assert_eq!(forecast.primary().len(), 10);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn ols_quadratic_trend_component() {
            // y = t^2
            let n = 50;
            let values: Vec<f64> = (0..n).map(|i| (i as f64).powi(2)).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_trend_component(TrendType::Quadratic)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            for (h, &pred) in forecast.primary().iter().enumerate() {
                let expected = ((n + h) as f64).powi(2);
                assert_relative_eq!(pred, expected, epsilon = 5.0);
            }
        }

        #[test]
        fn ols_multiple_fourier_periods() {
            // Two seasonal components: period 7 + period 12
            let n = 120;
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let t = i as f64;
                    5.0 * (2.0 * std::f64::consts::PI * t / 7.0).sin()
                        + 3.0 * (2.0 * std::f64::consts::PI * t / 12.0).cos()
                        + 0.5 * t
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .fourier(7, 2)
                    .fourier(12, 2)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let r2 = model.r_squared().unwrap();
            assert!(
                r2 > 0.85,
                "R² should be high for dual-seasonal data, got {}",
                r2
            );

            let forecast = model.predict(12).unwrap();
            assert_eq!(forecast.primary().len(), 12);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn ols_component_only_no_trend_index() {
            // Use TheilSen component as the only feature (no raw trend index)
            let ts = make_linear_ts(40);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_trend_component(TrendType::Linear)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(3).unwrap();
            assert_eq!(forecast.primary().len(), 3);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        // ── Feature safety classification tests ─────────────────────

        #[test]
        fn classify_features_deterministic_only() {
            let features = RegressionFeatures::new().trend().fourier(7, 2).no_exog();
            let classified = features.classify_features(&[]);
            // __trend + 4 Fourier columns = 5 total
            assert_eq!(classified.len(), 5);
            assert!(classified
                .iter()
                .all(|(_, s)| *s == FeatureSafety::Deterministic));
        }

        #[test]
        fn classify_features_mixed() {
            let features = RegressionFeatures::new()
                .trend()
                .lags(2)
                .with_trend_component(TrendType::TheilSen)
                .fourier(7, 1)
                .dummy_seasonal(12)
                .with_changepoint_steps(vec![10, 50])
                .no_exog();
            let exog: Vec<String> = vec![];
            let classified = features.classify_features(&exog);
            // __trend (D) + __lag_1 (D) + __lag_2 (D) + __theilsen (DD) +
            // __fourier_p7_sin_1 (D) + __fourier_p7_cos_1 (D) + __seasonal_12 (DD) +
            // __cp_step_1 (S) + __cp_step_2 (S) = 9
            assert_eq!(classified.len(), 9);
            assert_eq!(
                classified[0],
                ("__trend".into(), FeatureSafety::Deterministic)
            );
            assert_eq!(
                classified[1],
                ("__lag_1".into(), FeatureSafety::Deterministic)
            );
            assert_eq!(
                classified[3],
                ("__theilsen_trend".into(), FeatureSafety::DataDependent)
            );
            assert_eq!(
                classified[4],
                ("__fourier_p7_sin_1".into(), FeatureSafety::Deterministic)
            );
            assert_eq!(
                classified[6],
                ("__seasonal_12".into(), FeatureSafety::DataDependent)
            );
            assert_eq!(
                classified[7],
                ("__cp_step_1".into(), FeatureSafety::Structural)
            );
            assert_eq!(
                classified[8],
                ("__cp_step_2".into(), FeatureSafety::Structural)
            );
        }

        #[test]
        fn classify_features_with_exog() {
            let features = RegressionFeatures::new().trend();
            let exog = vec!["temperature".to_string()];
            let classified = features.classify_features(&exog);
            assert_eq!(classified.last().unwrap().1, FeatureSafety::External);
        }

        // ── Changepoint feature tests ───────────────────────────────

        #[test]
        fn changepoint_step_functions_multiple_cps() {
            let cp = ChangepointFeature::step_functions(vec![10, 50]);
            assert_eq!(cp.column_names(), vec!["__cp_step_1", "__cp_step_2"]);
            let cols = cp.compute(100);
            assert_eq!(cols.len(), 2);
            // CP at 10: first 10 are 0, rest are 1
            assert_eq!(cols[0][9], 0.0);
            assert_eq!(cols[0][10], 1.0);
            assert_eq!(cols[0][99], 1.0);
            // CP at 50: first 50 are 0, rest are 1
            assert_eq!(cols[1][49], 0.0);
            assert_eq!(cols[1][50], 1.0);
        }

        #[test]
        fn changepoint_regime_index() {
            let cp = ChangepointFeature::new(vec![10, 50], ChangepointEncoding::RegimeIndex);
            assert_eq!(cp.column_names(), vec!["__cp_regime"]);
            let cols = cp.compute(100);
            assert_eq!(cols.len(), 1);
            assert_eq!(cols[0][0], 0.0); // before any CP
            assert_eq!(cols[0][10], 1.0); // after first CP
            assert_eq!(cols[0][49], 1.0); // still in regime 1
            assert_eq!(cols[0][50], 2.0); // after second CP
            assert_eq!(cols[0][99], 2.0);
        }

        #[test]
        fn changepoint_cumulative_count() {
            let cp = ChangepointFeature::new(vec![10, 50], ChangepointEncoding::CumulativeCount);
            assert_eq!(cp.column_names(), vec!["__cp_count"]);
            let cols = cp.compute(80);
            // Same numeric values as RegimeIndex
            assert_eq!(cols[0][5], 0.0);
            assert_eq!(cols[0][30], 1.0);
            assert_eq!(cols[0][60], 2.0);
        }

        #[test]
        fn changepoint_forward_fill_step_functions() {
            let cp = ChangepointFeature::step_functions(vec![10, 50]);
            let fill = cp.fill_values(80);
            assert_eq!(fill.len(), 2);
            // Both CPs are before n_train=80
            assert_eq!(fill[0], 1.0);
            assert_eq!(fill[1], 1.0);
        }

        #[test]
        fn changepoint_forward_fill_regime_index() {
            let cp = ChangepointFeature::new(vec![10, 50], ChangepointEncoding::RegimeIndex);
            let fill = cp.fill_values(80);
            assert_eq!(fill, vec![2.0]); // 2 CPs before index 80
        }

        #[test]
        fn changepoint_cp_after_training_end() {
            let cp = ChangepointFeature::step_functions(vec![100]);
            let fill = cp.fill_values(80);
            // CP at 100 is beyond n_train=80 → forward-fill is 0
            assert_eq!(fill, vec![0.0]);
        }

        #[test]
        fn changepoint_feature_in_regression() {
            // y = 2*t + 10 with a level shift of +20 at t=30
            let n = 60;
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let base = 2.0 * i as f64 + 10.0;
                    if i >= 30 {
                        base + 20.0
                    } else {
                        base
                    }
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .with_changepoint_steps(vec![30])
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let r2 = model.r_squared().unwrap();
            assert!(
                r2 > 0.99,
                "R² should be high with changepoint step, got {}",
                r2
            );

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
                // Predictions should be on the shifted level (> 130 for t=60+)
                assert!(v > 120.0, "predicted {} should be on shifted level", v);
            }
        }

        #[test]
        fn changepoint_regime_index_in_regression() {
            // y = 5*regime + noise
            let n = 60;
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let regime = if i < 20 {
                        0.0
                    } else if i < 40 {
                        1.0
                    } else {
                        2.0
                    };
                    5.0 * regime + 0.01 * (i as f64 * 0.5).sin()
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_changepoints(vec![20, 40], ChangepointEncoding::RegimeIndex)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
                // Should predict near 10.0 (regime 2 * 5.0)
                assert!(
                    (v - 10.0).abs() < 1.0,
                    "predicted {} should be near 10.0",
                    v
                );
            }
        }

        #[test]
        fn structural_feature_trait_custom() {
            // Custom outlier indicator: fills with 0.0 during prediction
            #[derive(Debug)]
            struct OutlierIndicator {
                outlier_indices: Vec<usize>,
            }
            impl StructuralFeature for OutlierIndicator {
                fn column_names(&self) -> Vec<String> {
                    vec!["__outlier".into()]
                }
                fn compute(&self, n: usize) -> Vec<Vec<f64>> {
                    let mut col = vec![0.0; n];
                    for &idx in &self.outlier_indices {
                        if idx < n {
                            col[idx] = 1.0;
                        }
                    }
                    vec![col]
                }
                fn fill_values(&self, _n_train: usize) -> Vec<f64> {
                    vec![0.0] // constant fill, not forward-fill
                }
                fn name(&self) -> &str {
                    "OutlierIndicator"
                }
            }

            let n = 50;
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let base = 0.5 * i as f64;
                    if i == 25 {
                        base + 100.0
                    } else {
                        base
                    }
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let indicator = Arc::new(OutlierIndicator {
                outlier_indices: vec![25],
            });

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .with_structural(indicator)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
                // Should not predict an outlier spike
                assert!(v < 50.0, "predicted {} should not include outlier", v);
            }
        }

        // ── Backend tests ───────────────────────────────────────────

        #[test]
        fn backend_ridge_fit_predict() {
            let ts = make_linear_ts(50);
            let mut model =
                RegressionForecaster::ridge(0.1, RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Ridge");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            // Ridge on clean linear data should track closely
            for (h, &pred) in forecast.primary().iter().enumerate() {
                let expected = 2.0 * (50 + h) as f64 + 10.0;
                assert_relative_eq!(pred, expected, epsilon = 2.0);
            }
        }

        #[test]
        fn backend_ridge_fourier() {
            let ts = make_seasonal_ts(100, 7);
            let mut model = RegressionForecaster::ridge(
                0.01,
                RegressionFeatures::new().trend().fourier(7, 3).no_exog(),
            );
            model.fit(&ts).unwrap();
            let r2 = model.r_squared().unwrap();
            assert!(r2 > 0.90, "Ridge R² = {} should be high", r2);
        }

        #[test]
        fn backend_elastic_net_fit_predict() {
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::elastic_net(
                0.01,
                0.5,
                RegressionFeatures::new().trend().fourier(7, 2).no_exog(),
            );
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "ElasticNet");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn backend_quantile_median() {
            let ts = make_linear_ts(60);
            let mut model =
                RegressionForecaster::quantile(0.5, RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Quantile");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            // Median of near-linear data should track the trend
            for (h, &pred) in forecast.primary().iter().enumerate() {
                let expected = 2.0 * (60 + h) as f64 + 10.0;
                assert_relative_eq!(pred, expected, epsilon = 3.0);
            }
        }

        #[test]
        fn backend_quantile_upper() {
            let ts = make_linear_ts(60);
            let mut m50 =
                RegressionForecaster::quantile(0.5, RegressionFeatures::new().trend().no_exog());
            let mut m90 =
                RegressionForecaster::quantile(0.9, RegressionFeatures::new().trend().no_exog());
            m50.fit(&ts).unwrap();
            m90.fit(&ts).unwrap();

            let f50 = m50.predict(1).unwrap();
            let f90 = m90.predict(1).unwrap();
            // For near-linear data, q90 ≥ q50 (approximately)
            assert!(
                f90.primary()[0] >= f50.primary()[0] - 1.0,
                "q90={} should be ≥ q50={}",
                f90.primary()[0],
                f50.primary()[0],
            );
        }

        #[test]
        fn backend_wls_decay() {
            let ts = make_linear_ts(60);
            let mut model =
                RegressionForecaster::wls_decay(0.95, RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "WLS");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn backend_wls_logistic_fit_predict() {
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::wls_logistic(
                12.0,
                RegressionFeatures::new().trend().no_exog(),
            );
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "WLS");
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn wls_logistic_weights_sigmoid_centred_at_offset_from_end() {
            // For n = 36, offset = 12, the sigmoid is centred at i = n - offset = 24.
            // Weight at i = 24 should be ~0.5; at i = n - 1 = 35 it should be near 1.0;
            // at i = 0 it should be near 0.0.
            //
            // We can only observe weights indirectly via the fit, so use a series whose
            // tail dominates: y[t] = 0 for t < n/2, y[t] = 100 for t >= n/2.
            // ExponentialDecay with low decay would also fit; here we just check the
            // logistic-fit produces a finite, sensible intercept dominated by the tail.
            let n = 36;
            let values: Vec<f64> = (0..n)
                .map(|i| if i < n / 2 { 0.0 } else { 100.0 })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();
            let mut model = RegressionForecaster::wls_logistic(
                6.0,
                RegressionFeatures::new().trend().no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(3).unwrap();
            // Tail-weighted fit on a step function ending at 100 should predict near 100.
            for &v in forecast.primary() {
                assert!(
                    (v - 100.0).abs() < 30.0,
                    "logistic tail-weighted forecast should sit near 100, got {}",
                    v
                );
            }
        }

        #[test]
        fn wls_logistic_ridge_recovers_linear_trend_with_recency() {
            // Sanity: combined backend produces finite forecasts and
            // recovers the linear trend slope on a clean series.
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::wls_logistic_ridge(
                12.0,
                0.1,
                RegressionFeatures::new().trend().no_exog(),
            );
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "WlsLogisticRidge");
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            // Should continue the trend y = 2t + 10 within shrinkage tolerance.
            for (h, &pred) in forecast.primary().iter().enumerate() {
                assert!(pred.is_finite());
                let expected = 2.0 * (60 + h) as f64 + 10.0;
                assert!(
                    (pred - expected).abs() < 3.0,
                    "forecast {} drifted from expected {} at h={}",
                    pred,
                    expected,
                    h
                );
            }
        }

        #[test]
        fn wls_logistic_ridge_stays_stable_on_heavy_features_low_effective_n() {
            // The motivating case from issue #103: small effective N
            // (logistic offset = 6 on a 30-row series — only the last ~6
            // rows have appreciable weight) with a heavy feature design
            // would blow up under pure WLS. Ridge regularisation must
            // keep predictions bounded.
            let n = 30;
            let timestamps = make_timestamps(n);
            // Synthetic series with a small linear trend + noise.
            let values: Vec<f64> = (0..n)
                .map(|i| 5.0 + 0.3 * i as f64 + (i as f64 * 0.7).sin())
                .collect();
            let ts = TimeSeries::univariate(timestamps, values).unwrap();
            // Heavy features: trend + 5 AR lags (≥ 6 columns, n_effective ≈ 6).
            let features = RegressionFeatures::new().trend().lags(5).no_exog();
            let mut model = RegressionForecaster::wls_logistic_ridge(6.0, 1.0, features);
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            // Reasonable y values are in [5, 14] — predictions should
            // stay in a sane neighbourhood, not blow up by orders of
            // magnitude as pure WLS would.
            for &v in forecast.primary() {
                assert!(v.is_finite(), "forecast non-finite: {}", v);
                assert!(
                    v.abs() < 100.0,
                    "WlsLogisticRidge forecast blew up: {} — regularisation didn't kick in",
                    v
                );
            }
        }

        #[test]
        fn wls_logistic_ridge_populates_finite_standard_errors() {
            // Issue #115: SE on WlsLogisticRidge must come from the
            // original-scale weighted covariance, not from the cloned
            // augmented-fit residual. After this fix they must be
            // finite, positive, and length-aligned with coefficients.
            use crate::models::inspect::{Explanation, Inspectable};

            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::wls_logistic_ridge(
                12.0,
                1.0,
                RegressionFeatures::new().trend().no_exog(),
            );
            model.fit(&ts).unwrap();

            let Explanation::Regression(e) = model.explanation().unwrap() else {
                panic!("expected Explanation::Regression");
            };

            assert_eq!(
                e.coef_std_errors.len(),
                e.coefficients.len(),
                "coef_std_errors must align with coefficients (got {}, {})",
                e.coef_std_errors.len(),
                e.coefficients.len(),
            );
            assert!(
                !e.coef_std_errors.is_empty(),
                "coef_std_errors must be populated for WlsLogisticRidge",
            );
            for (i, &se) in e.coef_std_errors.iter().enumerate() {
                assert!(se.is_finite(), "SE[{}] not finite: {}", i, se);
                assert!(se > 0.0, "SE[{}] not positive: {}", i, se);
            }
            assert!(
                e.intercept_std_error.is_finite(),
                "intercept SE not finite: {}",
                e.intercept_std_error,
            );
            assert!(
                e.intercept_std_error > 0.0,
                "intercept SE not positive: {}",
                e.intercept_std_error,
            );
        }

        #[test]
        fn wls_logistic_ridge_rejects_negative_lambda() {
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::wls_logistic_ridge(
                12.0,
                -0.1,
                RegressionFeatures::new().trend().no_exog(),
            );
            assert!(model.fit(&ts).is_err());
        }

        #[test]
        fn wls_logistic_rejects_non_finite_offset() {
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::wls_logistic(
                f64::NAN,
                RegressionFeatures::new().trend().no_exog(),
            );
            assert!(model.fit(&ts).is_err());
        }

        #[test]
        fn backend_rls_fit_predict() {
            let ts = make_linear_ts(60);
            let mut model =
                RegressionForecaster::rls(0.99, RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "RLS");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn backend_tweedie_gaussian() {
            // var_power=0 is Gaussian — should behave like OLS
            let ts = make_linear_ts(50);
            let mut model =
                RegressionForecaster::tweedie(0.0, RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Tweedie");

            let forecast = model.predict(5).unwrap();
            for (h, &pred) in forecast.primary().iter().enumerate() {
                let expected = 2.0 * (50 + h) as f64 + 10.0;
                assert_relative_eq!(pred, expected, epsilon = 2.0);
            }
        }

        #[test]
        fn backend_poisson_count_data() {
            // y = exp(0.02*t + 1) — count-like data
            let n = 60;
            let values: Vec<f64> = (0..n)
                .map(|i| (0.02 * i as f64 + 1.0).exp().round().max(1.0))
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model =
                RegressionForecaster::poisson(RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Poisson");

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
                assert!(v > 0.0, "Poisson should predict positive, got {}", v);
            }
        }

        #[test]
        fn backend_bls_nonnegative() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::nnls(RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "BLS");

            let forecast = model.predict(5).unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn backend_new_generic_constructor() {
            let ts = make_linear_ts(40);
            let mut model = RegressionForecaster::new(
                RegressionBackend::Ridge { lambda: 0.5 },
                RegressionFeatures::new().trend().no_exog(),
            );
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Ridge");
            assert!(model.r_squared().unwrap() > 0.9);
        }

        #[test]
        fn backend_model_registry_integration() {
            use crate::models::{ModelRegistry, ModelSpec};

            let mut reg = ModelRegistry::new();
            reg.register(ModelSpec::new(
                "OLS",
                || Box::new(RegressionForecaster::linear_trend()),
                false,
            ));
            reg.register(ModelSpec::new(
                "Ridge(0.1)",
                || {
                    Box::new(RegressionForecaster::ridge(
                        0.1,
                        RegressionFeatures::new().trend().no_exog(),
                    ))
                },
                false,
            ));
            reg.register(ModelSpec::new(
                "Quantile(0.5)",
                || {
                    Box::new(RegressionForecaster::quantile(
                        0.5,
                        RegressionFeatures::new().trend().no_exog(),
                    ))
                },
                false,
            ));

            let ts = make_linear_ts(50);
            for spec in reg.iter() {
                let mut model = spec.create();
                model.fit(&ts).unwrap();
                if !model.has_exog() {
                    let fc = model.predict(5).unwrap();
                    assert_eq!(fc.primary().len(), 5);
                }
            }
        }

        #[test]
        fn backend_dynamic_fit_predict() {
            let ts = make_linear_ts(60);
            let mut model =
                RegressionForecaster::dynamic(RegressionFeatures::new().trend().no_exog());
            model.fit(&ts).unwrap();
            assert_eq!(model.name(), "Dynamic");

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn backend_dynamic_smoothed() {
            let ts = make_seasonal_ts(100, 7);
            let mut model = RegressionForecaster::dynamic_smoothed(
                0.3,
                RegressionFeatures::new().trend().fourier(7, 2).no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        // ── Differencing tests ───────────────────────────────────────

        #[test]
        fn differencing_d1_produces_finite_forecast() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().trend().differencing(1).no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Forecast should be finite, got {}", v);
            }
        }

        #[test]
        fn differencing_d1_continues_trend() {
            // Linear trend: values 1..=30
            let ts = make_linear_ts(30);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().trend().differencing(1).no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();

            // After d=1 differencing of a linear trend, the differenced series
            // is constant. Integration should continue the trend roughly.
            let last = ts.primary_values().last().copied().unwrap();
            for &v in forecast.primary() {
                assert!(
                    v > last * 0.5,
                    "Forecast {} should continue trend from {}",
                    v,
                    last
                );
            }
        }

        #[test]
        fn differencing_d1_with_lags() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .differencing(1)
                    .lags(2)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn seasonal_differencing_produces_finite_forecast() {
            let ts = make_seasonal_ts(60, 7);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .seasonal_differencing(1, 7)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Forecast should be finite, got {}", v);
            }
        }

        #[test]
        fn both_differencing_and_seasonal() {
            let ts = make_seasonal_ts(80, 7);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .differencing(1)
                    .seasonal_differencing(1, 7)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn multiple_seasonal_differencing_periods() {
            // Daily data with weekly seasonality (simplified — use enough observations)
            let n = 120;
            let ts = make_seasonal_ts(n, 7);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .seasonal_differencing(1, 7)
                    .seasonal_differencing(1, 14)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Forecast should be finite, got {}", v);
            }
        }

        #[test]
        fn seasonal_differencing_chains_not_overwrites() {
            let features = RegressionFeatures::new()
                .seasonal_differencing(1, 7)
                .seasonal_differencing(1, 365);
            assert_eq!(features.seasonal_diffs.len(), 2);
            assert_eq!(features.seasonal_diffs[0], (1, 7));
            assert_eq!(features.seasonal_diffs[1], (1, 365));
        }

        #[test]
        fn no_differencing_is_default() {
            let features = RegressionFeatures::new();
            assert_eq!(features.diff_order, 0);
            assert!(features.seasonal_diffs.is_empty());
        }

        #[test]
        fn differencing_preserves_exog_regressors() {
            // Regression with trend + exog + differencing should use exog
            let n = 60;
            let timestamps = make_timestamps(n);
            let values: Vec<f64> = (0..n)
                .map(|i| 10.0 + 2.0 * i as f64 + 3.0 * (i % 2) as f64)
                .collect();
            let exog: Vec<f64> = (0..n).map(|i| (i % 2) as f64).collect();

            let mut ts = TimeSeries::univariate(timestamps, values).unwrap();
            let cal =
                crate::core::CalendarAnnotations::new().with_regressor("pulse".to_string(), exog);
            ts.set_calendar(cal);

            // Without differencing: should use exog
            let mut model_no_diff =
                RegressionForecaster::ols(RegressionFeatures::new().trend().exog());
            model_no_diff.fit(&ts).unwrap();
            assert!(
                model_no_diff.has_exog(),
                "model without diff should have exog"
            );

            // With differencing: should ALSO use exog (this was the bug)
            let mut model_diff =
                RegressionForecaster::ols(RegressionFeatures::new().trend().differencing(1).exog());
            model_diff.fit(&ts).unwrap();
            assert!(
                model_diff.has_exog(),
                "model with differencing should preserve exog regressors"
            );
        }

        // ── Auto-lag selection tests ─────────────────────────────────

        #[test]
        fn auto_lags_selects_reasonable_order() {
            // AR(2) process: y[t] = 0.5*y[t-1] + 0.3*y[t-2] + noise
            let mut values = vec![0.0; 100];
            values[0] = 1.0;
            values[1] = 0.5;
            for i in 2..100 {
                let noise = ((i * 7 + 3) % 11) as f64 * 0.02 - 0.11;
                values[i] = 0.5 * values[i - 1] + 0.3 * values[i - 2] + noise;
            }
            let ts = TimeSeries::univariate(make_timestamps(100), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().no_trend().auto_lags(10).no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn auto_lags_with_aic() {
            let ts = make_linear_ts(50);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .auto_lags_with(5, LagSelectionCriterion::Aic)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn auto_lags_max_zero_selects_no_lags() {
            let ts = make_linear_ts(50);
            let mut model =
                RegressionForecaster::ols(RegressionFeatures::new().trend().auto_lags(0).no_exog());
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
        }

        #[test]
        fn auto_lags_default_is_none() {
            let features = RegressionFeatures::new();
            assert!(features.auto_lag_config.is_none());
        }

        #[test]
        fn auto_lags_with_differencing() {
            let ts = make_linear_ts(60);
            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .auto_lags(5)
                    .differencing(1)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        // ── Extended auto-lag selection tests ────────────────────────

        #[test]
        fn auto_lags_ar1_selects_order_ge_1() {
            // Pure AR(1): y[t] = 0.8 * y[t-1] + noise
            let n = 200;
            let mut values = vec![0.0_f64; n];
            values[0] = 1.0;
            for i in 1..n {
                let noise = ((i * 13 + 7) % 17) as f64 * 0.01 - 0.085;
                values[i] = 0.8 * values[i - 1] + noise;
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().no_trend().auto_lags(10).no_exog(),
            );
            model.fit(&ts).unwrap();

            assert!(
                model.features().max_lag >= 1,
                "AR(1) process: auto_lags should select order >= 1, got {}",
                model.features().max_lag
            );
        }

        #[test]
        fn auto_lags_ar3_selects_order_ge_3() {
            // Pure AR(3): y[t] = 0.5*y[t-1] + 0.2*y[t-2] + 0.15*y[t-3] + noise
            let n = 300;
            let mut values = vec![0.0_f64; n];
            values[0] = 1.0;
            values[1] = 0.5;
            values[2] = 0.7;
            for i in 3..n {
                let noise = ((i * 11 + 5) % 19) as f64 * 0.005 - 0.0475;
                values[i] =
                    0.5 * values[i - 1] + 0.2 * values[i - 2] + 0.15 * values[i - 3] + noise;
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().no_trend().auto_lags(10).no_exog(),
            );
            model.fit(&ts).unwrap();

            assert!(
                model.features().max_lag >= 3,
                "AR(3) process: auto_lags should select order >= 3, got {}",
                model.features().max_lag
            );
        }

        #[test]
        fn auto_lags_white_noise_selects_low_order() {
            // White noise (no autocorrelation): BIC should select a low lag order.
            // Use a deterministic PRNG (xorshift-style) to avoid sequential patterns.
            let n = 200;
            let mut state: u64 = 123456789;
            let values: Vec<f64> = (0..n)
                .map(|_| {
                    // xorshift64
                    state ^= state << 13;
                    state ^= state >> 7;
                    state ^= state << 17;
                    // Map to [-1, 1]
                    (state as f64 / u64::MAX as f64) * 2.0 - 1.0
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().no_trend().auto_lags(10).no_exog(),
            );
            model.fit(&ts).unwrap();

            // For white noise, BIC should select 0 or a very low order (<=2).
            // The exact result depends on the PRNG realization, but it should
            // not select a high order like 8-10 on uncorrelated data.
            assert!(
                model.features().max_lag <= 2,
                "White noise: auto_lags should select a low order (<= 2), got {}",
                model.features().max_lag
            );
        }

        #[test]
        fn auto_lags_aic_selects_ge_bic() {
            // AIC penalizes less than BIC, so AIC should select >= BIC order.
            // Use an AR(2) process with moderate signal.
            let n = 200;
            let mut values = vec![0.0_f64; n];
            values[0] = 1.0;
            values[1] = 0.5;
            for i in 2..n {
                let noise = ((i * 7 + 3) % 11) as f64 * 0.02 - 0.11;
                values[i] = 0.5 * values[i - 1] + 0.3 * values[i - 2] + noise;
            }

            let ts_bic = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();
            let ts_aic = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model_bic = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .auto_lags_with(10, LagSelectionCriterion::Bic)
                    .no_exog(),
            );
            model_bic.fit(&ts_bic).unwrap();

            let mut model_aic = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .auto_lags_with(10, LagSelectionCriterion::Aic)
                    .no_exog(),
            );
            model_aic.fit(&ts_aic).unwrap();

            let bic_order = model_bic.features().max_lag;
            let aic_order = model_aic.features().max_lag;
            assert!(
                aic_order >= bic_order,
                "AIC order ({}) should be >= BIC order ({})",
                aic_order,
                bic_order
            );
        }

        #[test]
        fn auto_lags_short_series_does_not_panic() {
            // Barely enough data: 6 observations, max_lag = 4
            let n = 6;
            let values: Vec<f64> = (0..n).map(|i| i as f64 * 2.0 + 1.0).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().no_trend().auto_lags(4).no_exog(),
            );
            // Should not panic, even if it selects order 0
            model.fit(&ts).unwrap();
            let forecast = model.predict(3).unwrap();
            assert_eq!(forecast.primary().len(), 3);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        // ── Extended differencing round-trip tests ───────────────────

        #[test]
        fn differencing_d1_then_integrate_continues_original_scale() {
            // Linear trend: 1, 2, 3, ..., 50
            let n = 50;
            let values: Vec<f64> = (1..=n).map(|i| i as f64).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().trend().differencing(1).no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(5).unwrap();

            // After integration, forecast should be on the original scale (around 51-55)
            let last = *values.last().unwrap(); // 50.0
            for (h, &v) in forecast.primary().iter().enumerate() {
                assert!(
                    v.is_finite(),
                    "Forecast step {} should be finite, got {}",
                    h,
                    v
                );
                // Should be roughly continuing the trend (within 10% of expected)
                let expected = last + (h + 1) as f64;
                assert!(
                    (v - expected).abs() < expected * 0.2,
                    "Forecast step {}: expected ~{}, got {}",
                    h,
                    expected,
                    v
                );
            }
        }

        #[test]
        fn linear_trend_differencing_gives_constant_diffs() {
            // 1, 2, 3, ..., 50 — first difference should be all 1.0
            use crate::models::arima::difference;
            let values: Vec<f64> = (1..=50).map(|i| i as f64).collect();
            let diffed = difference(&values, 1);
            assert_eq!(diffed.len(), 49);
            for &d in &diffed {
                assert_relative_eq!(d, 1.0, epsilon = 1e-12);
            }
        }

        #[test]
        fn forecast_after_integration_is_on_original_scale() {
            // Exponential-ish growth: values in 100..~250 range
            let n = 60;
            let values: Vec<f64> = (0..n).map(|i| 100.0 + 2.5 * i as f64).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new().trend().differencing(1).no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(10).unwrap();

            let last = *values.last().unwrap(); // 100 + 2.5*59 = 247.5
            for &v in forecast.primary() {
                // Forecasts should be on the original scale (> 200), not on
                // the differenced scale (which would be around 2.5)
                assert!(
                    v > last * 0.5,
                    "Forecast {} should be on original scale (last train = {}), not differenced",
                    v,
                    last
                );
            }
        }

        // ── Extended seasonal differencing round-trip tests ──────────

        #[test]
        fn seasonal_diff_weekly_integration_recovers_scale() {
            // Weekly pattern repeating over 10 weeks
            let period = 7;
            let n_weeks = 10;
            let n = period * n_weeks;
            let weekly_pattern = [10.0, 20.0, 15.0, 25.0, 30.0, 12.0, 8.0];
            let values: Vec<f64> = (0..n)
                .map(|i| weekly_pattern[i % period] + 0.5 * i as f64)
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .seasonal_differencing(1, 7)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(7).unwrap();

            // Forecasts should be on the original scale, not the differenced scale
            let _last = *values.last().unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Seasonal forecast should be finite");
                // Original scale is roughly 10-80 range; differenced scale is ~3.5
                assert!(
                    v > 5.0,
                    "Forecast {} should be on original scale (not seasonal-differenced)",
                    v
                );
            }
        }

        #[test]
        fn seasonal_diff_monthly_integration_recovers_scale() {
            // Monthly pattern (period=12) with linear trend
            let period = 12;
            let n_years = 5;
            let n = period * n_years;
            let monthly_pattern = [
                5.0, 8.0, 12.0, 18.0, 22.0, 25.0, 24.0, 22.0, 18.0, 12.0, 8.0, 5.0,
            ];
            let values: Vec<f64> = (0..n)
                .map(|i| monthly_pattern[i % period] + 1.0 * i as f64)
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .seasonal_differencing(1, 12)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(12).unwrap();

            let last = *values.last().unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Monthly seasonal forecast should be finite");
                // Should be on the original scale (last train ~ 64)
                assert!(
                    v > last * 0.3,
                    "Forecast {} should be on original scale (last train = {})",
                    v,
                    last
                );
            }
        }

        #[test]
        fn multiple_seasonal_diffs_chain_correctly() {
            // Two seasonal diffs: weekly (7) and bi-weekly (14)
            let features = RegressionFeatures::new()
                .seasonal_differencing(1, 7)
                .seasonal_differencing(1, 365);
            assert_eq!(features.seasonal_diffs.len(), 2);
            assert_eq!(features.seasonal_diffs[0], (1, 7));
            assert_eq!(features.seasonal_diffs[1], (1, 365));

            // Verify that chaining actually works end-to-end with a model
            // Use enough data for both periods (need > 365 + 7 for both diffs)
            let n = 400;
            let values: Vec<f64> = (0..n)
                .map(|i| {
                    let t = i as f64;
                    // Weekly component + small trend
                    5.0 * (2.0 * std::f64::consts::PI * t / 7.0).sin() + 0.1 * t + 50.0
                })
                .collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .trend()
                    .seasonal_differencing(1, 7)
                    .seasonal_differencing(1, 14)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();
            let forecast = model.predict(7).unwrap();
            assert_eq!(forecast.primary().len(), 7);
            for &v in forecast.primary() {
                assert!(
                    v.is_finite(),
                    "Multi-seasonal forecast should be finite, got {}",
                    v
                );
            }
        }

        // ── Combined differencing + auto_lags tests ─────────────────

        #[test]
        fn differencing_d1_auto_lags_fits_on_differenced_data() {
            // Linear trend with AR structure: after differencing, should have lag structure
            let n = 100;
            let mut values = vec![0.0_f64; n];
            values[0] = 10.0;
            for i in 1..n {
                let noise = ((i * 13 + 7) % 17) as f64 * 0.01 - 0.085;
                values[i] = values[i - 1] + 2.0 + 0.3 * noise;
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .differencing(1)
                    .auto_lags(5)
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            // auto_lag_config is preserved for re-fit; max_lag is set to the selected order
            assert!(
                model.features().max_lag <= 5,
                "selected lag should be <= max_lag"
            );

            let forecast = model.predict(10).unwrap();
            assert_eq!(forecast.primary().len(), 10);

            let last = *values.last().unwrap();
            for &v in forecast.primary() {
                assert!(v.is_finite(), "Forecast should be finite");
                // Forecast should be on the original scale, not differenced
                assert!(
                    v > last * 0.5,
                    "Forecast {} should be on original scale (last train = {})",
                    v,
                    last
                );
            }
        }

        // ── Rolling / recursive feature tests ────────────────────────

        #[test]
        fn rolling_feature_rejects_zero_lag() {
            let err = RollingFeature::with_lag(5, 0, RollingStatKind::Mean);
            assert!(err.is_err());
            if let Err(ForecastError::InvalidParameter(msg)) = err {
                assert!(msg.contains("lag"));
            }
        }

        #[test]
        fn rolling_feature_rejects_zero_window() {
            assert!(RollingFeature::new(0, RollingStatKind::Mean).is_err());
        }

        #[test]
        fn rolling_feature_rejects_invalid_ewm_alpha() {
            assert!(RollingFeature::new(5, RollingStatKind::EwmMean { alpha: 0.0 }).is_err());
            assert!(RollingFeature::new(5, RollingStatKind::EwmMean { alpha: 1.5 }).is_err());
            assert!(RollingFeature::new(5, RollingStatKind::EwmStd { alpha: -0.1 }).is_err());
            assert!(RollingFeature::new(5, RollingStatKind::EwmMean { alpha: 1.0 }).is_ok());
        }

        #[test]
        fn rolling_feature_warmup_is_window_plus_lag_minus_one() {
            let f = RollingFeature::new(7, RollingStatKind::Mean).unwrap();
            assert_eq!(f.warmup(), 7);
            let f = RollingFeature::with_lag(5, 3, RollingStatKind::Mean).unwrap();
            assert_eq!(f.warmup(), 7); // 5 + 3 - 1
        }

        #[test]
        fn rolling_mean_compute_fit_matches_hand_calculation() {
            // series = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
            // window=3, lag=1 → at target_idx=3, block is values[0..3] = [1,2,3], mean = 2
            let values: Vec<f64> = (1..=10).map(|i| i as f64).collect();
            let feat = RollingFeature::new(3, RollingStatKind::Mean).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 3, &mut out);
            assert_relative_eq!(out[0], 2.0, epsilon = 1e-12);
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 4.0, epsilon = 1e-12); // mean(3,4,5)
            feat.compute_fit(&values, 9, &mut out);
            assert_relative_eq!(out[0], 8.0, epsilon = 1e-12); // mean(7,8,9)
        }

        #[test]
        fn rolling_mean_compute_predict_anchors_correctly() {
            // recent = [8, 9, 10] (last 3 training values), about to predict t=10
            // window=3, lag=1 → block is recent[0..3] = [8,9,10], mean = 9
            let recent = vec![8.0, 9.0, 10.0];
            let feat = RollingFeature::new(3, RollingStatKind::Mean).unwrap();
            let mut out = vec![0.0];
            feat.compute_predict(&recent, &mut out);
            assert_relative_eq!(out[0], 9.0, epsilon = 1e-12);

            // After predicting y_hat_0 = 11 and appending, next block is [9, 10, 11]
            let recent = vec![8.0, 9.0, 10.0, 11.0];
            feat.compute_predict(&recent, &mut out);
            assert_relative_eq!(out[0], 10.0, epsilon = 1e-12);
        }

        #[test]
        fn rolling_feature_with_lag_2_skips_a_step() {
            // window=3, lag=2 → at target=5, block is values[1..4] = [2,3,4], mean=3
            let values: Vec<f64> = (1..=10).map(|i| i as f64).collect();
            let feat = RollingFeature::with_lag(3, 2, RollingStatKind::Mean).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 3.0, epsilon = 1e-12);

            // recent = [..., 8, 9, 10] — for next prediction at t=10:
            // block is recent[len-2-3+1 .. len-2+1] = recent[len-4 .. len-1]
            let recent = vec![7.0, 8.0, 9.0, 10.0];
            feat.compute_predict(&recent, &mut out);
            // block is recent[0..3] = [7, 8, 9], mean = 8
            assert_relative_eq!(out[0], 8.0, epsilon = 1e-12);
        }

        #[test]
        fn rolling_std_and_var_match_definition() {
            // window = [1, 2, 3, 4, 5] — sample var = 2.5, std = sqrt(2.5)
            let values: Vec<f64> = (1..=6).map(|i| i as f64).collect();
            let feat = RollingFeature::new(5, RollingStatKind::Var).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 2.5, epsilon = 1e-10);

            let feat = RollingFeature::new(5, RollingStatKind::Std).unwrap();
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 2.5_f64.sqrt(), epsilon = 1e-10);
        }

        #[test]
        fn rolling_min_max_median_sum_basic() {
            let values: Vec<f64> = vec![3.0, 1.0, 4.0, 1.0, 5.0, 9.0];
            // At target=5, window=5, lag=1 → block is values[0..5] = [3,1,4,1,5]
            for (kind, expected) in [
                (RollingStatKind::Min, 1.0),
                (RollingStatKind::Max, 5.0),
                (RollingStatKind::Median, 3.0), // sorted: 1,1,3,4,5
                (RollingStatKind::Sum, 14.0),
            ] {
                let feat = RollingFeature::new(5, kind).unwrap();
                let mut out = vec![0.0];
                feat.compute_fit(&values, 5, &mut out);
                assert_relative_eq!(out[0], expected, epsilon = 1e-12);
            }
        }

        #[test]
        fn rolling_quantile_range_iqr() {
            let values = vec![3.0, 1.0, 4.0, 1.0, 5.0, 9.0];
            // window of values[0..5] = [3,1,4,1,5] → sorted [1,1,3,4,5]
            let cases = [
                (RollingStatKind::Quantile { tau: 0.0 }, 1.0),
                (RollingStatKind::Quantile { tau: 0.5 }, 3.0),
                (RollingStatKind::Quantile { tau: 1.0 }, 5.0),
                (RollingStatKind::Range, 4.0),
                // Q3 of [1,1,3,4,5]: pos = 0.75 * 4 = 3.0 → 4.0
                // Q1 of [1,1,3,4,5]: pos = 0.25 * 4 = 1.0 → 1.0
                (RollingStatKind::Iqr, 3.0),
            ];
            for (kind, expected) in cases {
                let feat = RollingFeature::new(5, kind).unwrap();
                let mut out = vec![0.0];
                feat.compute_fit(&values, 5, &mut out);
                assert_relative_eq!(out[0], expected, epsilon = 1e-12);
            }
        }

        #[test]
        fn rolling_skew_kurt_zero_on_symmetric_window() {
            let values = vec![-2.0, -1.0, 0.0, 1.0, 2.0, 99.0]; // last is irrelevant; window covers [0..5]
            let feat = RollingFeature::new(5, RollingStatKind::Skew).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 0.0, epsilon = 1e-12);

            let feat = RollingFeature::new(5, RollingStatKind::Kurt).unwrap();
            feat.compute_fit(&values, 5, &mut out);
            // Uniform discrete: excess kurtosis ≈ -1.3 for n=5.
            assert!(
                out[0] < 0.0,
                "uniform window should have negative excess kurt"
            );
        }

        #[test]
        fn rolling_slope_recovers_linear_trend() {
            // values[0..5] = [0, 2, 4, 6, 8] → slope = 2
            let values = vec![0.0, 2.0, 4.0, 6.0, 8.0, 999.0];
            let feat = RollingFeature::new(5, RollingStatKind::Slope).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 2.0, epsilon = 1e-12);
        }

        #[test]
        fn rolling_rank_and_zscore() {
            // window [1,2,3,4,5] — last value 5 is the max
            let values = vec![1.0, 2.0, 3.0, 4.0, 5.0, 999.0];
            let feat = RollingFeature::new(5, RollingStatKind::Rank).unwrap();
            let mut out = vec![0.0];
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 1.0, epsilon = 1e-12); // 4/4 strictly less

            let feat = RollingFeature::new(5, RollingStatKind::ZScore).unwrap();
            feat.compute_fit(&values, 5, &mut out);
            // mean = 3, std (ddof=1) = sqrt(2.5) ≈ 1.5811 → z = 2/1.5811 ≈ 1.265
            assert_relative_eq!(out[0], 2.0 / 2.5_f64.sqrt(), epsilon = 1e-12);
        }

        #[test]
        fn rolling_count_above_below() {
            let values = vec![1.0, 2.0, 3.0, 4.0, 5.0, 999.0];
            let mut out = vec![0.0];

            let feat =
                RollingFeature::new(5, RollingStatKind::CountAbove { threshold: 2.5 }).unwrap();
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 3.0, epsilon = 1e-12); // 3, 4, 5

            let feat =
                RollingFeature::new(5, RollingStatKind::CountBelow { threshold: 2.5 }).unwrap();
            feat.compute_fit(&values, 5, &mut out);
            assert_relative_eq!(out[0], 2.0, epsilon = 1e-12); // 1, 2
        }

        #[test]
        fn rolling_quantile_validates_tau() {
            assert!(RollingFeature::new(5, RollingStatKind::Quantile { tau: -0.1 }).is_err());
            assert!(RollingFeature::new(5, RollingStatKind::Quantile { tau: 1.5 }).is_err());
            assert!(RollingFeature::new(5, RollingStatKind::Quantile { tau: f64::NAN }).is_err());
        }

        #[test]
        fn event_distance_fit_since_until() {
            let feat = EventDistanceFeature::new(vec![5, 12, 30], EventDistanceMode::Both);
            let mut out = vec![0.0; 2];
            feat.compute_fit(&[], 7, &mut out);
            assert_relative_eq!(out[0], 2.0, epsilon = 1e-12); // since: 7 - 5
            assert_relative_eq!(out[1], 5.0, epsilon = 1e-12); // until: 12 - 7

            // exact event hit → since = 0, until = next event - t
            feat.compute_fit(&[], 12, &mut out);
            assert_relative_eq!(out[0], 0.0, epsilon = 1e-12);
            assert_relative_eq!(out[1], 18.0, epsilon = 1e-12);

            // before any event → since = cap (default = u32::MAX as f64)
            feat.compute_fit(&[], 0, &mut out);
            assert_eq!(out[0], f64::from(u32::MAX));
            assert_relative_eq!(out[1], 5.0, epsilon = 1e-12);
        }

        #[test]
        fn event_distance_predict_uses_absolute_index() {
            let feat =
                EventDistanceFeature::new(vec![100, 110], EventDistanceMode::Since).with_cap(999);
            let mut out = vec![0.0];
            // n_train = 95, h = 4 → t = 99 → since = cap (no prior event)
            feat.compute_predict_at(&[], 4, 95, &mut out);
            assert_relative_eq!(out[0], 999.0, epsilon = 1e-12);
            // n_train = 95, h = 5 → t = 100 → since = 0
            feat.compute_predict_at(&[], 5, 95, &mut out);
            assert_relative_eq!(out[0], 0.0, epsilon = 1e-12);
            // n_train = 95, h = 8 → t = 103 → since = 3
            feat.compute_predict_at(&[], 8, 95, &mut out);
            assert_relative_eq!(out[0], 3.0, epsilon = 1e-12);
        }

        #[test]
        fn event_distance_dedupes_and_sorts() {
            let feat = EventDistanceFeature::new(vec![30, 10, 30, 20], EventDistanceMode::Until);
            assert_eq!(feat.events(), &[10, 20, 30]);
        }

        #[test]
        fn ewm_mean_converges_on_constant_input() {
            let feat = RollingFeature::new(20, RollingStatKind::EwmMean { alpha: 0.3 }).unwrap();
            let values = vec![7.0; 25];
            let mut out = vec![0.0];
            feat.compute_fit(&values, 20, &mut out);
            // Constant input → EWM should equal the constant exactly.
            assert_relative_eq!(out[0], 7.0, epsilon = 1e-12);
        }

        #[test]
        fn rolling_feature_end_to_end_fit_and_predict() {
            // y[t] = rolling_mean(y[t-3..t]) + 2
            // so after warmup, the pattern is predictable.
            let n = 60;
            let mut values = vec![10.0_f64, 10.0, 10.0];
            for _ in 3..n {
                let k = values.len();
                let rm = (values[k - 3] + values[k - 2] + values[k - 1]) / 3.0;
                values.push(rm + 2.0);
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values.clone()).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .with_rolling_mean(3)
                    .unwrap()
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(5).unwrap();
            assert_eq!(forecast.primary().len(), 5);
            for &v in forecast.primary() {
                assert!(v.is_finite());
                assert!(
                    v > *values.last().unwrap() * 0.5,
                    "forecast should track the growing pattern"
                );
            }

            // Recursive semantics: h=1 rolling-mean feature should equal
            // mean of (tail[-2], tail[-1], y_hat_0).
            // We can't introspect feature columns from Forecast, but we can
            // check that forecasts are monotonically increasing under this
            // pattern (each step adds ~2).
            let preds = forecast.primary();
            for i in 1..preds.len() {
                assert!(
                    preds[i] > preds[i - 1],
                    "pattern forecasts should grow; got {:?}",
                    preds
                );
            }
        }

        #[test]
        fn rolling_feature_combined_with_lags() {
            // Mix lag-1 and rolling_mean(5) features on a noisy AR(1)
            let n = 120;
            let mut values = vec![0.0_f64; n];
            values[0] = 1.0;
            for i in 1..n {
                let noise = ((i * 13 + 7) % 17) as f64 * 0.01 - 0.085;
                values[i] = 0.8 * values[i - 1] + noise;
            }
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let mut model = RegressionForecaster::ols(
                RegressionFeatures::new()
                    .no_trend()
                    .lags(1)
                    .with_rolling_mean(5)
                    .unwrap()
                    .no_exog(),
            );
            model.fit(&ts).unwrap();

            let forecast = model.predict(10).unwrap();
            assert_eq!(forecast.primary().len(), 10);
            for &v in forecast.primary() {
                assert!(v.is_finite());
            }
        }

        #[test]
        fn rolling_feature_warmup_absorbs_lag_offset() {
            // window=10, lag=1 → warmup=10. max_lag=3. offset should be max(3,10)=10.
            let feats = RegressionFeatures::new()
                .no_trend()
                .lags(3)
                .with_rolling_mean(10)
                .unwrap()
                .no_exog();
            assert_eq!(feats.max_effective_lag(), 3);
            assert_eq!(feats.max_recursive_warmup(), 10);
            assert_eq!(feats.lag_offset(), 10);
        }

        #[test]
        fn rolling_feature_classify_as_deterministic() {
            let feats = RegressionFeatures::new()
                .no_trend()
                .with_rolling_mean(5)
                .unwrap()
                .no_exog();
            let classified = feats.classify_features(&[]);
            assert_eq!(classified.len(), 1);
            assert_eq!(classified[0].1, FeatureSafety::Deterministic);
            assert!(classified[0].0.starts_with("__rolling_mean"));
        }

        #[test]
        fn rolling_feature_cross_validation_round_trip() {
            // The feature must survive rolling-origin cross-validation — each
            // fold creates a fresh forecaster via the factory closure.
            use crate::utils::cross_validation::{cross_validate, CVConfig};

            let n = 80;
            let values: Vec<f64> = (0..n).map(|i| (i as f64 * 0.1).sin() + 5.0).collect();
            let ts = TimeSeries::univariate(make_timestamps(n), values).unwrap();

            let factory = || {
                RegressionForecaster::ols(
                    RegressionFeatures::new()
                        .no_trend()
                        .with_rolling_mean(5)
                        .unwrap()
                        .no_exog(),
                )
            };

            let cv_config = CVConfig::expanding(50, 3);
            let result = cross_validate(&cv_config, &ts, factory);
            assert!(result.is_ok(), "CV failed: {:?}", result.err());
            let results = result.unwrap();
            assert!(results.n_folds >= 1);
        }
    }
}

#[cfg(feature = "postprocess")]
pub use ols_impl::{
    CategoricalStrategy, ChangepointEncoding, ChangepointFeature, EventDistanceFeature,
    EventDistanceMode, ExogFeatureSpec, FeatureSafety, LagSelectionCriterion, RecursiveFeature,
    RegressionBackend, RegressionFeatures, RegressionForecaster, RollingFeature, RollingStatKind,
    SeasonalSpec, StructuralFeature, TrendType, WeightStrategy,
};
// Re-export InformationCriterion so users can configure Dynamic backend without
// depending on anofox-regression directly.
#[cfg(feature = "postprocess")]
pub use anofox_regression::solvers::InformationCriterion;
