<!-- AUTO-GENERATED FILE. DO NOT EDIT MANUALLY. -->
<!-- Source: tools/ai/llms-base.txt + skills/ -->
<!-- Regenerate with: python tools/ai/generate_ai_context_files.py -->

# Skforecast

> Python library for time series forecasting using scikit-learn compatible models, statistical methods, and foundation models

This document is for skforecast v0.23.0+. If you are using an older version, check the documentation at skforecast.org.

Skforecast is a Python library for time series forecasting using scikit-learn compatible models, statistical methods, and foundation models. It works with any estimator compatible with the scikit-learn API (LightGBM, XGBoost, CatBoost, Keras, etc.).

## Quick Info

- Version: 0.23.0
- License: BSD-3-Clause
- Python: 3.10, 3.11, 3.12, 3.13, 3.14
- Repository: https://github.com/skforecast/skforecast
- Documentation: https://skforecast.org
- PyPI: https://pypi.org/project/skforecast/

## Installation

```bash
pip install skforecast
```

Optional dependencies:
```bash
pip install skforecast[stats]        # For ARIMA, SARIMAX, ETS models
pip install skforecast[plotting]     # For visualization
pip install skforecast[deeplearning] # For RNN/LSTM models
```

## Project Structure

```
skforecast/
├── base/                    # ForecasterBase - abstract parent class for all forecasters
├── recursive/               # ForecasterRecursive, ForecasterRecursiveMultiSeries,
│                            # ForecasterRecursiveClassifier, ForecasterStats, ForecasterEquivalentDate
├── direct/                  # ForecasterDirect, ForecasterDirectMultiVariate
├── deep_learning/           # ForecasterRnn, create_and_compile_model
├── foundation/              # FoundationModel, ForecasterFoundation
│                            # (zero-shot: Chronos-2, TimesFM 2.5, Moirai-2, TabICL, TabPFN-TS, TFC-T0)
├── stats/                   # Arima, Sarimax, Ets, Arar, acf, pacf, calculate_lag_autocorrelation
├── preprocessing/           # TimeSeriesDifferentiator, RollingFeatures, CalendarFeatures,
│                            # QuantileBinner, ConformalIntervalCalibrator, reshape_* functions
├── model_selection/         # backtesting_forecaster, grid/random/bayesian search, TimeSeriesFold
├── feature_selection/       # select_features, select_features_multiseries
├── metrics/                 # MASE, RMSSE, sMAPE, CRPS, coverage, pinball loss
├── datasets/                # 30+ built-in datasets (fetch_dataset, load_demo_dataset)
├── drift_detection/         # RangeDriftDetector, PopulationDriftDetector
├── utils/                   # Shared validation and transformation functions
├── exceptions/              # Custom warnings and exceptions
├── plot/                    # plot_residuals, plot_prediction_intervals, plot_prediction_distribution,
│                            # plot_multivariate_time_series_corr, set_dark_theme, backtesting_gif_creator
└── experimental/            # Experimental features (API may change)
```

### Module Relationships

- **Forecasters inheriting from `ForecasterBase`**: ForecasterRecursive, ForecasterRecursiveMultiSeries, ForecasterRecursiveClassifier, ForecasterDirect, ForecasterDirectMultiVariate, ForecasterRnn
- **Standalone forecasters (no inheritance)**: ForecasterStats, ForecasterEquivalentDate, ForecasterFoundation
- Statistical models in `stats/` are wrapped by `ForecasterStats` (in `recursive/`)
- `ForecasterFoundation` (in `foundation/`) wraps a `FoundationModel`, which delegates to an adapter class (`ChronosAdapter`, `TimesFMAdapter`, `MoiraiAdapter`, `TabICLAdapter`, `TabPFNAdapter`, `T0Adapter`) resolved from the HuggingFace `model_id`
- `model_selection/` functions work with all forecaster types
- `preprocessing/` classes can be passed to forecasters via `transformer_y`, `transformer_exog`, `window_features`

## Core Forecasters

| Forecaster | Use Case |
|------------|----------|
| ForecasterRecursive | Single series, recursive multi-step forecasting |
| ForecasterDirect | Single series, direct multi-step forecasting |
| ForecasterRecursiveMultiSeries | Multiple series forecasting (global model) |
| ForecasterDirectMultiVariate | Multivariate forecasting (multiple series as features) |
| ForecasterRnn | Deep learning (RNN/LSTM) forecasting |
| ForecasterStats | Statistical models (ARIMA, SARIMAX, ETS, ARAR) |
| ForecasterFoundation | Zero-shot forecasting with pre-trained foundation models (Chronos-2, TimesFM 2.5, Moirai-2, TabICL, TabPFN-TS, TFC-T0) |
| ForecasterRecursiveClassifier | Classification-based forecasting |
| ForecasterEquivalentDate | Baseline forecaster using equivalent past dates |

## Basic Usage Example

```python
# Single series forecasting with ForecasterRecursive
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

# Load data: y must be a pandas Series with DatetimeIndex and frequency set
data = pd.read_csv('data.csv', index_col='date', parse_dates=True)
data = data.asfreq('h')  # IMPORTANT: always set frequency before using skforecast

# Create and train forecaster
forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=24  # Use last 24 observations as features
)
forecaster.fit(y=data['target'])

# Predict next 10 steps
predictions = forecaster.predict(steps=10)

# Define cross-validation strategy for backtesting
cv = TimeSeriesFold(
    steps=10,
    initial_train_size=len(data) - 100,
    refit=False,
    fixed_train_size=False
)

# Backtesting for model evaluation
metric, predictions_backtest = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error'
)
```

## Multi-Series Forecasting Example

```python
# Multiple series with global model
from skforecast.recursive import ForecasterRecursiveMultiSeries
from lightgbm import LGBMRegressor

# Data: DataFrame with multiple series as columns
# series = pd.DataFrame({'series_1': [...], 'series_2': [...], ...})

forecaster = ForecasterRecursiveMultiSeries(
    estimator=LGBMRegressor(n_estimators=100, random_state=123),
    lags=24,
    encoding='ordinal'  # 'ordinal', 'ordinal_category', 'onehot', or None
)
forecaster.fit(series=series)

# Predict all series
predictions = forecaster.predict(steps=10)

# Predict specific series
predictions = forecaster.predict(steps=10, levels=['series_1', 'series_2'])
```

## With Exogenous Variables

```python
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24
)

# Fit with exogenous variables
forecaster.fit(y=y_train, exog=exog_train)

# Predict - exog must cover the forecast horizon
predictions = forecaster.predict(steps=10, exog=exog_test)
```

## Categorical Exogenous Variables

All ML forecasters (ForecasterRecursive, ForecasterDirect, ForecasterRecursiveMultiSeries, ForecasterDirectMultiVariate, ForecasterRecursiveClassifier) include a `categorical_features` parameter to handle categorical exogenous variables automatically.

```python
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    categorical_features='auto',  # Default: auto-detect non-numeric columns after transformer_exog
)
forecaster.fit(y=y_train, exog=exog_train)
```

**`categorical_features` options:**
- `'auto'` (default): Non-numeric columns (after `transformer_exog`) are automatically detected and encoded using an internal `OrdinalEncoder`. Native categorical support is configured automatically for compatible estimators (LightGBM, CatBoost, XGBoost, HistGradientBoostingRegressor).
- `list`: Explicit list of column names to treat as categorical (including numeric columns).
- `None`: No categorical encoding is applied.

**Important:** When `categorical_features` is not `None` do not set categorical features directly on the estimator or via `fit_kwargs`. The forecaster manages categorical configuration internally and overwrites any estimator-level settings.

**Choosing an encoding strategy:**

| Method | API | Best for |
|--------|-----|----------|
| Built-in `categorical_features` | `categorical_features='auto'` or `list` | Gradient boosting (LightGBM, XGBoost, CatBoost, HistGBR), simplest workflow |
| One-hot / Ordinal encoding | `transformer_exog` | Linear models, SVMs, non-gradient-boosting trees |
| Target encoding | Outside forecaster | High-cardinality features (applied manually to avoid leakage) |

**`categorical_features` and `transformer_exog` interaction:**
`transformer_exog` is applied **before** `categorical_features` detection. They can be used together, e.g., scale numeric columns with `transformer_exog` while `categorical_features='auto'` handles the categorical ones. Avoid applying both mechanisms to the same columns.

```python
from sklearn.compose import make_column_transformer
from sklearn.preprocessing import StandardScaler

# Scale numeric columns, leave categorical columns untouched
transformer_exog = make_column_transformer(
    (StandardScaler(), ['temp', 'hum']),
    remainder='passthrough',
    verbose_feature_names_out=False
).set_output(transform='pandas')

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    transformer_exog=transformer_exog,
    categorical_features='auto',  # Detects remaining non-numeric columns
)
```

## Handling Missing Values

All ML forecasters (ForecasterRecursive, ForecasterDirect, ForecasterRecursiveMultiSeries, ForecasterDirectMultiVariate, ForecasterRecursiveClassifier) include a `dropna_from_series` parameter to control how NaN values in the training data are handled. Not applicable to ForecasterStats, ForecasterEquivalentDate, or ForecasterRnn.

When `y`, `series`, or `exog` contain interspersed NaN values, the training matrices (`X_train`, `y_train`) will also contain NaNs. The `dropna_from_series` parameter determines what happens next.

**`dropna_from_series` options:**
- `False` (default): NaN rows are kept in the training matrices. A `MissingValuesWarning` is issued. Only use with NaN-tolerant estimators.
- `True`: Rows containing NaN in `X_train` (and the corresponding rows in `y_train`) are dropped before fitting.

**NaN-tolerant estimators** (support `dropna_from_series=False`):
- LightGBM (`LGBMRegressor`, `LGBMClassifier`)
- CatBoost (`CatBoostRegressor`, `CatBoostClassifier`)
- HistGradientBoosting (`HistGradientBoostingRegressor`, `HistGradientBoostingClassifier`)
- XGBoost (`XGBRegressor`, `XGBClassifier`) with `tree_method='hist'`

**Choosing a strategy:**

| Scenario | Recommended strategy |
|:---------|:---------------------|
| Using LightGBM, CatBoost, XGBoost (hist), or HistGradientBoosting | `dropna_from_series=False`: let the estimator handle NaNs. No data loss. |
| Estimator does not support NaN (RandomForest, LinearRegression, SVR...) | `dropna_from_series=True`: drop incomplete rows before fitting. |
| Few or small gaps in the series | `dropna_from_series=True`: simple and data loss is negligible. |
| Large gaps, need to preserve all training data | Imputation + `weight_func`: fill NaNs and down-weight imputed observations. |
| Multi-series with different lengths | `ForecasterRecursiveMultiSeries` with `dropna_from_series=True`. |

```python
from lightgbm import LGBMRegressor

# NaN-tolerant estimator: keep all rows including NaN
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(random_state=123, verbose=-1),
    lags=14,
    dropna_from_series=False,  # Default: NaN rows kept for NaN-tolerant estimators
)
forecaster.fit(y=y_train, suppress_warnings=True)  # Suppress MissingValuesWarning

# Non-NaN-tolerant estimator: drop rows with NaN
forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(random_state=123),
    lags=14,
    dropna_from_series=True,  # Drop NaN rows before fitting
)
forecaster.fit(y=y_train)
```

## Window Features (Rolling Statistics)

```python
from skforecast.preprocessing import RollingFeatures

# Create rolling features
rolling_features = RollingFeatures(
    stats=['mean', 'std', 'min', 'max'],
    window_sizes=7  # int applies to all stats, or list with same length as stats
)

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    window_features=rolling_features
)
```

## Prediction Intervals

```python
# Predict with confidence intervals
# Since v0.23.0 `interval` is expressed as quantiles in (0, 1), not percentiles
predictions = forecaster.predict_interval(
    steps=10,
    interval=[0.1, 0.9],  # 80% prediction interval
    method='bootstrapping',  # or 'conformal'
    n_boot=500
)
# Returns DataFrame with columns: pred, lower_bound, upper_bound
```

## Backtesting

Backtesting evaluates forecaster performance using time series cross-validation.

```python
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

# Define cross-validation strategy
cv = TimeSeriesFold(
    steps=10,
    initial_train_size=len(data) - 100,
    refit=False,
    fixed_train_size=False
)

metric, predictions = backtesting_forecaster(
    forecaster=forecaster,                      # Forecaster object to evaluate
    y=data['target'],                           # Time series data (pandas Series with DatetimeIndex)
    cv=cv,                                      # TimeSeriesFold with CV configuration
    metric='mean_absolute_error',               # Metric(s): str, callable, or list
    exog=exog,                                  # Exogenous variables (optional)
    interval=[0.1, 0.9],                        # Prediction intervals as quantiles in (0, 1) (optional)
    interval_method='bootstrapping',            # 'bootstrapping' or 'conformal'
    n_boot=250,                                 # Bootstrap iterations (only if method='bootstrapping')
    use_in_sample_residuals=True,               # Use training residuals for intervals
    use_binned_residuals=True,                  # Select residuals based on predicted values
    random_state=123,                           # Seed for reproducibility
    return_predictors=False,                    # Return predictor values used in each fold
    n_jobs='auto',                              # Parallel jobs (-1 for all cores, 'auto' for automatic)
    verbose=False,                              # Print fold information
    show_progress=True,                         # Show progress bar
    suppress_warnings=False                     # Suppress skforecast warnings
)
```

## Cross-Validation Strategies

Skforecast provides two cross-validation strategies for time series: `TimeSeriesFold` for multi-step ahead validation and `OneStepAheadFold` for one-step ahead validation.

### TimeSeriesFold

Class to split time series data into train and test folds for backtesting and hyperparameter search.

```python
from skforecast.model_selection import TimeSeriesFold

cv = TimeSeriesFold(
    steps=12,                        # (required) Number of observations to predict in each fold (forecast horizon)
    initial_train_size=100,          # Number of observations for initial training. Can be int, str (date), or pd.Timestamp
    fold_stride=None,                # Observations between consecutive test set starts. If None, equals steps (no overlap)
    window_size=None,                # Observations needed for autoregressive predictors (set automatically by forecaster)
    differentiation=None,            # Differencing order to extend last_window (set automatically by forecaster)
    refit=False,                     # Whether to refit forecaster each fold: True, False, or int (refit every n folds)
    fixed_train_size=True,           # If True, training size is fixed; if False, expands each fold
    gap=0,                           # Observations between end of training and start of test set
    skip_folds=None,                 # Folds to skip: int (every n-th fold) or list of fold indexes to skip
    allow_incomplete_fold=True,      # Whether to allow last fold with fewer observations than steps
    return_all_indexes=False,        # Whether to return all indexes or only start/end of each fold
    verbose=True                     # Whether to print information about generated folds
)

# View the folds
folds = cv.split(X=data, as_pandas=True)
print(folds)
```

**Key behaviors:**
- If `fold_stride == steps`: test sets are back-to-back without overlap
- If `fold_stride < steps`: test sets overlap (multiple forecasts for same observations)
- If `fold_stride > steps`: gaps between consecutive test sets

### OneStepAheadFold

Class for one-step-ahead forecasting validation. Faster than `TimeSeriesFold` as it doesn't require recursive predictions.

```python
from skforecast.model_selection import OneStepAheadFold

cv = OneStepAheadFold(
    initial_train_size=100,          # (required) Number of observations for initial training. Can be int, str (date), or pd.Timestamp
    window_size=None,                # Observations needed for autoregressive predictors (set automatically by forecaster)
    differentiation=None,            # Differencing order to extend last_window (set automatically by forecaster)
    return_all_indexes=False,        # Whether to return all indexes or only start/end of each fold
    verbose=True                     # Whether to print information about generated folds
)

# View the fold
fold = cv.split(X=data, as_pandas=True)
print(fold)
```

**When to use:**
- `TimeSeriesFold`: When you need to evaluate multi-step forecasting performance (realistic backtesting)
- `OneStepAheadFold`: When you need fast hyperparameter tuning (less accurate but much faster)

## Hyperparameter Tuning

Three search strategies are available: `grid_search_forecaster`, `random_search_forecaster`, and `bayesian_search_forecaster` (Optuna-based). All accept `TimeSeriesFold` or `OneStepAheadFold`. Multi-series variants: `grid_search_forecaster_multiseries`, `random_search_forecaster_multiseries`, `bayesian_search_forecaster_multiseries`.

```python
from skforecast.model_selection import bayesian_search_forecaster, TimeSeriesFold

cv = TimeSeriesFold(steps=12, initial_train_size=len(data) - 100, refit=False)

# Bayesian Search: lags can be included in search_space
def search_space(trial):
    return {
        'lags': trial.suggest_categorical('lags', [3, 5, [1, 2, 3, 20]]),
        'n_estimators': trial.suggest_int('n_estimators', 50, 200),
        'max_depth': trial.suggest_int('max_depth', 3, 15),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True)
    }

results, study = bayesian_search_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    search_space=search_space,
    metric='mean_absolute_error',
    n_trials=50,
    random_state=123,
    return_best=True,
    n_jobs='auto',
    show_progress=True
)
# Access the best trial with study.best_trial
```

## Statistical Models (ARIMA, ETS, ARAR)

Statistical models are wrapped by `ForecasterStats` for a sklearn-compatible interface. Available models: `Arima`, `Sarimax`, `Ets`, `Arar`.
Autocorrelation utilities are also available in `skforecast.stats`: `acf`, `pacf`, and `calculate_lag_autocorrelation`.

```python
from skforecast.recursive import ForecasterStats
from skforecast.stats import Arima, Ets, Sarimax, Arar
from skforecast.stats import acf, pacf, calculate_lag_autocorrelation

# ARIMA model (order=(p,d,q), seasonal_order=(P,D,Q), m=seasonal_period)
forecaster = ForecasterStats(estimator=Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12))
forecaster.fit(y=data['target'])
predictions = forecaster.predict(steps=10)

# Auto ARIMA (automatic order selection) - set order=None and seasonal_order=None
forecaster = ForecasterStats(estimator=Arima(order=None, seasonal_order=None, m=12))

# ETS model (model string: 1st=Error, 2nd=Trend, 3rd=Seasonal; A=Add, M=Mult, N=None, Z=Auto)
forecaster = ForecasterStats(estimator=Ets(m=12, model='AAA'))
```

## Foundation Models (Zero-Shot)

Pre-trained time series foundation models that forecast without task-specific training. Each model requires its own backend library installed separately (`chronos-forecasting`, `timesfm`, `uni2ts`, `tabicl`, `tabpfn-time-series`, `tfc-t0`). Models are downloaded from HuggingFace on first use.

`FoundationModel` is the low-level interface; `ForecasterFoundation` wraps it to integrate with the rest of the skforecast ecosystem (backtesting, model selection, uniform `predict` / `predict_interval` / `predict_quantiles` API).

```python
from skforecast.foundation import FoundationModel, ForecasterFoundation

# Zero-shot single-series forecasting with Chronos-2
model = FoundationModel(
    model_id='autogluon/chronos-2-small',
    context_length=2048,
    device_map='auto',
)
forecaster = ForecasterFoundation(estimator=model)
forecaster.fit(series=data['target'])  # Only stores context; no training
predictions = forecaster.predict(steps=24)  # Long-format: columns ['level', 'pred']

# Prediction intervals (native quantile output, no bootstrapping needed)
predictions = forecaster.predict_interval(steps=24, interval=[0.1, 0.9])
predictions = forecaster.predict_quantiles(steps=24, quantiles=[0.1, 0.5, 0.9])

# Multi-series (global zero-shot model): pass a wide DataFrame
forecaster.fit(series=series_df)
predictions = forecaster.predict(steps=24, levels=['series_1', 'series_2'])
```

Supported adapters (selected automatically from `model_id`):

| Adapter | `model_id` prefix | Exog | Default `context_length` | Quantiles |
|---------|-------------------|------|--------------------------|-----------|
| ChronosAdapter (Amazon) | `autogluon/chronos` | Yes (past & future covariates) | 8192 | Any in `(0, 1)` |
| TimesFMAdapter (Google) | `google/timesfm` | No | 512 | `[0.1, 0.2, ..., 0.9]` |
| MoiraiAdapter (Salesforce) | `Salesforce/moirai` | No | 2048 | `[0.1, 0.2, ..., 0.9]` |
| TabICLAdapter (Soda-INRIA) | `soda-inria/tabicl` | Yes (past & future covariates) | 4096 | Any in `(0, 1)` |
| TabPFNAdapter (Prior Labs) | `priorlabs/tabpfn` | Yes (known-future covariates) | 32768 | Any in `(0, 1)` |
| T0Adapter (The Forecasting Company) | `theforecastingcompany/t0` | Yes (future-known covariates) | 8192 | Any in `(0, 1)` |

Key points:
- `fit()` only stores the last `context_length` observations and metadata. It does **not** train the model.
- The index must have a frequency (`data.asfreq(...)`), same requirement as other skforecast forecasters.
- `predict(..., context=...)` lets you override the stored context (used internally by backtesting).
- Use `backtesting_foundation` (not `backtesting_forecaster`) to evaluate a `ForecasterFoundation`. Refit is always disabled internally, model weights are preserved across folds since there is no per-fold training.

## Feature Selection

Use sklearn selectors (RFECV, SelectFromModel, etc.) to identify relevant lags, window features, and exogenous variables. Multi-series variant: `select_features_multiseries`.

```python
from sklearn.feature_selection import RFECV
from skforecast.feature_selection import select_features

selected_lags, selected_window_features, selected_exog = select_features(
    forecaster=forecaster,
    selector=RFECV(estimator=RandomForestRegressor(), step=1, cv=3),
    y=y_train,
    exog=exog_train,
    select_only=None,              # 'autoreg', 'exog', or None (all features)
    force_inclusion=None,          # Features to always include (list or regex str)
    subsample=0.5,
    random_state=123,
    verbose=True
)
forecaster.set_lags(selected_lags)
```

## Drift Detection

Two drift detection tools for monitoring data distribution changes during deployment.

```python
from skforecast.drift_detection import RangeDriftDetector, PopulationDriftDetector

# RangeDriftDetector: lightweight, checks out-of-range values vs training ranges
detector = RangeDriftDetector()
detector.fit(series=y_train, exog=exog_train)
flag_out_of_range, out_of_range_series, out_of_range_exog = detector.predict(
    last_window=new_data, exog=new_exog
)

# PopulationDriftDetector: statistical tests (KS, Chi-Square, Jensen-Shannon)
detector = PopulationDriftDetector(chunk_size=100, threshold=3, threshold_method='std')
detector.fit(X=reference_data)
results, summary = detector.predict(X=new_data)
```

## Key Classes and Imports

**Note:** In versions prior to 0.14.0, `ForecasterRecursive` was named `ForecasterAutoreg` and `ForecasterRecursiveMultiSeries` was named `ForecasterAutoregMultiSeries`. Always use the current names shown below.

```python
# Forecasters
from skforecast.recursive import ForecasterRecursive
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.recursive import ForecasterRecursiveClassifier
from skforecast.recursive import ForecasterStats
from skforecast.recursive import ForecasterEquivalentDate
from skforecast.direct import ForecasterDirect
from skforecast.direct import ForecasterDirectMultiVariate
from skforecast.deep_learning import ForecasterRnn
from skforecast.deep_learning import create_and_compile_model
from skforecast.foundation import FoundationModel
from skforecast.foundation import ForecasterFoundation

# Model Selection
from skforecast.model_selection import backtesting_forecaster
from skforecast.model_selection import backtesting_forecaster_multiseries
from skforecast.model_selection import backtesting_stats
from skforecast.model_selection import backtesting_foundation
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import grid_search_forecaster_multiseries
from skforecast.model_selection import random_search_forecaster
from skforecast.model_selection import random_search_forecaster_multiseries
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import bayesian_search_forecaster_multiseries
from skforecast.model_selection import grid_search_stats
from skforecast.model_selection import random_search_stats
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import OneStepAheadFold

# Preprocessing
from skforecast.preprocessing import RollingFeatures
from skforecast.preprocessing import RollingFeaturesClassification
from skforecast.preprocessing import TimeSeriesDifferentiator
from skforecast.preprocessing import CalendarFeatures
from skforecast.preprocessing import create_calendar_features
from skforecast.preprocessing import calculate_distance_from_holiday
from skforecast.preprocessing import QuantileBinner
from skforecast.preprocessing import ConformalIntervalCalibrator
# Data reshaping utilities
from skforecast.preprocessing import reshape_series_wide_to_long
from skforecast.preprocessing import reshape_series_long_to_dict
from skforecast.preprocessing import reshape_exog_long_to_dict
from skforecast.preprocessing import reshape_series_exog_dict_to_long

# Feature Selection
from skforecast.feature_selection import select_features
from skforecast.feature_selection import select_features_multiseries

# Datasets
from skforecast.datasets import fetch_dataset
from skforecast.datasets import load_demo_dataset
from skforecast.datasets import show_datasets_info

# Metrics
from skforecast.metrics import mean_absolute_scaled_error
from skforecast.metrics import root_mean_squared_scaled_error
from skforecast.metrics import symmetric_mean_absolute_percentage_error
from skforecast.metrics import add_y_train_argument
from skforecast.metrics import crps_from_predictions
from skforecast.metrics import crps_from_quantiles
from skforecast.metrics import calculate_coverage
from skforecast.metrics import create_mean_pinball_loss

# Statistical models (used with ForecasterStats)
from skforecast.stats import Arima, Ets, Sarimax, Arar
from skforecast.stats import acf, pacf, calculate_lag_autocorrelation

# Drift Detection
from skforecast.drift_detection import RangeDriftDetector
from skforecast.drift_detection import PopulationDriftDetector

# Plotting
from skforecast.plot import plot_residuals
from skforecast.plot import plot_multivariate_time_series_corr
from skforecast.plot import plot_prediction_distribution
from skforecast.plot import plot_prediction_intervals
from skforecast.plot import backtesting_gif_creator
from skforecast.plot import set_dark_theme

# Exceptions and Warnings
from skforecast.exceptions import DataTypeWarning
from skforecast.exceptions import DataTransformationWarning
from skforecast.exceptions import ExogenousInterpretationWarning
from skforecast.exceptions import FeatureOutOfRangeWarning
from skforecast.exceptions import IgnoredArgumentWarning
from skforecast.exceptions import InputTypeWarning
from skforecast.exceptions import LongTrainingWarning
from skforecast.exceptions import MissingExogWarning
from skforecast.exceptions import MissingValuesWarning
from skforecast.exceptions import OneStepAheadValidationWarning
from skforecast.exceptions import ResidualsUsageWarning
from skforecast.exceptions import SaveLoadSkforecastWarning
from skforecast.exceptions import SkforecastVersionWarning
from skforecast.exceptions import UnknownLevelWarning
from skforecast.exceptions import set_warnings_style
from skforecast.exceptions import warn_skforecast_categories
from skforecast.exceptions import runtime_deprecated

```

## Available Datasets

30+ built-in datasets covering multiple frequencies (15min, 30min, hourly, daily, monthly, quarterly). Most commonly used:

```python
from skforecast.datasets import fetch_dataset, show_datasets_info

data = fetch_dataset(name='h2o')             # Monthly, single series
data = fetch_dataset(name='items_sales')     # Daily, 3 series (multi-series)
data = fetch_dataset(name='bike_sharing')    # Hourly, with exogenous variables
data = fetch_dataset(name='store_sales')     # Daily, 50 products x 10 stores
data = fetch_dataset(name='h2o_exog')        # Monthly, with exogenous variables

# See all available datasets
show_datasets_info()
```

## Tips for Best Results

1. **Always set frequency**: Use `data.asfreq('h')` or similar, skforecast requires a DatetimeIndex with frequency
2. **Handle missing values**: Use `dropna_from_series=True` to drop NaN rows, or `dropna_from_series=False` (default) with NaN-tolerant estimators (LightGBM, CatBoost, HistGradientBoosting)
3. **Scale data**: Use `transformer_y` for better model performance
4. **Use backtesting**: Always validate with realistic train/test splits
5. **Consider differentiation**: For non-stationary series, use `differentiation` parameter
6. **Start simple**: Begin with ForecasterRecursive before trying complex models

## Documentation

- Quick Start: https://skforecast.org/latest/quick-start/quick-start-skforecast.html
- User Guides: https://skforecast.org/latest/user_guides/table-of-contents.html
- API Reference: https://skforecast.org/latest/api/forecasterrecursive.html
- Examples: https://skforecast.org/latest/examples/examples_english.html
- Release Notes: https://skforecast.org/latest/releases/releases.html

## Citation

```
Amat Rodrigo, J., & Escobar Ortiz, J. (2026). skforecast (Version 0.23.0) [Computer software]. https://doi.org/10.5281/zenodo.8382787
```

================================================================================
# SKILL: forecasting-single-series
================================================================================

# Forecasting a Single Time Series

## When to Use

Use this workflow when you have **one time series** and want to predict its future values.

- **ForecasterRecursive**: Default choice. Uses its own predictions as inputs for multi-step forecasting. Works with any sklearn-compatible regressor.
- **ForecasterDirect**: Trains one model per forecast step. Better when the relationship between lags and target changes across the horizon.

### Related skills

- **Before**: `choosing-a-forecaster` (decide between Recursive and Direct based on the data)
- **Before**: `autocorrelation-and-lag-selection` (pick the `lags` argument from ACF/PACF analysis)
- **Before**: `feature-engineering` (assemble the rolling, calendar, and exogenous features)
- **After**: `hyperparameter-optimization` (tune the forecaster once a baseline is trained)
- **After**: `prediction-intervals` (add bootstrap or conformal intervals on top of the point forecasts)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| Set the index frequency before fitting | `ValueError: ... must be a DatetimeIndex with frequency` | Call `data = data.asfreq('h')` (or the correct alias) on `y` and `exog` |
| `exog` passed to `predict()` must cover every future step | `exog` length / index error, or the forecast stops short | Slice exog to the full horizon: one row per step, dates matching the forecast |
| Fit with `store_in_sample_residuals=True` before `predict_interval(method='bootstrapping')` | `No in-sample residuals stored` / empty residuals | Refit: `forecaster.fit(y=y_train, store_in_sample_residuals=True)` |
| Split chronologically, never shuffle | Over-optimistic metrics, leakage | Use a time-ordered split (`TimeSeriesFold`); do not random-split |

## Complete Workflow

```python
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from skforecast.recursive import ForecasterRecursive
from skforecast.preprocessing import RollingFeatures
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

# 1. Load and prepare data (MUST have DatetimeIndex with frequency)
data = pd.read_csv('data.csv', index_col='date', parse_dates=True)
data = data.asfreq('h')  # Set frequency — required

# 2. Train/test split (time-based, never random)
end_train = '2023-01-01'
y_train = data.loc[:end_train, 'target']
y_test = data.loc[end_train:, 'target']

# 3. Create forecaster with optional rolling features
rolling_features = RollingFeatures(
    stats=['mean', 'std'],
    window_sizes=24
)

forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=24,
    window_features=rolling_features,
    transformer_y=None,          # e.g., StandardScaler() for scaling
    categorical_features='auto', # Auto-detect and encode non-numeric exog columns
    differentiation=None,        # e.g., 1 for first-order differencing
    dropna_from_series=False,    # True to drop NaN rows; False to keep (NaN-tolerant estimators)
)

# 4. Train
forecaster.fit(y=y_train)

# 5. Predict
predictions = forecaster.predict(steps=10)

# 6. Backtesting (proper evaluation)
cv = TimeSeriesFold(
    steps=10,
    initial_train_size=len(y_train),
    refit=False,
    fixed_train_size=False,
)

metric, predictions_bt = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
)
print(f"MAE: {metric}")

# 7. Prediction intervals
# Default interval is [0.05, 0.95] (90%). Here [0.1, 0.9] creates an 80% interval.
forecaster.fit(y=y_train, store_in_sample_residuals=True)
predictions_interval = forecaster.predict_interval(
    steps=10,
    interval=[0.1, 0.9],        # 80% prediction interval (quantiles, 0-1 scale)
    method='bootstrapping',
    n_boot=500,
)
```

## With Exogenous Variables

```python
# Exogenous variables must cover the forecast horizon during prediction
forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=24,
)
forecaster.fit(y=y_train, exog=exog_train)
predictions = forecaster.predict(steps=10, exog=exog_test)
```

## Using ForecasterDirect

```python
from skforecast.direct import ForecasterDirect

# Must specify `steps` at creation — trains one model per step
forecaster = ForecasterDirect(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=24,
    steps=10,
    categorical_features='auto',  # Auto-detect and encode non-numeric exog columns
)
forecaster.fit(y=y_train, exog=exog_train)
predictions = forecaster.predict(exog=exog_test)
```

## Common Mistakes

1. **Missing frequency on index**: Always call `data.asfreq('h')` (or `'D'`, `'MS'`, etc.).
2. **NaN in data**: Forecasters reject NaN by default. Use `dropna_from_series=True` to drop incomplete rows, or keep `dropna_from_series=False` (default) with NaN-tolerant estimators (LightGBM, CatBoost, HistGradientBoosting, XGBoost hist). Alternatively, impute missing values first.
3. **Exog not covering forecast horizon**: The exogenous DataFrame for `predict()` must have rows for every future step.
4. **Random train/test split**: Time series must be split chronologically, never shuffled.
5. **Forgetting `store_in_sample_residuals=True`**: Required before calling `predict_interval()` with `method='bootstrapping'` on a standalone forecaster. During backtesting, residuals are computed automatically.

================================================================================
# SKILL: forecasting-multiple-series
================================================================================

# Forecasting Multiple Time Series

## When to Use

- **ForecasterRecursiveMultiSeries**: A single global model learns patterns across many series. Each series is predicted independently but the model shares parameters. Best default for multi-series.
- **ForecasterDirectMultiVariate**: Uses values from multiple series as input features to predict one target series. Use when series are strongly correlated and influence each other.

### Related skills

- **Before**: `choosing-a-forecaster` (decide between MultiSeries, MultiVariate, Rnn, or Foundation for multi-series problems)
- **Before**: `autocorrelation-and-lag-selection` (analyse representative series to inform the shared `lags` argument)
- **Before**: `feature-engineering` (build per-series exogenous and rolling features)
- **After**: `hyperparameter-optimization` (tune the global model across series)
- **After**: `prediction-intervals` (add intervals to the multi-series forecasts)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| Use `backtesting_forecaster_multiseries` and the `*_multiseries` search functions | `backtesting_forecaster` / `grid_search_forecaster` raises on a multi-series forecaster | Call the `_multiseries` variant with `series=` instead of `y=` |
| `ForecasterDirectMultiVariate` defaults to `transformer_series=StandardScaler()` | Series are scaled unexpectedly (other forecasters default to `None`) | Pass `transformer_series=None` explicitly if you do not want scaling |
| `exog` format must match the `series` format (both wide, or both dict) | Index / format mismatch during fit or predict | Convert exog to the same layout as `series` before fitting |
| Regressors with native categorical support need `encoding='ordinal_category'` | Categoricals silently encoded as plain ordinals, degrading the model | Set `encoding='ordinal_category'` for LightGBM / CatBoost / XGBoost / HistGBR |

## Data Formats

ForecasterRecursiveMultiSeries accepts three input formats:

1. **Wide DataFrame** — columns are series, index is datetime:
   ```python
   #            series_1  series_2  series_3
   # 2020-01-01     1.0       2.5       3.1
   # 2020-01-02     1.2       2.3       3.4
   ```

2. **Dictionary** — `{series_id: pd.Series}`:
   ```python
   {'series_1': pd.Series([1.0, 1.2, ...]), 'series_2': pd.Series([2.5, 2.3, ...])}
   ```

> **Note:** Long-format DataFrames are not directly accepted. Use `reshape_series_long_to_dict()` to convert long format to a dictionary first (see Data Reshaping Utilities below).

## Complete Workflow

```python
import pandas as pd
from lightgbm import LGBMRegressor
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import backtesting_forecaster_multiseries, TimeSeriesFold

# 1. Load data as wide DataFrame (columns = series)
series = pd.read_csv('data.csv', index_col='date', parse_dates=True)
series = series.asfreq('D')

# 2. Create forecaster
forecaster = ForecasterRecursiveMultiSeries(
    estimator=LGBMRegressor(n_estimators=200, random_state=123),
    lags=24,
    encoding='ordinal',       # 'ordinal', 'ordinal_category', 'onehot', or None
    transformer_series=None,  # Apply same transformer to all series
    # Per-series options (dict):
    # transformer_series={'series_1': StandardScaler(), 'series_2': MinMaxScaler()},
    # weight_func={'series_1': custom_weights_fn, '_default': None},
    # differentiation={'series_1': 1, 'series_2': None},
    categorical_features='auto',  # Auto-detect and encode non-numeric exog columns
    differentiation=None,
    dropna_from_series=False,     # True to drop NaN rows; False to keep (NaN-tolerant estimators)
)

# 3. Train
forecaster.fit(series=series)

# 4. Predict all series
predictions = forecaster.predict(steps=10)

# 5. Predict specific series only
predictions = forecaster.predict(steps=10, levels=['series_1', 'series_2'])

# 6. Backtesting (multi-series)
cv = TimeSeriesFold(
    steps=10,
    initial_train_size=len(series) - 100,
    refit=False,
)

metric, predictions_bt = backtesting_forecaster_multiseries(
    forecaster=forecaster,
    series=series,
    cv=cv,
    metric='mean_absolute_error',
    levels=None,  # Evaluate all series
)
print(metric)  # Shows per-series and aggregated metrics
```

## With Exogenous Variables

```python
# Exog can also be wide DataFrame, long DataFrame, or dict
forecaster.fit(series=series, exog=exog_df)
predictions = forecaster.predict(steps=10, exog=exog_test)
```

## ForecasterDirectMultiVariate

```python
from skforecast.direct import ForecasterDirectMultiVariate

# Predicts ONE target series using lags from ALL series as features
# Note: transformer_series defaults to StandardScaler() (unlike other forecasters)
forecaster = ForecasterDirectMultiVariate(
    level='target_series',    # Name of the series to predict
    steps=10,
    estimator=LGBMRegressor(n_estimators=100, random_state=123),
    lags=24,                  # Or dict: {'series_a': 12, 'series_b': 24}
    transformer_series=StandardScaler(),  # Default — set None to disable scaling
    categorical_features='auto',  # Auto-detect and encode non-numeric exog columns
    dropna_from_series=False,     # True to drop NaN rows; False to keep (NaN-tolerant estimators)
)
forecaster.fit(series=series_df)
predictions = forecaster.predict()
```

## Data Reshaping Utilities

```python
from skforecast.preprocessing import (
    reshape_series_wide_to_long,
    reshape_series_long_to_dict,
    reshape_exog_long_to_dict,
    reshape_series_exog_dict_to_long,
)

# Wide → Long
series_long = reshape_series_wide_to_long(series_wide)

# Long → Dict (freq is required)
series_dict = reshape_series_long_to_dict(series_long, freq='D')
exog_dict = reshape_exog_long_to_dict(exog_long, freq='D')
```

## Common Mistakes

1. **Mismatched series lengths**: ForecasterRecursiveMultiSeries handles different-length series if `dropna_from_series=True`.
2. **Wrong encoding for categorical regressor**: Use `encoding='ordinal_category'` with regressors that natively handle categoricals (LightGBM, CatBoost).
3. **Exog format mismatch**: Exog format (wide/dict) must match the series format.
4. **Forgetting `levels` parameter**: By default `predict()` forecasts all series. Use `levels` to limit predictions.
5. **Unexpected scaling in ForecasterDirectMultiVariate**: `transformer_series` defaults to `StandardScaler()`, unlike other forecasters that default to `None`. Set `transformer_series=None` explicitly if you don't want automatic scaling.

================================================================================
# SKILL: statistical-models
================================================================================

# Statistical Models (ARIMA, ETS, SARIMAX, ARAR)

## References

See [references/model-parameters.md](references/model-parameters.md) for
complete constructor signatures of all statistical models (Arima, Sarimax,
Ets, Arar), the Ets model string format, Auto-ARIMA parameters, seasonal_order
differences between Arima and Sarimax, and grid search param_grid examples.

## When to Use

Use statistical models when:
- The series is short (< 200 observations)
- Interpretability is important (ARIMA coefficients, ETS components)
- You need built-in prediction intervals without residual bootstrapping
- As a baseline to compare against ML models

### Related skills

- **Before**: `autocorrelation-and-lag-selection` (read ACF/PACF to identify ARIMA orders `p`, `d`, `q` before fitting)
- **After**: `prediction-intervals` (`ForecasterStats` provides built-in parametric intervals via the `interval_method` argument)
- **After**: `hyperparameter-optimization` (tune ARIMA `order` / `seasonal_order` via grid search)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| Backtest and tune with `backtesting_stats` and `grid_search_stats`, not the ML variants | `backtesting_forecaster` / `grid_search_forecaster` raises on `ForecasterStats` | Call `backtesting_stats` / `grid_search_stats` |
| `Arima` takes a 3-tuple `seasonal_order=(P, D, Q)` plus `m`; `Sarimax` takes a 4-tuple `(P, D, Q, m)` | Wrong model order or `TypeError` | Use `Arima(order=(p,d,q), seasonal_order=(P,D,Q), m=12)` |
| `Ets` uses `model='AAA'` + `m`, not `error=` / `trend=` / `seasonal=` | Deprecated-argument error | Use `Ets(model='AAA', m=12)` (A/M/N/Z per position) |
| Seasonal models require `m` | Seasonality silently ignored | Pass `m=<seasonal period>` to `Arima` / `Ets` |

## Available Models

| Model | Class | Description |
|-------|-------|-------------|
| **ARIMA** | `Arima` | AutoRegressive Integrated Moving Average |
| **Auto-ARIMA** | `Arima(order=None)` | Automatic order selection |
| **SARIMAX** | `Sarimax` | ARIMA with exogenous variables (seasonal) |
| **ETS** | `Ets` | Exponential Smoothing (Error-Trend-Seasonal) |
| **ARAR** | `Arar` | Autoregressive model with memory shortening |

## Complete Workflow: ARIMA

```python
import pandas as pd
from skforecast.recursive import ForecasterStats
from skforecast.stats import Arima
from skforecast.model_selection import backtesting_stats, TimeSeriesFold

# 1. Load data
data = pd.read_csv('data.csv', index_col='date', parse_dates=True)
data = data.asfreq('MS')  # Monthly Start frequency

# 2. Manual ARIMA: specify order and seasonal_order
arima_model = Arima(
    order=(1, 1, 1),              # (p, d, q)
    seasonal_order=(1, 1, 1),     # (P, D, Q)
    m=12,                         # Seasonal period
)
forecaster = ForecasterStats(estimator=arima_model)
forecaster.fit(y=data['target'])
predictions = forecaster.predict(steps=12)

# 3. Prediction intervals (all stat models support this natively,
#    no bootstrapping needed). Accepts both `interval` and `alpha`.
predictions_interval = forecaster.predict_interval(
    steps=12,
    interval=[0.1, 0.9],  # quantiles (0-1). Or use alpha=0.2 for 80% interval
)
```

## Auto-ARIMA (Automatic Order Selection)

```python
# Set order=None and seasonal_order=None to enable automatic order selection
auto_arima = Arima(order=None, seasonal_order=None, m=12)
forecaster = ForecasterStats(estimator=auto_arima)
forecaster.fit(y=data['target'])

# Check selected order
print(forecaster.estimator.best_params_['order'])
print(forecaster.estimator.best_params_['seasonal_order'])

predictions = forecaster.predict(steps=12)
```

## ETS (Exponential Smoothing)

```python
from skforecast.stats import Ets

# Model string: 1st=Error, 2nd=Trend, 3rd=Seasonal
# A=Additive, M=Multiplicative, N=None, Z=Auto-select
ets_model = Ets(model='AAA', m=12)
forecaster = ForecasterStats(estimator=ets_model)
forecaster.fit(y=data['target'])
predictions = forecaster.predict(steps=12)
```

## Auto-ETS (Automatic Model Selection)

```python
# Use model='ZZZ' (or model=None) to let ETS automatically select
# the best Error, Trend, and Seasonal components
auto_ets = Ets(model='ZZZ', m=12)
forecaster = ForecasterStats(estimator=auto_ets)
forecaster.fit(y=data['target'])

# Check the selected model configuration
print(forecaster.estimator.best_params_)

predictions = forecaster.predict(steps=12)
```

## SARIMAX (with Exogenous Variables)

```python
from skforecast.stats import Sarimax

sarimax_model = Sarimax(
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),  # (P, D, Q, seasonal_period)
)
forecaster = ForecasterStats(estimator=sarimax_model)
forecaster.fit(y=data['target'], exog=exog_train)

# For prediction, exog must cover the forecast horizon
predictions = forecaster.predict(steps=12, exog=exog_test)
```

## ARAR

```python
from skforecast.stats import Arar

arar_model = Arar()
forecaster = ForecasterStats(estimator=arar_model)
forecaster.fit(y=data['target'])
predictions = forecaster.predict(steps=12)
```

## Backtesting Statistical Models

```python
cv = TimeSeriesFold(
    steps=12,
    initial_train_size=len(data) - 60,
    refit=False,
)

metric, predictions_bt = backtesting_stats(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
    freeze_params=True,  # Params from first fit reused in refits (avoids re-running auto selection)
)
# If freeze_params=False, auto selection runs independently each fold and output
# includes an extra 'estimator_params' column with the parameters selected per fold.
```

## Multiple Models Simultaneously

```python
# ForecasterStats accepts a list of models — fits each independently
from skforecast.stats import Arima, Ets

models = [
    Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12),
    Ets(model='AAA', m=12),
]
forecaster = ForecasterStats(estimator=models)
forecaster.fit(y=data['target'])

# predict returns DataFrame with one column per model
predictions = forecaster.predict(steps=12)
```

## Common Mistakes

1. **Using deprecated `Ets(error=, trend=, seasonal=)` syntax**: Use `Ets(model='AAA', m=12)` with a model string instead.
2. **Forgetting `m` parameter**: ARIMA and ETS seasonal models require `m` (seasonal period).
3. **Not using `backtesting_stats`**: Use `backtesting_stats()` for statistical models, NOT `backtesting_forecaster()`.
4. **Using grid_search_forecaster for stats**: Use `grid_search_stats()` or `random_search_stats()` instead.
5. **Passing `seasonal_order=(1,1,1,12)` to `Arima`**: `Arima` uses a 3-tuple `seasonal_order=(P,D,Q)` plus a separate `m=12` parameter. The 4-tuple `(P,D,Q,s)` format is only for `Sarimax`.

---
### Reference: model-parameters

# Statistical Models — Parameter Reference

## Constructor Comparison

| Parameter | `Arima` | `Sarimax` | `Ets` | `Arar` |
|-----------|:-------:|:---------:|:-----:|:------:|
| `order` | ✓ `(1,0,0)` | ✓ `(1,0,0)` | — | — |
| `seasonal_order` | ✓ `(0,0,0)` 3-tuple | ✓ `(0,0,0,0)` **4-tuple** | — | — |
| `m` | ✓ `1` | — (in seasonal_order) | ✓ `1` | — |
| `model` | — | — | ✓ `'ZZZ'` | — |
| Auto-selection | `order=None` | — | `model='ZZZ'` or `None` | — |

> **Critical difference:** `Arima` uses 3-tuple `seasonal_order=(P,D,Q)` + separate `m`.
> `Sarimax` uses 4-tuple `seasonal_order=(P,D,Q,s)` with season period embedded.

## Arima

```python
Arima(
    # --- Manual ARIMA (set order explicitly) ---
    order=(1, 0, 0),                  # tuple(p, d, q) | None → auto ARIMA
    seasonal_order=(0, 0, 0),         # tuple(P, D, Q) — 3 elements, NOT 4
    m=1,                              # int, seasonal period (1 = no seasonality)
    include_mean=True,                # bool
    transform_pars=True,              # bool
    method='CSS-ML',                  # str, fitting method
    n_cond=None,                      # int | None
    SSinit='Gardner1980',             # str
    optim_method='BFGS',              # str
    optim_kwargs=None,                # dict | None
    kappa=1e6,                        # float

    # --- Auto ARIMA (only used when order=None) ---
    max_p=5,                          # int
    max_q=5,                          # int
    max_P=2,                          # int
    max_Q=2,                          # int
    max_order=5,                      # int
    max_d=2,                          # int
    max_D=1,                          # int
    start_p=2,                        # int
    start_q=2,                        # int
    start_P=1,                        # int
    start_Q=1,                        # int
    stationary=False,                 # bool
    seasonal=True,                    # bool
    ic='aicc',                        # 'aic' | 'aicc' | 'bic'
    stepwise=True,                    # bool
    nmodels=94,                       # int
    trace=False,                      # bool
    approximation=None,               # bool | None
    truncate=None,                    # int | None
    test='kpss',                      # str
    test_kwargs=None,                 # dict | None
    seasonal_test='seas',             # str
    seasonal_test_kwargs=None,        # dict | None
    allowdrift=True,                  # bool
    allowmean=True,                   # bool

    # --- Box-Cox transformation ---
    lambda_bc=None,                   # float | str | None
    biasadj=False,                    # bool
)
```

### Arima usage modes

| Mode | How to activate | Description |
|------|----------------|-------------|
| **Manual** | `order=(p,d,q)` | User specifies exact order |
| **Auto** | `order=None` | Automatic order selection via stepwise algorithm |

After fitting (auto mode), selected order available in:
- `forecaster.estimator.best_params_['order']`
- `forecaster.estimator.best_params_['seasonal_order']`

## Sarimax

```python
Sarimax(
    order=(1, 0, 0),                  # tuple(p, d, q)
    seasonal_order=(0, 0, 0, 0),      # tuple(P, D, Q, s) — 4 elements
    trend=None,                       # str | None
    measurement_error=False,          # bool
    time_varying_regression=False,    # bool
    mle_regression=True,              # bool
    simple_differencing=False,        # bool
    enforce_stationarity=True,        # bool
    enforce_invertibility=True,       # bool
    hamilton_representation=False,    # bool
    concentrate_scale=False,          # bool
    trend_offset=1,                   # int
    use_exact_diffuse=False,          # bool
    dates=None,
    freq=None,
    missing='none',                   # str
    validate_specification=True,      # bool
    method='lbfgs',                   # str, fitting method
    maxiter=50,                       # int
    start_params=None,                # np.ndarray | None
    disp=False,                       # bool
    sm_init_kwargs={},                # dict, extra kwargs for statsmodels init
    sm_fit_kwargs={},                 # dict, extra kwargs for statsmodels fit
    sm_predict_kwargs={},             # dict, extra kwargs for statsmodels predict
)
```

> **Note:** `Sarimax` is the only model that supports exogenous variables natively.

## Ets

```python
Ets(
    m=1,                              # int, seasonal period (1 = no seasonality)
    model='ZZZ',                      # str | None, model specification string
    damped=None,                      # bool | None
    alpha=None,                       # float | None, smoothing level
    beta=None,                        # float | None, smoothing trend
    gamma=None,                       # float | None, smoothing seasonal
    phi=None,                         # float | None, damping parameter
    lambda_param=None,                # float | None, Box-Cox lambda
    lambda_auto=False,                # bool
    bias_adjust=True,                 # bool
    bounds='both',                    # str
    seasonal=True,                    # bool
    trend=None,                       # bool | None
    ic='aicc',                        # 'aic' | 'aicc' | 'bic'
    allow_multiplicative=True,        # bool
    allow_multiplicative_trend=False, # bool
)
```

### Ets model string format

The `model` parameter is a 3-character string: `Error`, `Trend`, `Seasonal`.

| Character | Meaning | Position |
|-----------|---------|----------|
| `A` | Additive | Any |
| `M` | Multiplicative | Any |
| `N` | None (not present) | Trend or Seasonal |
| `Z` | Auto-select | Any |

| Example | Error | Trend | Seasonal | Description |
|---------|-------|-------|----------|-------------|
| `'AAN'` | Additive | Additive | None | Simple exponential smoothing with trend |
| `'AAA'` | Additive | Additive | Additive | Holt-Winters additive |
| `'MAM'` | Multiplicative | Additive | Multiplicative | Holt-Winters multiplicative seasonal |
| `'ANA'` | Additive | None | Additive | Seasonal without trend |
| `'ZZZ'` | Auto | Auto | Auto | Auto-ETS (selects best) |

After fitting (auto mode), selected configuration available in:
- `forecaster.estimator.best_params_`

## Arar

```python
Arar(
    max_ar_depth=None,                # int | None
    max_lag=None,                     # int | None
    safe=True,                        # bool
)
```

> `Arar` is the simplest model — no configuration of order or seasonality needed.
> The algorithm automatically determines the best AR structure via memory shortening.

## ForecasterStats Constructor

```python
ForecasterStats(
    estimator=None,           # Arima | Sarimax | Ets | Arar | list of these
    transformer_y=None,       # sklearn transformer for target variable
    transformer_exog=None,    # sklearn transformer | ColumnTransformer for exog
    forecaster_id=None,       # str | int, optional identifier
)
```

### Multiple models simultaneously

```python
# Pass a list to fit multiple models at once
models = [
    Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12),
    Ets(model='AAA', m=12),
    Arar(),
]
forecaster = ForecasterStats(estimator=models)
forecaster.fit(y=data)

# predict() returns DataFrame with one column per model
predictions = forecaster.predict(steps=12)
```

## Grid Search param_grid Examples

```python
# Arima manual orders
param_grid = {
    'order': [(1, 0, 0), (1, 1, 0), (1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0), (1, 1, 1)],
    'm': [12],
}

# Ets models
param_grid = {
    'model': ['AAA', 'ANA', 'MAM', 'MNM', 'ZZZ'],
    'm': [12],
}

# Sarimax
param_grid = {
    'order': [(1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0, 12), (1, 1, 1, 12)],
}
```

## Backtesting with freeze_params

```python
backtesting_stats(
    forecaster=forecaster,
    y=data,
    cv=cv,
    metric='mean_absolute_error',
    freeze_params=True,     # Default. Reuse params from first fit in all folds.
    # freeze_params=False,  # Re-run auto selection each fold (slow but more realistic).
    #                       # Output includes extra 'estimator_params' column.
)
```

## Exogenous Variables Support

| Model | Exog in fit/predict | Notes |
|-------|:--:|------|
| `Arima` | — | No exog support |
| `Sarimax` | ✓ | Full exog support via statsmodels SARIMAX |
| `Ets` | — | No exog support |
| `Arar` | — | No exog support |

For `Sarimax` with exog, pass `exog` to both `fit()` and `predict()`:

```python
forecaster = ForecasterStats(estimator=Sarimax(order=(1,1,1), seasonal_order=(1,1,1,12)))
forecaster.fit(y=y_train, exog=exog_train)
predictions = forecaster.predict(steps=12, exog=exog_test)
```

================================================================================
# SKILL: metric-selection
================================================================================

# Metric Selection

## When to Use This Skill

Use this skill when the user needs help choosing an evaluation metric, comparing
metrics, understanding trade-offs between error measures, or configuring the
`metric` parameter in `backtesting_forecaster`, `bayesian_search_forecaster`,
or any other model selection function.

### Related skills

- **After**: `choosing-a-forecaster` (the forecaster type determines which metrics apply)
- **Before**: `hyperparameter-optimization` (the chosen metric drives the search objective)
- **Before**: `prediction-intervals` (probabilistic metrics evaluate interval quality)

## Quick Recommendations

If unsure, start here:

**Point forecast metrics** (pass as `metric=` in backtesting/search):

| Task | Default metric | Why |
|------|---------------|-----|
| General purpose | `'mean_absolute_error'` | Interpretable, robust to outliers, works at any scale |
| Compare across series | `'mean_absolute_scaled_error'` | Scale-independent — the only fair comparison across differently-scaled series |
| Classification | `'balanced_accuracy_score'` | Handles class imbalance common in time series classification |

**Probabilistic metrics** (computed post-hoc on interval/quantile predictions):

| Task | Metric | Why |
|------|--------|-----|
| Interval calibration | `calculate_coverage` | Check if actual coverage matches the nominal level |
| Interval quality (overall) | `crps_from_predictions` | Evaluates sharpness and calibration together |
| Quantile quality (foundation models) | `crps_from_quantiles` | Proper scoring rule for quantile predictions |
| Single-quantile optimization | `create_mean_pinball_loss(alpha)` | Can be passed as `metric=` to optimize for a specific quantile |

## Step 1 — What Are You Evaluating?

```
What is your forecaster producing?
│
├─► Point forecasts (single value per step)
│   └─► Go to Step 2
│
├─► Prediction intervals (lower_bound, upper_bound)
│   └─► Go to Step 3
│
├─► Quantile predictions (multiple quantile levels)
│   └─► Go to Step 3
│
├─► Class labels (ForecasterRecursiveClassifier)
│   └─► Go to Step 4
│
└─► Multi-series with aggregation needs
    └─► Go to Step 5 (then return to Step 2 or 3 for the base metric)
```

## Step 2 — Point Forecast Metrics

### Decision Table

| Criterion | Recommended metric | Avoid |
|-----------|-------------------|-------|
| **General purpose** | MAE (`'mean_absolute_error'`) | — |
| **Penalize large errors more** | MSE (`'mean_squared_error'`) | — |
| **Robust to outliers** | MAE or MedAE (`'median_absolute_error'`) | MSE (inflated by outliers) |
| **Need percentage interpretation** | SMAPE (`'symmetric_mean_absolute_percentage_error'`) | MAPE if data has zeros |
| **Data contains zeros or near-zero values** | MAE, MASE, RMSSE | MAPE (divides by y_true → infinite) |
| **Compare across different-scale series** | MASE (`'mean_absolute_scaled_error'`) or RMSSE (`'root_mean_squared_scaled_error'`) | MAE, MSE (scale-dependent) |
| **Target is always positive** | MSLE (`'mean_squared_log_error'`) | — (use when relative errors matter more than absolute) |
| **Interpretable baseline comparison** | MASE (value < 1 means better than naive forecast) | — |

### Metric Properties

| Metric | Scale-independent | Robust to outliers | Handles zeros | Requires `y_train` | Range |
|--------|:-----------------:|:------------------:|:-------------:|:-------------------:|-------|
| MAE | — | ✓ | ✓ | — | [0, ∞) |
| MSE | — | — | ✓ | — | [0, ∞) |
| MedAE | — | ✓✓ | ✓ | — | [0, ∞) |
| MAPE | ✓ | — | — | — | [0, ∞) |
| SMAPE | ✓ | — | ✓ | — | [0, 200] % |
| MSLE | — | — | ✓ (if ≥ 0) | — | [0, ∞) |
| MASE | ✓ | ✓ | ✓ | ✓ | [0, ∞) |
| RMSSE | ✓ | — | ✓ | ✓ | [0, ∞) |

### Using Metrics in Code

#### Direct computation (train/test split)

For a simple train/test evaluation, import the metric and call it directly.
Sklearn metrics are imported from `sklearn.metrics`; skforecast-specific metrics
(`mean_absolute_scaled_error`, `root_mean_squared_scaled_error`,
`symmetric_mean_absolute_percentage_error`) from `skforecast.metrics`.

```python
from sklearn.metrics import mean_squared_error
from skforecast.metrics import mean_absolute_scaled_error

forecaster.fit(y=data_train)
predictions = forecaster.predict(steps=36)

# Sklearn metrics: func(y_true, y_pred)
error_mse = mean_squared_error(y_true=data_test, y_pred=predictions)

# Skforecast metrics that need y_train: func(y_true, y_pred, y_train)
error_mase = mean_absolute_scaled_error(
    y_true=data_test, y_pred=predictions, y_train=data_train
)
```

#### Inside backtesting

When using backtesting or hyperparameter search, pass metrics as strings —
no import needed, and `y_train` is handled automatically:

```python
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

cv = TimeSeriesFold(steps=12, initial_train_size=len(y_train), refit=False)

# Single metric (pass as string)
metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
)

# Multiple metrics at once
metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric=['mean_absolute_error', 'mean_absolute_scaled_error'],
)
# metric is a DataFrame with one column per metric
```

**MASE and RMSSE** require training data to compute the naive forecast baseline.
When passed as a string or callable, skforecast handles `y_train` automatically —
it detects whether the function signature includes a `y_train` parameter and
provides the training data if so.

```python
# Both options work identically in backtesting — y_train is handled internally
metric = 'mean_absolute_scaled_error'            # Option 1: as string
metric = mean_absolute_scaled_error              # Option 2: as callable (import from skforecast.metrics)
```

## Step 3 — Probabilistic Forecast Metrics

Use these metrics when evaluating prediction intervals or quantile forecasts.

| Metric | Evaluates | Input | Use case |
|--------|-----------|-------|----------|
| `calculate_coverage` | Calibration | y_true, lower_bound, upper_bound | Check if actual coverage matches nominal level |
| `crps_from_predictions` | Calibration + sharpness | y_true (scalar), y_pred (array of bootstrap samples) | Evaluate bootstrapped interval quality |
| `crps_from_quantiles` | Calibration + sharpness | y_true (scalar), pred_quantiles, quantile_levels | Evaluate quantile predictions (foundation models) |
| `create_mean_pinball_loss(alpha)` | Single-quantile accuracy | y_true, y_pred (at quantile alpha) | Evaluate a specific quantile forecast |

### Coverage

Coverage measures the proportion of true values that fall within the predicted
interval. **Target: match the nominal level** (e.g., 90% interval should have ~90% coverage).

```python
from skforecast.metrics import calculate_coverage

coverage = calculate_coverage(
    y_true=y_test,
    lower_bound=predictions['lower_bound'],
    upper_bound=predictions['upper_bound'],
)
# Ideal: coverage ≈ 0.90 for a 90% interval
# coverage >> 0.90 → intervals too wide (not sharp)
# coverage << 0.90 → intervals too narrow (miscalibrated)
```

### CRPS (Continuous Ranked Probability Score)

CRPS is a proper scoring rule that rewards both calibration and sharpness.
Lower is better.

```python
from skforecast.metrics import crps_from_predictions, crps_from_quantiles
import numpy as np

# From bootstrap samples (e.g., predict_bootstrapping output)
crps = crps_from_predictions(
    y_true=100.0,
    y_pred=np.array([98.2, 101.5, 99.8, 102.1, 97.5])  # Bootstrap predictions
)

# From quantile predictions (e.g., foundation model output)
crps = crps_from_quantiles(
    y_true=100.0,
    pred_quantiles=np.array([90.0, 95.0, 100.5, 105.0, 110.0]),
    quantile_levels=np.array([0.1, 0.25, 0.5, 0.75, 0.9]),
)
```

### Pinball Loss (Quantile Loss)

Evaluates a single quantile forecast. Use `create_mean_pinball_loss(alpha)` to
create a metric function for a specific quantile level.

```python
from skforecast.metrics import create_mean_pinball_loss

# Create metric for the 90th percentile
pinball_90 = create_mean_pinball_loss(alpha=0.9)

# Use in backtesting (as callable)
metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric=pinball_90,
)
```

## Step 4 — Classification Metrics

For `ForecasterRecursiveClassifier` only. These metrics evaluate predicted
class labels, not continuous values.

| Metric | Best for |
|--------|----------|
| `'accuracy_score'` | Balanced classes (equal importance per class) |
| `'balanced_accuracy_score'` | Imbalanced classes (common in time series — e.g., rare events) |
| `'f1_score'` | When both precision and recall matter |
| `'precision_score'` | When false positives are costly |
| `'recall_score'` | When false negatives are costly (e.g., missing a spike) |

```python
from skforecast.recursive import ForecasterRecursiveClassifier
from sklearn.ensemble import RandomForestClassifier
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

forecaster = ForecasterRecursiveClassifier(
    estimator=RandomForestClassifier(n_estimators=100, random_state=123),
    lags=24,
)

# y contains class labels (e.g., 'low', 'medium', 'high'); encoding is handled internally
cv = TimeSeriesFold(steps=10, initial_train_size=len(y) - 100, refit=False)

metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=y,
    cv=cv,
    metric='balanced_accuracy_score',
)
```

## Step 5 — Multi-Series Metric Aggregation

When predicting multiple series, skforecast computes per-series metrics and
can aggregate them into a single score. The aggregation options are:

| Aggregation | Formula | Use when |
|-------------|---------|----------|
| `'average'` | Arithmetic mean of per-series metrics | All series equally important |
| `'weighted_average'` | Weighted by number of predicted values per level | Series with more observations contribute more |
| `'pooling'` | Metric computed on all predictions pooled together | Want a single global error regardless of series identity |

### In backtesting

Use `add_aggregated_metric=True` (default) to include all three aggregations:

```python
from skforecast.model_selection import backtesting_forecaster_multiseries

metric, predictions = backtesting_forecaster_multiseries(
    forecaster=forecaster_multi,
    series=series_df,
    cv=cv,
    metric='mean_absolute_error',
    add_aggregated_metric=True,  # Default: includes average, weighted_average, pooling
)
# metric DataFrame has rows for each level + aggregated rows
```

### In hyperparameter search

Use the `aggregate_metric` parameter to select which aggregations to report:

```python
from skforecast.model_selection import bayesian_search_forecaster_multiseries

results, study = bayesian_search_forecaster_multiseries(
    forecaster=forecaster_multi,
    series=series_df,
    cv=cv,
    search_space=search_space,
    metric='mean_absolute_error',
    aggregate_metric=['weighted_average', 'average', 'pooling'],  # Select which to report
)
```

**Tip:** Use `'pooling'` when you care about overall accuracy regardless
of which series contributes the error. Use `'average'` when all series are
equally important regardless of their length.

## Custom Metrics

Any callable with signature `func(y_true, y_pred)` can be passed as a metric.
If your metric also needs training data, include a `y_train` parameter:

```python
# Custom metric without y_train
def mean_bias_error(y_true, y_pred):
    return np.mean(y_pred - y_true)

# Custom metric WITH y_train (detected automatically)
def relative_mase(y_true, y_pred, y_train):
    naive_mae = np.mean(np.abs(np.diff(y_train)))
    return np.mean(np.abs(y_true - y_pred)) / naive_mae

# Both work directly in backtesting
metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric=[mean_bias_error, relative_mase],
)
```

Skforecast automatically wraps callables via `add_y_train_argument`. If the function
already has a `y_train` parameter, it will receive the training data. If not, the
argument is silently ignored.

## Common Mistakes

| Mistake | Problem | Fix |
|---------|---------|-----|
| Using MAPE with data containing zeros | Division by zero → infinite error | Use SMAPE, MAE, or MASE instead |
| Using MSE when data has outliers | A few large errors dominate the metric | Use MAE or MedAE for robustness |
| Comparing MAE across differently-scaled series | MAE = 10 on sales vs MAE = 10 on temperature means nothing | Use MASE or RMSSE for cross-series comparison |
| Interpreting high coverage as "good" | 99% coverage with a 90% interval means intervals are too wide | Target coverage ≈ nominal level; combine with CRPS for sharpness |
| Using only point metrics for probabilistic forecasts | Ignores uncertainty quality entirely | Add `calculate_coverage` and/or CRPS alongside point metrics |
| Passing MASE as callable without understanding y_train | Works fine — skforecast detects the `y_train` parameter automatically | Just pass the callable or string; no extra work needed |

## Metric Compatibility Reference

See [references/metric-compatibility.md](references/metric-compatibility.md) for the
complete matrix of all 18 metrics with their properties, compatible forecasters,
and usage guidance.

---
### Reference: metric-compatibility

# Metric Compatibility Matrix

Complete reference for all metrics available in skforecast: properties, compatibility, and usage guidance.

## Point Forecast Metrics (Regression)

| Metric | String name | Source | Requires `y_train` | Scale-independent | Robust to outliers | Handles zeros | Range | Interpretation |
|--------|-------------|--------|:-------------------:|:-----------------:|:------------------:|:-------------:|-------|----------------|
| Mean Absolute Error | `'mean_absolute_error'` | sklearn | — | — | ✓ | ✓ | [0, ∞) | Average absolute deviation from truth |
| Mean Squared Error | `'mean_squared_error'` | sklearn | — | — | — | ✓ | [0, ∞) | Average squared deviation; penalizes large errors |
| Median Absolute Error | `'median_absolute_error'` | sklearn | — | — | ✓✓ | ✓ | [0, ∞) | Median of absolute deviations; very robust |
| Mean Absolute Percentage Error | `'mean_absolute_percentage_error'` | sklearn | — | ✓ | — | — | [0, ∞) | Percentage error; undefined when y_true = 0 |
| Mean Squared Log Error | `'mean_squared_log_error'` | sklearn | — | — | — | ✓* | [0, ∞) | Log-scale MSE; penalizes under-predictions more |
| Symmetric MAPE | `'symmetric_mean_absolute_percentage_error'` | skforecast | — | ✓ | — | ✓ | [0, 200] % | Symmetric percentage; avoids MAPE's asymmetry |
| Mean Absolute Scaled Error | `'mean_absolute_scaled_error'` | skforecast | ✓ | ✓ | ✓ | ✓ | [0, ∞) | Ratio vs naive forecast; < 1 = better than naive |
| Root Mean Squared Scaled Error | `'root_mean_squared_scaled_error'` | skforecast | ✓ | ✓ | — | ✓ | [0, ∞) | RMSE ratio vs naive; < 1 = better than naive |

*MSLE requires y_true ≥ 0 and y_pred ≥ 0 (typically used with positive data).

## Probabilistic Forecast Metrics

| Metric | Function | Input signature | Use with | Measures |
|--------|----------|-----------------|----------|----------|
| Coverage | `calculate_coverage` | `(y_true, lower_bound, upper_bound)` | Any forecaster with intervals | Calibration (proportion of truths in interval) |
| CRPS (bootstrap) | `crps_from_predictions` | `(y_true: float, y_pred: ndarray)` | Bootstrapped predictions | Calibration + sharpness (single observation) |
| CRPS (quantiles) | `crps_from_quantiles` | `(y_true: float, pred_quantiles: ndarray, quantile_levels: ndarray)` | Foundation models, quantile predictions | Calibration + sharpness (single observation) |
| Pinball Loss | `create_mean_pinball_loss(alpha)` | Returns `func(y_true, y_pred)` | Any quantile forecast at level alpha | Single-quantile accuracy |

## Classification Metrics

| Metric | String name | Source | Best for |
|--------|-------------|--------|----------|
| Accuracy | `'accuracy_score'` | sklearn | Balanced classes |
| Balanced Accuracy | `'balanced_accuracy_score'` | sklearn | Imbalanced classes (recommended default) |
| F1 Score | `'f1_score'` | sklearn | Balance of precision and recall |
| Precision | `'precision_score'` | sklearn | Minimizing false positives |
| Recall | `'recall_score'` | sklearn | Minimizing false negatives |

## Forecaster × Metric Type Compatibility

| Forecaster | Point metrics | Probabilistic metrics | Classification metrics |
|------------|:-------------:|:---------------------:|:---------------------:|
| `ForecasterRecursive` | ✓ | ✓ (bootstrapping, conformal) | — |
| `ForecasterDirect` | ✓ | ✓ (bootstrapping, conformal) | — |
| `ForecasterRecursiveMultiSeries` | ✓ | ✓ (bootstrapping, conformal) | — |
| `ForecasterDirectMultiVariate` | ✓ | ✓ (bootstrapping, conformal) | — |
| `ForecasterRnn` | ✓ | ✓ (conformal) | — |
| `ForecasterStats` | ✓ | ✓ (native parametric) | — |
| `ForecasterFoundation` | ✓ | ✓ (native quantiles) | — |
| `ForecasterEquivalentDate` | ✓ | ✓ (conformal) | — |
| `ForecasterRecursiveClassifier` | — | — | ✓ |

## When to Use / When to Avoid

| Metric | Use when | Avoid when |
|--------|----------|------------|
| MAE | General purpose; robust default; easy to explain | Need to heavily penalize large errors |
| MSE | Large errors should be penalized disproportionally | Data has outliers (a few bad predictions dominate) |
| MedAE | Data is noisy with frequent outliers | Need sensitivity to all errors (median ignores tails) |
| MAPE | Need percentage interpretation; all values are non-zero | Data contains zeros or near-zero values |
| SMAPE | Need percentage interpretation; data may have zeros | Need strict mathematical properties (SMAPE has known biases) |
| MSLE | Positive data where relative errors matter more than absolute | Data has zeros or negative values |
| MASE | Comparing across series with different scales; need baseline reference | Very short training series (naive denominator unreliable) |
| RMSSE | Same as MASE but want to penalize large errors more | Same as MASE |
| Coverage | Evaluating prediction interval calibration | Want to measure interval width/sharpness (use CRPS) |
| CRPS | Comprehensive interval/quantile evaluation (calibration + sharpness) | Only have point forecasts |
| Pinball loss | Evaluating a specific quantile (e.g., 90th percentile for risk) | Need overall distributional assessment (use CRPS) |
| Balanced accuracy | Classification with imbalanced classes | Perfectly balanced classes (accuracy is simpler) |
| F1 | Classification where both false positives and negatives matter | Need per-class detail (use precision/recall separately) |

## Multi-Series Aggregation Options

Used in `backtesting_forecaster_multiseries` (via `add_aggregated_metric=True`) and
multi-series search functions (via `aggregate_metric` parameter).

| Aggregation | Behavior | Best for |
|-------------|----------|----------|
| `'average'` | Arithmetic mean of per-series metric values | All series equally important |
| `'weighted_average'` | Weighted by number of predicted values per level | Series with more observations contribute more |
| `'pooling'` | Metric computed on all predictions concatenated | Global accuracy regardless of series identity |

In backtesting, set `add_aggregated_metric=True` (default) to include all three.
In search functions, use `aggregate_metric` to select which to report (default: all three).

================================================================================
# SKILL: backtesting-configuration
================================================================================

# Backtesting Configuration

## When to Use

Use this skill to translate a deployment scenario (retraining cadence, forecast
horizon, data budget, ingestion delay) into `TimeSeriesFold` parameters.
`TimeSeriesFold` is the cross-validation strategy passed via the `cv` argument to
`backtesting_forecaster` (and its multi-series / stats variants) and to the
hyperparameter search functions.

```python
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

cv = TimeSeriesFold(
    steps=7,                       # forecast horizon
    initial_train_size=365,        # first training window (required)
    refit=False,                   # train once (default)
    fixed_train_size=True,         # rolling window (default)
    gap=0,
)

metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
)
```

For fast hyperparameter tuning where multi-step realism is not required, use
`OneStepAheadFold` instead: it validates one step ahead (no recursive
prediction), so it is much faster but less representative of multi-step
performance. Use `TimeSeriesFold` for realistic multi-step backtesting.

### Related skills

- **Before**: `forecasting-single-series` / `forecasting-multiple-series` (have a fitted forecaster before backtesting)
- **With**: `metric-selection` (choose the metric(s) backtesting reports)
- **With**: `hyperparameter-optimization` (the same `cv` object drives the search functions)
- **After**: `prediction-intervals` (add `interval=` / `interval_method=` to backtest uncertainty)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| `initial_train_size` must be provided (the API default is `None`, which only works when reusing an already-fitted forecaster) | `ValueError` / no training in the first fold | Pass an int, date string, or Timestamp, e.g. `initial_train_size=len(y) - 100` |
| `initial_train_size` does not accept a float fraction | `ValueError`: must be int, date string, Timestamp, or None | Convert fractions to an int, e.g. `int(len(data) * 0.7)` |
| Configuration must yield at least 2 folds | Single-fold or empty backtest | Ensure `initial_train_size + gap + 2 * steps <= n_observations` |
| `refit=True` retrains every fold (slow path); the default is `refit=False` | Backtest much slower than expected | Use `refit=False` or an int cadence (e.g. `refit=7`) unless per-fold retraining is required |
| `gap`, `steps`, and `fold_stride` count observations, not calendar time | Off-by-frequency delays | Convert the delay to the series frequency (e.g. 2 days at daily freq = `gap=2`) |

## TimeSeriesFold Parameters

| Parameter | Type | Default | Description |
|-----------|------|---------|-------------|
| `initial_train_size` | int, str, pd.Timestamp | None (required) | Observations for initial training. Int = count, str/Timestamp = last training date. A common starting point is `int(len(data) * 0.7)`. |
| `refit` | bool, int | False | Refit every fold (True), train once (False, default), or every n folds (int). |
| `fixed_train_size` | bool | True | Rolling fixed window (True, default) vs expanding window (False). |
| `gap` | int | 0 | Observations between training end and test start. |
| `fold_stride` | int, None | None (= steps) | Distance between consecutive test set starts. |
| `skip_folds` | int, list, None | None | Skip folds to reduce compute. Int = keep every n-th. |
| `allow_incomplete_fold` | bool | True | Allow final fold with fewer observations than steps. |

## Constraints

- Must produce **at least 2 folds**: `initial_train_size + gap + 2 * steps <= n_observations`
- `initial_train_size` must be large enough for the model to learn patterns (at minimum 2× the lag order for ML models, or 2× steps for statistical models)
- `gap` simulates real-world delay between data availability and forecast usage
- When `fixed_train_size=True`, the window rolls forward (oldest data discarded). Use for concept drift or when old data is less relevant.
- When `fixed_train_size=False`, the window expands (all history retained). Use when more data always helps.

## Business Scenario Mapping

### Retraining frequency

| Scenario | Configuration |
|----------|--------------|
| Retrain every time new data arrives | `refit=True` |
| Retrain weekly (with daily forecasts, steps=1) | `refit=7` |
| Retrain monthly (with daily forecasts, steps=7) | `refit=4` (every 4 folds × 7 steps ≈ monthly) |
| Never retrain / train once, evaluate across time | `refit=False` (`fixed_train_size` has no effect without refit) |

### Data freshness vs volume

| Scenario | Configuration |
|----------|--------------|
| Recent data more relevant (concept drift) | `fixed_train_size=True` |
| All historical data valuable | `fixed_train_size=False` (expanding) |
| Limited compute budget | `refit=False` or `skip_folds=3` (keep every 3rd fold) |

### Deployment gap

| Scenario | Configuration |
|----------|--------------|
| Real-time predictions (no delay) | `gap=0` |
| 1-day delay between data collection and forecast | `gap=1` (if freq=daily) |
| Forecast must be ready before weekend (2-day gap) | `gap=2` |

### Initial training size

| Scenario | Configuration |
|----------|--------------|
| Default (balanced) | `int(len(data) * 0.7)` |
| Maximize evaluation coverage | Minimum viable: `2 * max_lag` or `2 * steps` |
| Maximize training data | `n_observations - gap - 2 * steps` (minimum 2 folds) |
| Start from specific date | `initial_train_size="2023-01-01"` |
| Conservative (large training set) | `int(len(data) * 0.8)` |

### Fold stride (test set overlap)

| Scenario | Configuration |
|----------|--------------|
| Non-overlapping evaluation (default) | `fold_stride=None` (equals steps) |
| Sliding window evaluation (1-step shift) | `fold_stride=1` |
| Sparse evaluation (save compute) | `fold_stride=steps * 2` |

## Examples

### "I retrain my model every Monday and forecast the next 7 days"
```
refit = True
fixed_train_size = False  # Keep all history
gap = 0  # No delay
fold_stride = None  # Non-overlapping weeks
```

### "I want to simulate deploying once and seeing how the model degrades"
```
refit = False  # Train once; fixed_train_size has no effect without refit
gap = 0
```

### "There's a 2-day lag between data ingestion and when forecasts are needed"
```
gap = 2
refit = True
```

### "I have limited compute — evaluate every other week"
```
refit = True
skip_folds = 2  # Keep every 2nd fold
```

### "I want maximum evaluation coverage with a 12-step horizon"
```
initial_train_size = <minimum viable>  # 2 * max_lag
refit = False  # Faster; trains once, fixed_train_size has no effect
fold_stride = 1  # Sliding window (many overlapping folds)
```

================================================================================
# SKILL: hyperparameter-optimization
================================================================================

# Hyperparameter Optimization

## References

See [references/search-parameters.md](references/search-parameters.md) for
the complete parameter comparison across all 9 search functions, function
routing by forecaster type, and `lags_grid` / `search_space` / `param_grid`
usage details.

## When to Use

Use hyperparameter search after establishing a baseline forecaster to improve prediction accuracy. Skforecast supports three strategies:

| Strategy | When to Use | Speed |
|----------|-------------|-------|
| **Bayesian Search** | **Recommended default.** Smart exploration via Optuna | Fastest to converge |
| **Random Search** | Large parameter space, limited compute budget | Medium |
| **Grid Search** | Small parameter space, exhaustive exploration | Slowest |

### Related skills

- **Before**: `autocorrelation-and-lag-selection` (narrow the `lags` search space to a statistically informed candidate set)
- **Before**: `feature-selection` (run the search on a reduced feature set to make it tractable)
- **After**: `prediction-intervals` (add uncertainty quantification once the configuration is fixed)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| The search refits the forecaster in place with the best params when `return_best=True` (the default) | With `return_best=False`, the forecaster keeps its pre-search params | Rely on the default, or refit with the best params from the results table |
| Use the `*_multiseries` / `*_stats` search variant matching the forecaster type | Search function raises on the wrong forecaster type | Call e.g. `bayesian_search_forecaster_multiseries` / `grid_search_stats` |
| Include `lags` in the Bayesian `search_space()` | Suboptimal search; the highest-impact parameter stays fixed | Add `trial.suggest_categorical('lags', [...])` to the search space |

## Bayesian Search (Recommended)

Always prefer Bayesian search as the default strategy. It uses Optuna to intelligently explore the search space.

```python
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import bayesian_search_forecaster, TimeSeriesFold
from lightgbm import LGBMRegressor

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(random_state=123),
    lags=24,
)

cv = TimeSeriesFold(
    steps=12,
    initial_train_size=len(data) - 100,
    refit=False,
)

# Define search space as a function — lags CAN be included here
def search_space(trial):
    return {
        'lags': trial.suggest_categorical('lags', [12, 24, [1, 2, 3, 23, 24]]),
        'n_estimators': trial.suggest_int('n_estimators', 50, 500),
        'max_depth': trial.suggest_int('max_depth', 3, 15),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 10.0, log=True),
    }

# n_trials=20 is the default. Increase for better results (50-200 recommended).
results, study = bayesian_search_forecaster(
    forecaster=forecaster,
    y=data['target'],
    exog=exog,
    cv=cv,
    search_space=search_space,
    metric='mean_absolute_error',
    n_trials=20,
    random_state=123,
    return_best=True,        # Automatically updates forecaster with best params
    n_jobs='auto',
    show_progress=True,
    output_file='search_results.csv',  # Save results incrementally
)
# results is a DataFrame sorted by metric (best first)
# study is the full Optuna Study; access the best trial with study.best_trial
```

## Grid Search

```python
from skforecast.model_selection import grid_search_forecaster

# Different lag configurations to try
lags_grid = [3, 10, 24, [1, 2, 3, 23, 24]]

param_grid = {
    'n_estimators': [50, 100, 200],
    'max_depth': [5, 10, 15],
    'learning_rate': [0.01, 0.1],
}

results = grid_search_forecaster(
    forecaster=forecaster,
    y=data['target'],
    exog=exog,
    cv=cv,
    lags_grid=lags_grid,
    param_grid=param_grid,
    metric='mean_absolute_error',
    return_best=True,
    n_jobs='auto',
    show_progress=True,
)
```

## Random Search

```python
from skforecast.model_selection import random_search_forecaster

# Note: uses param_distributions (not param_grid) and n_iter
param_distributions = {
    'n_estimators': [50, 100, 200, 500],
    'max_depth': [3, 5, 10, 15],
    'learning_rate': [0.01, 0.05, 0.1, 0.3],
}

results = random_search_forecaster(
    forecaster=forecaster,
    y=data['target'],
    exog=exog,
    cv=cv,
    lags_grid=lags_grid,
    param_distributions=param_distributions,
    n_iter=10,               # Number of random parameter combinations to try
    random_state=123,
    metric='mean_absolute_error',
    return_best=True,
    n_jobs='auto',
    show_progress=True,
)
```

## Multi-Series Search

```python
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import bayesian_search_forecaster_multiseries

forecaster = ForecasterRecursiveMultiSeries(
    estimator=LGBMRegressor(random_state=123),
    lags=24,
    encoding='ordinal',
)

cv = TimeSeriesFold(
    steps=12,
    initial_train_size=len(series) - 100,
    refit=False,
)

results, study = bayesian_search_forecaster_multiseries(
    forecaster=forecaster,
    series=series,
    exog=exog,
    cv=cv,
    search_space=search_space,
    metric='mean_absolute_error',
    aggregate_metric=['weighted_average', 'average', 'pooling'],  # Default
    levels=None,             # None = evaluate all series; or list of series names
    n_trials=20,
    return_best=True,
    n_jobs='auto',
    show_progress=True,
)
# Access the best trial with study.best_trial
```

## Statistical Models Search

```python
from skforecast.recursive import ForecasterStats
from skforecast.stats import Arima
from skforecast.model_selection import grid_search_stats

forecaster = ForecasterStats(estimator=Arima(order=(1, 1, 1)))

param_grid = {
    'order': [(1, 0, 0), (1, 1, 0), (1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0), (1, 1, 1)],
    'm': [12],
}

results = grid_search_stats(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    param_grid=param_grid,
    metric='mean_absolute_error',
    return_best=True,
)
```

## Fast Tuning with OneStepAheadFold

```python
from skforecast.model_selection import OneStepAheadFold

# Much faster than TimeSeriesFold — no recursive predictions needed
cv_fast = OneStepAheadFold(
    initial_train_size=len(data) - 100,
)

results, study = bayesian_search_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv_fast,
    search_space=search_space,
    metric='mean_absolute_error',
    n_trials=100,
    return_best=True,
)
# Access the best trial with study.best_trial
```

## Common Mistakes

1. **Not setting `return_best=True`**: The forecaster is not updated with the best parameters unless this is True.
2. **Too few trials in Bayesian search**: Start with at least 20-50 trials for meaningful exploration.
3. **Using TimeSeriesFold for initial tuning**: Use `OneStepAheadFold` first for fast screening, then validate the top candidates with `TimeSeriesFold`.
4. **Forgetting to include lags in search space**: For Bayesian search, lags can be included in `search_space()` — this is often the most impactful parameter.

---
### Reference: search-parameters

# Hyperparameter Search — Parameter Reference

## Function Routing

| Forecaster | Grid Search | Random Search | Bayesian Search |
|------------|-------------|---------------|-----------------|
| ForecasterRecursive | `grid_search_forecaster` | `random_search_forecaster` | `bayesian_search_forecaster` |
| ForecasterDirect | `grid_search_forecaster` | `random_search_forecaster` | `bayesian_search_forecaster` |
| ForecasterRecursiveMultiSeries | `grid_search_forecaster_multiseries` | `random_search_forecaster_multiseries` | `bayesian_search_forecaster_multiseries` |
| ForecasterDirectMultiVariate | `grid_search_forecaster_multiseries` | `random_search_forecaster_multiseries` | `bayesian_search_forecaster_multiseries` |
| ForecasterRnn | `grid_search_forecaster_multiseries` | `random_search_forecaster_multiseries` | `bayesian_search_forecaster_multiseries` |
| ForecasterStats | `grid_search_stats` | `random_search_stats` | N/A |
| ForecasterEquivalentDate | N/A | N/A | N/A |
| ForecasterRecursiveClassifier | `grid_search_forecaster` | `random_search_forecaster` | `bayesian_search_forecaster` |

## Parameter Comparison Across Search Functions

### Single-Series Functions

| Parameter | `grid_search_forecaster` | `random_search_forecaster` | `bayesian_search_forecaster` |
|-----------|:-:|:-:|:-:|
| `forecaster` | ✓ | ✓ | ✓ |
| `y` | ✓ | ✓ | ✓ |
| `cv` | TimeSeriesFold \| OneStepAheadFold | TimeSeriesFold \| OneStepAheadFold | TimeSeriesFold \| OneStepAheadFold |
| `param_grid` | ✓ | — | — |
| `param_distributions` | — | ✓ | — |
| `search_space` | — | — | ✓ (Callable) |
| `metric` | ✓ | ✓ | ✓ |
| `exog` | ✓ | ✓ | ✓ |
| `lags_grid` | ✓ | ✓ | — (included in `search_space`) |
| `n_iter` | — | ✓ (default: 10) | — |
| `n_trials` | — | — | ✓ (default: 20) |
| `random_state` | — | ✓ (default: 123) | ✓ (default: 123) |
| `return_best` | ✓ (default: True) | ✓ (default: True) | ✓ (default: True) |
| `n_jobs` | ✓ (default: 'auto') | ✓ (default: 'auto') | ✓ (default: 'auto') |
| `verbose` | ✓ | ✓ | ✓ |
| `show_progress` | ✓ | ✓ | ✓ |
| `suppress_warnings` | ✓ | ✓ | ✓ |
| `output_file` | ✓ | ✓ | ✓ |
| `kwargs_create_study` | — | — | ✓ (default: None) |
| `kwargs_study_optimize` | — | — | ✓ (default: None) |
| **Returns** | `pd.DataFrame` | `pd.DataFrame` | `tuple[pd.DataFrame, object]` |

### Multi-Series Functions (additional parameters)

These functions have all the parameters above plus:

| Parameter | grid | random | bayesian |
|-----------|:----:|:------:|:--------:|
| `series` (replaces `y`) | ✓ | ✓ | ✓ |
| `aggregate_metric` | ✓ | ✓ | ✓ |
| `levels` | ✓ | ✓ | ✓ |

Default `aggregate_metric = ['weighted_average', 'average', 'pooling']`

### Stats Functions (limited parameters)

| Parameter | `grid_search_stats` | `random_search_stats` |
|-----------|:---:|:---:|
| `forecaster` | ✓ | ✓ |
| `y` | ✓ | ✓ |
| `cv` | **TimeSeriesFold only** | **TimeSeriesFold only** |
| `param_grid` / `param_distributions` | ✓ | ✓ |
| `metric` | ✓ | ✓ |
| `exog` | ✓ | ✓ |
| `lags_grid` | — | — |
| `n_iter` | — | ✓ (default: 10) |
| `random_state` | — | ✓ (default: 123) |
| `return_best` | ✓ | ✓ |
| `n_jobs` | ✓ | ✓ |
| **Returns** | `pd.DataFrame` | `pd.DataFrame` |

> **Note:** Stats search does NOT support `OneStepAheadFold`, `lags_grid`,
> or Bayesian search.

## How `lags_grid` Works

For grid and random search, `lags_grid` is a list of lag configurations to try:

```python
# List format — each element is a configuration
lags_grid = [3, 10, 24, [1, 2, 3, 23, 24]]
# Tries: lags=3, lags=10, lags=24, lags=[1,2,3,23,24]

# Dict format — keys become labels in the results
lags_grid = {
    'short': 3,
    'medium': 12,
    'long': 24,
    'custom': [1, 2, 3, 23, 24],
}
```

For Bayesian search, lags are included in the `search_space` function:

```python
def search_space(trial):
    return {
        'lags': trial.suggest_categorical('lags', [3, 12, 24, [1, 2, 3, 23, 24]]),
        'n_estimators': trial.suggest_int('n_estimators', 50, 500),
    }
```

## How `param_grid` vs `param_distributions` vs `search_space` Work

### Grid Search: `param_grid`

All combinations are evaluated (Cartesian product):

```python
param_grid = {
    'n_estimators': [50, 100, 200],
    'max_depth': [5, 10],
}
# Evaluates: 3 × 2 = 6 combinations
```

### Random Search: `param_distributions`

Random sample of `n_iter` combinations:

```python
param_distributions = {
    'n_estimators': [50, 100, 200, 500],
    'max_depth': [3, 5, 10, 15],
    'learning_rate': [0.01, 0.05, 0.1, 0.3],
}
# Evaluates n_iter=10 random combinations (default)
```

### Bayesian Search: `search_space`

Optuna trial function with suggest methods:

```python
def search_space(trial):
    return {
        'lags': trial.suggest_categorical('lags', [12, 24]),
        'n_estimators': trial.suggest_int('n_estimators', 50, 500),
        'max_depth': trial.suggest_int('max_depth', 3, 15),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 10.0, log=True),
    }
# Optuna methods: suggest_int, suggest_float, suggest_categorical
# Evaluates n_trials=20 trials (default), guided by TPE sampler
```

## Stats Model param_grid

Parameters in `param_grid` for stats models are passed to the model constructor:

```python
# Arima
param_grid = {
    'order': [(1, 0, 0), (1, 1, 0), (1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0), (1, 1, 1)],
    'm': [12],
}

# Ets
param_grid = {
    'model': ['AAA', 'ANA', 'MAM', 'ZZZ'],
    'm': [12],
}
```

## Optuna kwargs

```python
# Advanced: customize Optuna study
results, study = bayesian_search_forecaster(
    ...,
    kwargs_create_study={
        'sampler': optuna.samplers.TPESampler(seed=123),
        'direction': 'minimize',
    },
    kwargs_study_optimize={
        'timeout': 600,  # seconds
        'gc_after_trial': True,
    },
)
# Access the best trial with study.best_trial
```

## Return Values

| Function | Returns | Study object |
|----------|---------|:-:|
| `grid_search_*` | `pd.DataFrame` sorted by metric | — |
| `random_search_*` | `pd.DataFrame` sorted by metric | — |
| `bayesian_search_*` | `tuple[pd.DataFrame, optuna Study]` | ✓ |
| `*_stats` | `pd.DataFrame` sorted by metric | — |

When `return_best=True`, the forecaster is automatically updated with the
best parameters found. The results DataFrame always has rows sorted by
metric (best first). Access the best Optuna trial with `study.best_trial`.

================================================================================
# SKILL: prediction-intervals
================================================================================

# Prediction Intervals

## When to Use

Use prediction intervals to quantify forecast uncertainty. Skforecast offers three methods:

| Method | Forecasters | Description |
|--------|-------------|-------------|
| **Bootstrapping** | Recursive, Direct, RecursiveMultiSeries, DirectMultiVariate | Resample from training residuals |
| **Conformal** | All ML forecasters (incl. RNN, EquivalentDate) | Distribution-free intervals via conformal prediction |
| **Built-in** | ForecasterStats (ARIMA, ETS) | Parametric intervals from the statistical model |

### Choosing a method

- **Conformal**: distribution-free, fast, and gives a single symmetric coverage
  guarantee. The default for multi-series forecasters and the only option for
  `ForecasterRnn` / `ForecasterEquivalentDate`. Prefer it when you need a quick,
  well-calibrated interval and do not need the full predictive distribution.
- **Bootstrapping**: resamples residuals to build the full predictive
  distribution, so it also powers `predict_quantiles()` and `predict_dist()`.
  Prefer it when you need arbitrary quantiles, a fitted distribution, or
  asymmetric intervals. Needs enough stored residuals to be reliable.
- **Built-in**: parametric intervals from the statistical model itself
  (`ForecasterStats`). No residual management required.

### Related skills

- **Before**: `forecasting-single-series` / `forecasting-multiple-series` (have a fitted forecaster before adding intervals)
- **Before**: `hyperparameter-optimization` (interval calibration assumes a tuned point forecaster)
- **After**: `drift-detection` (monitor whether residual assumptions still hold once the model is in production)

## Intervals Are Quantiles (Changed in 0.23.0)

Since skforecast 0.23.0, `interval` is expressed as **quantiles in `[0, 1]`**, not
percentiles in `[0, 100]`. For example, an 80% interval is `interval=[0.1, 0.9]`
(the default is `[0.05, 0.95]` = 90%).

- `interval` also accepts a single `float` as nominal coverage: `interval=0.95`
  is equivalent to `interval=[0.025, 0.975]`.
- Passing percentiles (e.g. `[10, 90]`) still works but emits a `FutureWarning`
  and will be removed in skforecast 0.25.0.
- Mixing scales (e.g. `[0.1, 90]`) raises a `ValueError` (ambiguous scale).
- `predict_quantiles(quantiles=...)` was already on the 0-1 scale and is unchanged.

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| Fit with `store_in_sample_residuals=True` before `predict_interval(method='bootstrapping')` | `No in-sample residuals stored` | Refit: `forecaster.fit(y=y_train, store_in_sample_residuals=True)` |
| Multi-series forecasters default to `method='conformal'`, not `'bootstrapping'` | Bootstrapping silently not applied on `ForecasterRecursiveMultiSeries` / `ForecasterDirectMultiVariate` | Pass `method='bootstrapping'` explicitly when that is what you want |
| ML forecasters accept `interval=`; only `ForecasterStats` accepts `alpha=` | `TypeError` / unexpected argument | Use `interval=[lower, upper]` for ML forecasters |
| `interval` must be quantiles in `[0, 1]`, not percentiles | `FutureWarning` (percentiles deprecated, removed in 0.25.0) | Use `interval=[0.1, 0.9]` instead of `interval=[10, 90]` |
| Do not mix quantile and percentile scales in one `interval` | `ValueError`: scale is ambiguous | Use a single scale, e.g. `interval=[0.05, 0.95]` |

## Bootstrapping Method

```python
from skforecast.recursive import ForecasterRecursive
from sklearn.ensemble import RandomForestRegressor

forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=24,
)

# IMPORTANT: store_in_sample_residuals=True is required for bootstrapping
forecaster.fit(y=y_train, store_in_sample_residuals=True)

predictions = forecaster.predict_interval(
    steps=10,
    interval=[0.1, 0.9],             # Quantiles → 80% interval (default is [0.05, 0.95] = 90%)
    method='bootstrapping',
    n_boot=250,                      # Number of bootstrap samples (default)
    use_in_sample_residuals=True,    # Use training residuals
    use_binned_residuals=True,       # Better calibration: residuals binned by prediction level
    random_state=123,
)
# Returns: DataFrame with columns [pred, lower_bound, upper_bound]
```

## Conformal Prediction

```python
# IMPORTANT: store_in_sample_residuals=True is also required for conformal
forecaster.fit(y=y_train, store_in_sample_residuals=True)

predictions = forecaster.predict_interval(
    steps=10,
    interval=[0.1, 0.9],
    method='conformal',
    use_in_sample_residuals=True,    # Uses in_sample_residuals_ stored during fit
    use_binned_residuals=True,
)
```

## Statistical Model Intervals

```python
from skforecast.recursive import ForecasterStats
from skforecast.stats import Arima

forecaster = ForecasterStats(
    estimator=Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12)
)
forecaster.fit(y=y_train)

# Uses parametric intervals from statsmodels — different interface
predictions = forecaster.predict_interval(
    steps=12,
    interval=[0.1, 0.9],       # Quantiles → 80% interval. Equivalently, alpha=0.2
)
```

## Foundation Model Intervals

`ForecasterFoundation` returns intervals and quantiles directly from the underlying foundation model's native quantile output — no bootstrapping or conformal calibration required.

```python
from skforecast.foundation import FoundationModel, ForecasterFoundation

forecaster = ForecasterFoundation(
    estimator=FoundationModel(model_id='autogluon/chronos-2-small')
)
forecaster.fit(series=y_train)

predictions = forecaster.predict_interval(steps=24, interval=[0.1, 0.9])
predictions = forecaster.predict_quantiles(
    steps=24, quantiles=[0.1, 0.5, 0.9]
)
```

TimesFM 2.5 and Moirai-2 restrict quantiles to `[0.1, 0.2, …, 0.9]`; Chronos-2, TabICL, TabPFN-TS and TFC-T0 accept any quantile in `(0, 1)`. See the `foundation-forecasting` skill for details.

## During Backtesting

```python
from skforecast.model_selection import backtesting_forecaster, TimeSeriesFold

cv = TimeSeriesFold(
    steps=10,
    initial_train_size=len(y_train),
    refit=False,
)

metric, predictions = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
    interval=[0.1, 0.9],
    interval_method='bootstrapping',   # or 'conformal'
    n_boot=250,
    use_in_sample_residuals=True,
    use_binned_residuals=True,
)
# predictions has columns: pred, lower_bound, upper_bound
```

## Multi-Series Intervals

```python
from skforecast.recursive import ForecasterRecursiveMultiSeries

# NOTE: default method for multi-series is 'conformal', not 'bootstrapping'
predictions = forecaster_multi.predict_interval(
    steps=10,
    levels=['series_1', 'series_2'],
    interval=[0.1, 0.9],
    method='conformal',
)
```

## Out-of-Sample Residuals (Better Calibration)

```python
# For better interval calibration, use out-of-sample residuals
# First, compute them via backtesting
metric, predictions_bt = backtesting_forecaster(
    forecaster=forecaster,
    y=data['target'],
    cv=cv,
    metric='mean_absolute_error',
)

# Set out-of-sample residuals on the forecaster
forecaster.set_out_sample_residuals(
    y_true=data['target'].loc[predictions_bt.index],
    y_pred=predictions_bt['pred'],
)

# Now use them for intervals
predictions = forecaster.predict_interval(
    steps=10,
    interval=[0.1, 0.9],
    method='bootstrapping',
    use_in_sample_residuals=False,  # Use out-of-sample residuals
)
```

## Evaluating Interval Quality

```python
from skforecast.metrics import calculate_coverage

coverage = calculate_coverage(
    y_true=y_test,
    lower_bound=predictions['lower_bound'],
    upper_bound=predictions['upper_bound'],
)
print(f"Coverage: {coverage:.2%}")  # Should be close to 0.80 for [0.1, 0.9] interval
```

## Beyond Intervals

For bootstrapping-capable ML forecasters (Recursive, Direct, RecursiveMultiSeries,
DirectMultiVariate), two richer probabilistic outputs build on the same residual
resampling:

```python
# Arbitrary quantiles — returns one column per quantile (q_0.05, q_0.5, q_0.95)
predictions = forecaster.predict_quantiles(
    steps=10,
    quantiles=[0.05, 0.5, 0.95],
    n_boot=250,
)

# Fit a scipy distribution to the bootstrapped predictions
from scipy.stats import norm

predictions = forecaster.predict_dist(
    steps=10,
    distribution=norm,   # any scipy.stats distribution
    n_boot=250,
)
```

See [references/interval-compatibility.md](references/interval-compatibility.md) for the full method/forecaster support matrix.

## Common Mistakes

1. **Passing percentiles instead of quantiles**: Since 0.23.0, `interval` is on the 0-1 scale. Use `interval=[0.1, 0.9]`, not `[10, 90]` (percentiles are deprecated and emit a `FutureWarning`, removed in 0.25.0). Mixing scales, e.g. `[0.1, 90]`, raises a `ValueError`.
2. **Forgetting `store_in_sample_residuals=True`**: Required in `fit()` before using `predict_interval(method='bootstrapping')`.
3. **Wrong default method for multi-series**: `ForecasterRecursiveMultiSeries` and `ForecasterDirectMultiVariate` default to `method='conformal'`, not `'bootstrapping'`.
4. **Mixing `alpha` and `interval`**: `ForecasterStats` supports both `alpha` (e.g., `alpha=0.05` for 95% interval) and `interval=[lo, hi]` (quantiles). ML forecasters only support `interval`.
5. **Not evaluating coverage**: Always check if actual coverage matches nominal interval width.

## References

See [references/interval-compatibility.md](references/interval-compatibility.md) for the complete compatibility matrix: which forecaster supports which method, parameter differences, binned residuals support, and common error patterns.

---
### Reference: interval-compatibility

# Prediction Intervals — Compatibility Reference

## Interval Scale: Quantiles (Changed in 0.23.0)

Since skforecast 0.23.0, `interval` is expressed as **quantiles in `[0, 1]`**,
not percentiles in `[0, 100]`.

| Input | Behavior |
|-------|----------|
| All values in `[0, 1]`, e.g. `[0.05, 0.95]` | Treated as quantiles (recommended). |
| Single `float`, e.g. `0.9` | Nominal coverage; expands to symmetric quantiles `[0.05, 0.95]`. |
| All values in `(1, 100]`, e.g. `[5, 95]` | Legacy percentiles. Converted to quantiles with a `FutureWarning`. Removed in 0.25.0. |
| Mixed scale, e.g. `[0.05, 95]` | `ValueError` — ambiguous scale. |

`predict_quantiles(quantiles=...)` already used the 0-1 scale and is unaffected.

## Method Compatibility by Forecaster

| Forecaster | `'bootstrapping'` | `'conformal'` | Built-in (parametric) | Default method |
|------------|:------------------:|:-------------:|:---------------------:|:--------------|
| `ForecasterRecursive` | ✓ | ✓ | — | `'bootstrapping'` |
| `ForecasterDirect` | ✓ | ✓ | — | `'bootstrapping'` |
| `ForecasterRecursiveMultiSeries` | ✓ | ✓ | — | `'conformal'` |
| `ForecasterDirectMultiVariate` | ✓ | ✓ | — | `'conformal'` |
| `ForecasterEquivalentDate` | — | ✓ | — | `'conformal'` |
| `ForecasterRnn` | — | ✓ | — | `'conformal'` |
| `ForecasterStats` | — | — | ✓ | N/A (uses `alpha` or `interval`) |
| `ForecasterRecursiveClassifier` | — | — | — | N/A (use `predict_proba()`) |

## Parameter Differences by Forecaster

### ML Forecasters (Recursive, Direct, RecursiveMultiSeries, DirectMultiVariate)

```python
forecaster.predict_interval(
    steps,
    method='bootstrapping',          # or 'conformal'
    interval=[0.05, 0.95],           # quantiles → 90% interval (or a float, e.g. 0.9)
    n_boot=250,                      # only used with method='bootstrapping'
    use_in_sample_residuals=True,
    use_binned_residuals=True,
    random_state=123,
)
```

### ForecasterStats (different interface)

```python
forecaster.predict_interval(
    steps,
    alpha=0.05,                      # significance level → 95% interval
    interval=None,                   # list[float] | None — quantiles, alternative to alpha
)
# NOTE: No `method`, `n_boot`, `use_in_sample_residuals`, `use_binned_residuals`,
#       or `random_state` parameters.
# Use EITHER `alpha` OR `interval`, not both.
```

### ForecasterEquivalentDate (conformal only)

```python
forecaster.predict_interval(
    steps,
    method='conformal',              # only valid value
    interval=[0.05, 0.95],
    use_in_sample_residuals=True,
    use_binned_residuals=True,
    random_state=None,               # Any, accepted but ignored
    exog=None,                       # Any, accepted but ignored
    n_boot=None,                     # Any, accepted but ignored
)
# NOTE: `random_state`, `exog`, and `n_boot` exist for API compatibility but are ignored.
```

### ForecasterRnn (conformal only)

```python
forecaster.predict_interval(
    steps=None,
    levels=None,
    method='conformal',              # only valid value
    interval=[0.05, 0.95],
    use_in_sample_residuals=True,
    use_binned_residuals=True,       # bool, selects residuals by predicted value level
    n_boot=None,                     # Any, accepted but ignored
    random_state=None,               # Any, accepted but ignored
)
# NOTE: `n_boot` and `random_state` exist for API compatibility but are ignored.
```

## Backtesting with Intervals

| Backtesting function | `interval_method` default | Supports `n_boot` | Supports `alpha` |
|---------------------|:------------------------:|:-----------------:|:----------------:|
| `backtesting_forecaster` | `'bootstrapping'` | ✓ | — |
| `backtesting_forecaster_multiseries` | `'conformal'` | ✓ | — |
| `backtesting_stats` | N/A | — | ✓ |

## Prerequisites for Bootstrapping

Bootstrapping requires stored residuals. Two approaches:

### In-sample residuals (simpler)

```python
forecaster.fit(y=y_train, store_in_sample_residuals=True)
forecaster.predict_interval(
    steps=10,
    method='bootstrapping',
    use_in_sample_residuals=True,
)
```

### Out-of-sample residuals (better calibration)

```python
# Step 1: Get predictions via backtesting
metric, preds = backtesting_forecaster(forecaster=forecaster, y=data, cv=cv, metric='mae')

# Step 2: Store residuals
forecaster.set_out_sample_residuals(
    y_true=data.loc[preds.index],
    y_pred=preds['pred'],
)

# Step 3: Use out-of-sample residuals
forecaster.predict_interval(
    steps=10,
    method='bootstrapping',
    use_in_sample_residuals=False,  # Use out-of-sample
)
```

## Binned Residuals

When `use_binned_residuals=True`, residuals are selected based on the predicted
value level (using KBinsDiscretizer). This produces better-calibrated intervals
because residual variance often depends on the prediction magnitude.

| Forecaster | Supports `use_binned_residuals` |
|------------|:-------------------------------:|
| ForecasterRecursive | ✓ |
| ForecasterDirect | ✓ |
| ForecasterRecursiveMultiSeries | ✓ |
| ForecasterDirectMultiVariate | ✓ |
| ForecasterEquivalentDate | ✓ |
| ForecasterRnn | ✓ |
| ForecasterStats | — |

## Probabilistic Prediction Methods Beyond Intervals

| Method | Available in | Description |
|--------|-------------|-------------|
| `predict_interval()` | All except Classifier | Lower/upper bounds at given quantiles |
| `predict_quantiles()` | Recursive, Direct, RecursiveMultiSeries, DirectMultiVariate | Arbitrary quantile predictions |
| `predict_dist()` | Recursive, Direct, RecursiveMultiSeries, DirectMultiVariate | Fit a scipy distribution to bootstrapped predictions |
| `predict_proba()` | RecursiveClassifier only | Class probabilities |

### predict_quantiles signature

```python
forecaster.predict_quantiles(
    steps,
    quantiles=[0.05, 0.5, 0.95],    # any list of quantiles between 0 and 1
    n_boot=250,
    use_in_sample_residuals=True,
    use_binned_residuals=True,
    random_state=123,
)
# Returns DataFrame with one column per quantile: q_0.05, q_0.5, q_0.95
```

### predict_dist signature

```python
from scipy.stats import norm

forecaster.predict_dist(
    steps,
    distribution=norm,               # any scipy.stats distribution
    n_boot=250,
    use_in_sample_residuals=True,
    use_binned_residuals=True,
    random_state=123,
)
# Returns DataFrame with fitted distribution parameters (loc, scale, etc.)
```

## Common Error Patterns

| Error | Cause | Fix |
|-------|-------|-----|
| "No in-sample residuals stored" | `store_in_sample_residuals=False` in `fit()` | Set `store_in_sample_residuals=True` |
| `FutureWarning` about percentiles | `interval` passed as percentiles, e.g. `[5, 95]` | Use quantiles, e.g. `interval=[0.05, 0.95]` (removed in 0.25.0) |
| `ValueError` about ambiguous scale | `interval` mixes quantile and percentile values, e.g. `[0.05, 95]` | Use a single scale, e.g. `interval=[0.05, 0.95]` |
| `method='bootstrapping'` on ForecasterRnn | Bootstrapping not supported | Use `method='conformal'` |
| `method='bootstrapping'` on ForecasterEquivalentDate | Bootstrapping not supported | Use `method='conformal'` |
| Using `alpha` on ML forecasters | Only `ForecasterStats` accepts `alpha` | Use `interval=[lo, hi]` instead |
| Using `method` on ForecasterStats | Stats models use built-in parametric intervals | Remove `method`, use `alpha` or `interval` |

================================================================================
# SKILL: autocorrelation-and-lag-selection
================================================================================

# Autocorrelation and Lag Selection

## When to Use This Skill

Use this skill when the user wants to:

- Inspect the ACF / PACF of a series before training a forecaster.
- Pick a sensible starting set of lags for `ForecasterRecursive`, `ForecasterDirect`, or any multi-series forecaster.
- Detect seasonality (peaks at multiples of the period) from autocorrelation.
- Replace `statsmodels.tsa.stattools.acf` / `pacf` calls with a faster equivalent on long series.

### Related skills

- **After**: `feature-engineering` (turn the candidate lags into a full feature pipeline with rolling and calendar features)
- **After**: `hyperparameter-optimization` (refine the candidate set via cross-validation with `bayesian_search_forecaster`)
- **After**: `feature-selection` (prune redundant lags or exog variables with `select_features` once the forecaster is configured)

## Overview

| Function | Module | Purpose |
|----------|--------|---------|
| `acf` | `skforecast.stats` | Biased ACF via FFT; optional Bartlett confidence intervals |
| `pacf` | `skforecast.stats` | PACF via FFT + Levinson&ndash;Durbin; white-noise CI |
| `calculate_lag_autocorrelation` | `skforecast.stats` | Tabular ACF + PACF, ranked by &#124;PACF&#124; |

> The skforecast functions mirror the public API of `statsmodels.tsa.stattools.acf` / `pacf` but are reimplemented from scratch (FFT for the ACF, Levinson&ndash;Durbin for the PACF) to **maximise speed on long series**. Only the most common configuration options are exposed.
>
> If the user needs the wider catalogue offered by statsmodels — alternative estimators (`'ywunbiased'`, `'ols'`, `'ld'`, `'burg'`), Ljung&ndash;Box statistics, etc. — recommend calling `statsmodels.tsa` directly.
>
> For plotting, keep the statsmodels helpers with their fastest configuration: `plot_acf(..., alpha=0.05, fft=True)` and `plot_pacf(..., alpha=0.05, method='burg')`.

## Computing the ACF

```python
from skforecast.stats import acf

# Just the values (lag 0 included; output length = nlags + 1)
acf_vals = acf(y, nlags=20)

# With 95% Bartlett confidence intervals
acf_vals, confint = acf(y, nlags=20, alpha=0.05)

# Unbiased estimator (denominator n - k)
acf_vals = acf(y, nlags=20, adjusted=True)
```

Notes:
- `nlags` must satisfy `0 < nlags < len(x)`. If `None`, defaults to `min(int(10 * log10(n)), n - 1)`.
- `acf_vals[0]` is always `1.0`.
- Default is the **biased estimator** (`adjusted=False`, denominator `n`). It guarantees a positive semi-definite Toeplitz matrix and is required for downstream PACF computation.

## Computing the PACF

```python
from skforecast.stats import pacf

pacf_vals = pacf(y, nlags=20)
pacf_vals, confint = pacf(y, nlags=20, alpha=0.05)
```

Notes:
- `nlags` must satisfy `0 < nlags < len(x) // 2` — Levinson&ndash;Durbin becomes unreliable when the AR order approaches half the sample size; the function raises explicitly.
- Confidence intervals are asymptotic under the white-noise null: `±z_{α/2} / sqrt(n)` for all lags `≥ 1`. Lag 0 has no uncertainty.
- The implementation uses the biased ACF internally, so values can differ slightly (~2e-2 for typical series) from `statsmodels.pacf(method='yw')`, which uses the unbiased estimator.

## Reading ACF / PACF Patterns

| Pattern | Suggests |
|---------|----------|
| ACF cuts off after lag *q* | MA(q) component |
| ACF decays slowly | AR or mixed process — use PACF |
| PACF cuts off after lag *p* | AR(p) component (lags 1..*p* are the relevant ones) |
| PACF decays slowly | MA or mixed process |
| ACF peaks at lags *k*, *2k*, *3k* | Seasonality of period *k* |
| Persistent slow decay touching every lag | Trend — difference the series before reading orders |

```python
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, axes = plt.subplots(1, 2, figsize=(11, 3))
plot_acf(y, lags=36, alpha=0.05, fft=True, ax=axes[0])
plot_pacf(y, lags=36, alpha=0.05, method='burg', ax=axes[1])
plt.tight_layout()
```

## Ranking Lags with `calculate_lag_autocorrelation`

```python
from skforecast.stats import calculate_lag_autocorrelation

# Default: rows ordered by |PACF| descending, lags 1..n_lags
lag_table = calculate_lag_autocorrelation(
    data=y,
    n_lags=24,
    sort_by='partial_autocorrelation_abs',
)

# Other valid sort keys
calculate_lag_autocorrelation(data=y, n_lags=24, sort_by='lag')
calculate_lag_autocorrelation(data=y, n_lags=24, sort_by='autocorrelation_abs')
```

Returns a DataFrame with columns
`['lag', 'partial_autocorrelation_abs', 'partial_autocorrelation', 'autocorrelation_abs', 'autocorrelation']`
covering lags 1..`n_lags` (lag 0 dropped because it is always 1.0).

`last_n_samples` keeps only the most recent observations — useful when the series is very long and only the recent dynamics matter:

```python
calculate_lag_autocorrelation(data=y, n_lags=24, last_n_samples=2000)
```

`data` accepts a `pd.Series` or a single-column `pd.DataFrame`. `n_lags` must be strictly less than `len(data) // 2`.

## End-to-End: From PACF to a Forecaster

```python
import numpy as np
from skforecast.stats import calculate_lag_autocorrelation
from skforecast.recursive import ForecasterRecursive
from sklearn.ensemble import HistGradientBoostingRegressor

# 1. Rank lags by |PACF|
lag_table = calculate_lag_autocorrelation(data=y_train, n_lags=24)

# 2. Threshold against the asymptotic 95% white-noise band (1.96 / sqrt(n))
threshold = 1.96 / np.sqrt(len(y_train))
selected_lags = (
    lag_table
    .loc[lag_table['partial_autocorrelation_abs'] > threshold, 'lag']
    .astype(int)
    .sort_values()
    .tolist()
)

# 3. Feed directly into any forecaster
forecaster = ForecasterRecursive(
    estimator=HistGradientBoostingRegressor(random_state=963),
    lags=selected_lags,
)
forecaster.fit(y=y_train)
predictions = forecaster.predict(steps=12)
```

The selected lags are an *informed starting point*, not a final answer — refine with `bayesian_search_forecaster` (see `hyperparameter-optimization`) or prune with `select_features` (see `feature-selection`).

## Common Mistakes

1. **Using `nlags >= n // 2` for PACF.** Levinson&ndash;Durbin loses stability as the AR order approaches half the sample size; the function raises explicitly. Reduce `nlags` or use a longer series.
2. **Reading the ACF to choose AR order.** ACF is contaminated by indirect dependencies through earlier lags and rarely cuts off cleanly for AR processes — use the PACF for AR order, the ACF for MA order.
3. **Treating the confidence bands as joint tests.** Bartlett (ACF) and white-noise (PACF) bands are *pointwise*. Several lags will cross the band by chance; use them to spot candidates, not to formally reject white noise.
4. **Skipping detrending / differencing.** A strong trend dominates the ACF and masks higher-lag structure. If the series shows visible trend, read ACF / PACF on `y.diff().dropna()` (or use `differentiation=1` in the forecaster).
5. **Switching to `adjusted=True` by default.** The biased estimator is preferred for forecasting workflows; the unbiased one can produce `|ACF| > 1` at high lags and breaks the positive semi-definite property required by Levinson&ndash;Durbin.
6. **Expecting statsmodels-only options.** `method='ols'`, `method='burg'`, alternative confidence intervals, etc. are *not* exposed by `skforecast.stats.pacf`. For those, call `statsmodels.tsa.stattools.pacf` directly — the API is the same.
7. **Passing a multi-column DataFrame to `calculate_lag_autocorrelation`.** It accepts a `Series` or a single-column `DataFrame` and raises otherwise. Select the target column first.

================================================================================
# SKILL: feature-engineering
================================================================================

# Feature Engineering

## References

See [references/calendar-features-reference.md](references/calendar-features-reference.md)
for the complete `CalendarFeatures` constructor, all supported features and
encodings, the delegated (`calendar_features` parameter) vs. manual (`exog`)
workflows, per-forecaster support, and gotchas.

See [references/rolling-stats-reference.md](references/rolling-stats-reference.md) for
the complete `RollingFeatures` constructor, all 9 available statistics,
feature name generation formula, window behavior, and `kwargs_stats` usage.

## When to Use This Skill

Use this skill when the user wants to create features to improve forecasting
accuracy: calendar/datetime features, cyclical / onehot / spline encoding,
holiday distance features, rolling statistics, differencing, or data scaling.

### Related skills

- **Before**: `autocorrelation-and-lag-selection` (use ACF/PACF to choose a candidate set of lags before adding rolling and calendar features)
- **After**: `feature-selection` (prune redundant features with `select_features` after engineering)
- **After**: `hyperparameter-optimization` (jointly tune the engineered configuration and the estimator hyperparameters)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| Do not mix the delegated (`calendar_features=`) and manual (`CalendarFeatures` as exog) paths for the same features | `ValueError` on duplicate feature names | Choose one path per feature set |
| Manual calendar / holiday exog must cover the full forecast horizon | `predict()` exog error or truncated forecast | Generate exog for all future dates before predicting |
| Cyclical calendar features need sin/cos (or spline) encoding | Hour 23 treated as far from hour 0, degrading the model | Use `encoding='cyclical'` (or spline) for periodic features |
| `window_features` are computed on the differenced series when `differentiation` is set | `roll_mean_7` is a mean of changes, not of raw values | Account for differencing when interpreting or scaling window features |

## Overview

| Tool | Module | Purpose |
|------|--------|---------|
| `calendar_features` param | skforecast ML forecasters | **Delegated** calendar features: pass a `CalendarFeatures` instance and the forecaster generates them at train and predict — no manual exog. **New in 0.23.0.** |
| `CalendarFeatures` | `skforecast.preprocessing` | Sklearn-compatible transformer: extract calendar features from a `DatetimeIndex` and (optionally) encode them. Pass to `calendar_features`, use as `transformer_exog`, or in a `Pipeline`. |
| `create_calendar_features` | `skforecast.preprocessing` | Function form of the same logic, for one-shot use without a transformer. |
| `calculate_distance_from_holiday` | `skforecast.preprocessing` | Periods to next / since last holiday |
| `RollingFeatures` | `skforecast.preprocessing` | Rolling window statistics (mean, std, min, max, etc.) |
| `differentiation` param | skforecast forecasters | Make non-stationary series stationary |

## Calendar Features

There are **two ways** to add calendar features (month, day of week, hour, …).
Both use the `CalendarFeatures` class; they differ in who builds the features.

| Workflow | How | When to use |
|----------|-----|-------------|
| **Delegated** (new in 0.23.0) | Pass a `CalendarFeatures` instance to the forecaster's `calendar_features` parameter. The forecaster generates the features at train **and** predict — no manual exog. | The 4 supported forecasters: `ForecasterRecursive`, `ForecasterRecursiveMultiSeries`, `ForecasterDirect`, `ForecasterDirectMultiVariate`. |
| **Manual** | Build features with `CalendarFeatures.fit_transform` / `create_calendar_features` and pass them as `exog` (or wire `CalendarFeatures` as `transformer_exog`). | Forecasters without `calendar_features` support (`ForecasterRecursiveClassifier`, `ForecasterRnn`, `ForecasterStats`, `ForecasterFoundation`, `ForecasterEquivalentDate`), or when the features live inside a `Pipeline` / `ColumnTransformer`. |

> Full constructor, all features/encodings, per-forecaster support, and gotchas:
> [references/calendar-features-reference.md](references/calendar-features-reference.md).

**Supported features:** `'year'`, `'month'`, `'week'`, `'day_of_week'`,
`'day_of_month'`, `'day_of_year'`, `'weekend'`, `'hour'`, `'minute'`,
`'second'`, `'quarter'`. **Encodings:** `'cyclical'` (default, sin/cos pair),
`'onehot'`, `'spline'`, `None`. `'year'` and `'weekend'` are never encoded.

### Delegated — `calendar_features` parameter (preferred when supported)

The forecaster builds the calendar features from the datetime index during both
training and prediction. You do **not** build or pass a calendar `exog`, and you
do **not** need to cover the forecast horizon manually.

```python
from lightgbm import LGBMRegressor
from skforecast.preprocessing import CalendarFeatures
from skforecast.recursive import ForecasterRecursive

calendar = CalendarFeatures(
    features=['month', 'day_of_week', 'hour'],
    encoding='cyclical',                 # 'cyclical' | 'onehot' | 'spline' | None
    keep_original_columns=False,
)

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    calendar_features=calendar,          # requires a pandas DatetimeIndex
)
forecaster.fit(y=y_train)                # calendar features built automatically
predictions = forecaster.predict(steps=24)   # calendar features built automatically
```

Requires a `DatetimeIndex` (otherwise `TypeError`). Calendar features are added
as predictors alongside lags, window features, and any `exog`.

### Manual — build calendar features as exog

For forecasters without the `calendar_features` parameter, or for use in a
`Pipeline`:

```python
import pandas as pd
from skforecast.preprocessing import CalendarFeatures

data = data.asfreq('h')                  # DatetimeIndex with frequency set

calendar_transformer = CalendarFeatures(
    features=['month', 'week', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)
exog_calendar = calendar_transformer.fit_transform(data)
# Columns: month_sin, month_cos, week_sin, week_cos, ...

# The exog passed to predict() must cover the entire forecast horizon.
```

The function form `create_calendar_features(X=data, ...)` does the same in one
shot without a transformer. See the reference for `max_values`, `spline_kwargs`,
and `features_to_encode` options.

## Holiday Distance Features

`calculate_distance_from_holiday` computes the number of periods to the next
and since the last holiday. The time unit is inferred from the index frequency
(days, hours, minutes, …) when `date_column=None`, and is always days when a
date column is used.

```python
from skforecast.preprocessing import calculate_distance_from_holiday

# Index-based (preferred): unit follows the index frequency
holiday_dist = calculate_distance_from_holiday(
    X=data[['is_holiday']],          # boolean / 0-1 column
    holiday_column='is_holiday',
    fill_na=0,                       # value for rows before first / after last holiday
)
# Columns: time_to_holiday, time_since_holiday

# Series form: values used directly as the holiday indicator
holiday_dist = calculate_distance_from_holiday(
    X=data['is_holiday'],
    fill_na=0,
)
```

Combine with calendar features into a single `exog`:

```python
exog = pd.concat([exog_calendar, holiday_dist], axis=1)
```

## Rolling Features (Window Statistics)

```python
from skforecast.preprocessing import RollingFeatures
from skforecast.recursive import ForecasterRecursive
from lightgbm import LGBMRegressor

# Single window size for all stats
rolling = RollingFeatures(
    stats=['mean', 'std', 'min', 'max'],
    window_sizes=7,                  # int applies same window to all stats
)

# Different window sizes per statistic
rolling = RollingFeatures(
    stats=['mean', 'std', 'min', 'max'],
    window_sizes=[7, 7, 14, 14],     # Must match length of stats
)

# Multiple RollingFeatures objects
rolling_short = RollingFeatures(stats=['mean', 'std'], window_sizes=7)
rolling_long = RollingFeatures(stats=['mean', 'std'], window_sizes=30)

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    window_features=[rolling_short, rolling_long],   # List of RollingFeatures
)
```

### Available Rolling Statistics

Standard: `'mean'`, `'std'`, `'min'`, `'max'`, `'sum'`, `'median'`, `'ratio_min_max'`, `'coef_variation'`

Exponential weighted: `'ewm'` — requires `kwargs_stats`:
```python
rolling = RollingFeatures(
    stats=['ewm'],
    window_sizes=7,
    kwargs_stats={'ewm': {'alpha': 0.3}},
)
```

### The `window_features` parameter

The forecaster's `window_features` argument accepts **any object (or list of objects)** that implements the contract:

```python
def transform_batch(y: pd.Series) -> pd.DataFrame:
    # Returns one row per training sample, indexed to align with y[max_window_size:]
    ...
```

`RollingFeatures` is the built-in implementation; user-defined classes that follow the same contract are accepted. See [references/rolling-stats-reference.md](references/rolling-stats-reference.md#custom-window-feature-classes-protocol) for the full protocol and a custom-class example.

**Effective `window_size`.** The forecaster sets `window_size = max(max_lag, max_size_window_features) (+ differentiation)`. A `RollingFeatures(window_sizes=30)` combined with `lags=7` raises `window_size` to 30 — `last_window` must hold ≥ 30 values and the first 30 rows of training data are dropped. This is a frequent source of "not enough data to train" errors when the rolling window is much larger than the lags.

**Lags-free configuration.** `lags=None, window_features=...` is valid (and useful for purely smoothed-history models). At least one of `lags` / `window_features` must be non-`None` — passing both as `None` raises `ValueError`.

**Multiple `RollingFeatures` instances.** When passing a list, each instance is processed independently and their feature names are concatenated. Names must be unique across instances; collisions silently overwrite columns. Use the `features_names=` argument of `RollingFeatures` to disambiguate when reusing the same `stat`/`window_size` combination.

**Runtime mutation.** After construction, use `forecaster.set_window_features(new_window_features)` to swap the configuration (analogous to `set_lags`). This recomputes `window_size`, `max_size_window_features`, `window_features_names`, and the differentiator's window. Useful inside manual hyperparameter loops.

**Interaction with `differentiation`.** Window features are computed on the **differenced** series, not the original one. With `differentiation=1`, `roll_mean_7` is the mean of seven consecutive *changes*, not seven raw values. Take this into account when interpreting feature importances.

**Multi-series.** In `ForecasterRecursiveMultiSeries`, the same `RollingFeatures` instance is applied independently to each series — configuration is global but the computed values are per-series (no leakage between series).

**Classifier variant.** For `ForecasterRecursiveClassifier`, use `RollingFeaturesClassification` (also in `skforecast.preprocessing`) instead of `RollingFeatures`.

### Choosing window features

**When window features help:**
- Noisy series — a smoothed signal (`mean`, `ewm`) reduces variance the lags cannot.
- Capturing a longer scale than the lags without exploding column count: with hourly data and `lags=24`, a `roll_mean_168` represents the weekly level with one feature instead of adding 144 more lags.
- Heteroscedastic or regime-changing series (energy, finance, traffic) — `roll_std`, `coef_variation`.
- When extreme behaviour drives decisions (peak demand, drawdowns) — `roll_min`, `roll_max`.

**When they don't help:**
- Short series — each large window discards training rows (`window_size` grows).
- Highly periodic, low-noise series — lags + calendar features are usually enough.
- Redundant ranges — a `roll_mean_24` alongside `lags=range(1, 25)` adds little for a gradient boosting model that can already construct the average from splits.

**Which statistic for which signal:**

| Stat | Adds signal when… | Notes |
|------|------------------|-------|
| `mean` | Default smoothing — recent level. | Start here. |
| `ewm` | Recent observations should dominate. | Tune `alpha` in `kwargs_stats` (0.1 slow, 0.5 reactive). |
| `std` | Volatility / heteroscedasticity matters. | Near-mandatory in finance and energy demand. |
| `median` | Frequent outliers. | Replaces `mean`, do not combine. |
| `min` / `max` | Threshold crossings, peaks, operational floors. | Useful in capacity planning. |
| `sum` | Cumulative quantities (rainfall, sales) rather than levels. | Scale-sensitive — pair with `transformer_y` or scaling. |
| `ratio_min_max` | Bounded regime indicator. | Often noisy on stationary series. |
| `coef_variation` | Scale-free volatility. | Useful in `ForecasterRecursiveMultiSeries` to compare series of different magnitudes. |

**Window sizes:**
- Tie them to actual seasonality: hourly → 24, 168; daily → 7, 30; weekly → 4, 13, 52; monthly → 12.
- Multi-scale usually wins: combine one short (reactive) and one long (trend) window — e.g. `stats=['mean', 'std', 'mean', 'std'], window_sizes=[7, 7, 30, 30]`.
- Keep `max_window_size < ~25-30% of len(y_train)` to preserve enough training rows.

**Practical recipe:**
1. **Baseline**: lags only (from ACF/PACF) + calendar features if applicable. Measure via backtesting.
2. **If residuals track recent level** → add `roll_mean_<seasonal_period>`.
3. **If errors scale with magnitude** → add `roll_std_<seasonal_period>`.
4. **If outliers/peaks matter** → add `roll_max` or `roll_min` with the same window.
5. **One addition per iteration**, validated with backtesting. Drop it if MAE/RMSE doesn't improve on independent folds.
6. **Prune with `select_features`** at the end — typically 2-4 of 6-8 candidates survive.

**Combinations that tend to work:**
- **Hourly electricity demand**: `lags=[1..24, 168]` + `RollingFeatures(['mean', 'std'], [24, 24])` + cyclical calendar.
- **Daily sales with promotions**: `lags=[1..14]` + `RollingFeatures(['mean', 'median'], [7, 28])` + holiday distance + promo exog.
- **Financial series**: `lags=[1..5]` + `RollingFeatures(['std', 'coef_variation'], [5, 20])` + `differentiation=1` (window features compute on returns, which is usually what you want).
- **Multi-series with heterogeneous scales**: prefer `coef_variation` and `ratio_min_max` (scale-free), or apply per-series `transformer_series`.

**Anti-patterns:**
- Adding 6-8 statistics at once "just in case" — dilutes signal and increases overfitting risk.
- `roll_sum` without scaling the target — its magnitude dominates the lags and destabilises training.
- `min` / `max` with very long windows on trending series — the feature becomes near-constant or monotone.
- Duplicating information: keeping `lags=range(1, 31)` together with `roll_mean_30` rarely helps a gradient boosting model — shrink the lag set when you add the rolling stat.

## Differencing (Non-Stationary Series)

```python
# Built-in — forecaster handles differencing and inverse transform automatically
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    differentiation=1,               # First-order differencing (removes linear trend)
    # differentiation=2,             # Second-order (removes quadratic trend)
)
forecaster.fit(y=y_train)
predictions = forecaster.predict(steps=10)   # Auto inverse-transformed
```

## Categorical Exogenous Variables

All ML forecasters include a `categorical_features` parameter (default `'auto'`)
that automatically detects and encodes non-numeric exogenous columns using an
internal `OrdinalEncoder` (into float codes, since numpy arrays are used internally).
Native categorical support is configured automatically for LightGBM, CatBoost,
XGBoost, and HistGradientBoostingRegressor.

```python
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    categorical_features='auto',     # Default — auto-detect non-numeric columns
)
```

**`categorical_features` options:**
- `'auto'` (default): Auto-detect non-numeric columns after `transformer_exog`.
- `list`: Explicit column names to treat as categorical (including numeric columns).
- `None`: No internal categorical encoding.

**Important:** When `categorical_features` is not `None`, do not set categorical
features directly on the estimator or via `fit_kwargs`. The forecaster manages
the configuration internally and overwrites estimator-level settings.

**Choosing an encoding strategy:**

| Method | API | Best for |
|--------|-----|----------|
| Built-in `categorical_features` | `categorical_features='auto'` or `list` | Gradient boosting (LightGBM, XGBoost, CatBoost, HistGBR) — simplest workflow |
| One-hot / Ordinal encoding | `transformer_exog` | Linear models, SVMs, non-gradient-boosting trees |
| Target encoding | Outside forecaster | High-cardinality features (applied manually to avoid leakage) |

**Combining `transformer_exog` and `categorical_features`:**
`transformer_exog` is applied **before** `categorical_features` detection. Scale
numeric columns with `transformer_exog` while `categorical_features='auto'` handles
the rest. Avoid applying both mechanisms to the same columns.

```python
from sklearn.compose import make_column_transformer
from sklearn.preprocessing import StandardScaler

transformer_exog = make_column_transformer(
    (StandardScaler(), ['temp', 'hum']),
    remainder='passthrough',
    verbose_feature_names_out=False
).set_output(transform='pandas')

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    transformer_exog=transformer_exog,
    categorical_features='auto',     # Detects remaining non-numeric columns
)
```

## Combining Features — Full Example

```python
import pandas as pd
from skforecast.preprocessing import (
    CalendarFeatures,
    calculate_distance_from_holiday,
    RollingFeatures,
)
from skforecast.recursive import ForecasterRecursive
from sklearn.preprocessing import StandardScaler
from lightgbm import LGBMRegressor

# 1. Calendar features with cyclical encoding (no extra deps)
calendar_transformer = CalendarFeatures(
    features=['month', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)
exog_calendar = calendar_transformer.fit_transform(data)

# 2. Holiday distance features (unit inferred from index frequency)
exog_holiday = calculate_distance_from_holiday(
    X=data[['is_holiday']],
    holiday_column='is_holiday',
    fill_na=0,
)

# 3. Combine with other exogenous variables
exog = pd.concat([exog_external, exog_calendar, exog_holiday], axis=1)

# 4. Rolling features + lags + differencing
rolling = RollingFeatures(stats=['mean', 'std'], window_sizes=[7, 14])

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=[1, 2, 3, 7, 14, 24],
    window_features=rolling,
    transformer_y=StandardScaler(),
    differentiation=1,
)
forecaster.fit(y=y_train, exog=exog.loc[y_train.index])
predictions = forecaster.predict(steps=10, exog=exog.loc[forecast_index])
```

### Variant — delegate calendar features to the forecaster

On a forecaster that supports `calendar_features`, drop the manual calendar exog
and let the forecaster build it at train and predict time. Only the
non-calendar exog (holiday + external) is passed manually:

```python
calendar = CalendarFeatures(
    features=['month', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)

exog = pd.concat([exog_external, exog_holiday], axis=1)   # no calendar columns here

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=[1, 2, 3, 7, 14, 24],
    window_features=rolling,
    calendar_features=calendar,          # built automatically at fit/predict
    transformer_y=StandardScaler(),
    differentiation=1,
)
forecaster.fit(y=y_train, exog=exog.loc[y_train.index])
predictions = forecaster.predict(steps=10, exog=exog.loc[forecast_index])
```

## Common Mistakes

1. **Not encoding cyclical features**: Using raw integers for hour/month/day_of_week loses the cyclical relationship (hour 23 appears far from hour 0). Use `encoding='cyclical'` (or `'onehot'` / `'spline'`).
2. **Forgetting frequency on index**: `CalendarFeatures` (and the `calendar_features` parameter) requires a `DatetimeIndex`; `calculate_distance_from_holiday` (index mode) infers the time unit from the frequency. Always set `data.asfreq('h')` (or similar) first.
3. **Not covering forecast horizon with calendar/holiday exog (manual path only)**: When building calendar features manually as `exog`, the `predict()` exog must include future dates covering the entire forecast horizon. With the delegated `calendar_features` parameter this is handled automatically.
4. **Mixing the delegated and manual calendar paths for the same features**: Passing `calendar_features=` *and* also adding manually-built calendar columns with the same names to `exog` raises `ValueError: Duplicated feature names detected in X_train`. Pick one path per feature.
5. **Using `calendar_features` on an unsupported forecaster**: Only `ForecasterRecursive`, `ForecasterRecursiveMultiSeries`, `ForecasterDirect`, and `ForecasterDirectMultiVariate` accept it. For `ForecasterRecursiveClassifier`, `ForecasterRnn`, `ForecasterStats`, `ForecasterFoundation`, and `ForecasterEquivalentDate`, build calendar features manually as `exog`.
6. **Overriding `max_values` for `'week'` or `'day_of_year'`**: The defaults (53, 366) are intentional — they handle ISO week 53 and leap years correctly. Use the smaller value (52 / 365) only if you have verified your data never reaches the maximum.
7. **Over-engineering features**: Start with lags only, then add rolling features and calendar features incrementally. Validate each addition with backtesting.
8. **Rolling window larger than lags without checking training size**: The forecaster's `window_size` is `max(max_lag, max_size_window_features) (+ differentiation)`. A 30-step rolling window with `lags=7` drops the first 30 rows of training data and requires `last_window` of length ≥ 30 at predict time.
9. **Colliding feature names across multiple `RollingFeatures`**: When passing `window_features=[wf1, wf2]`, names like `roll_mean_7` produced by both instances overwrite each other. Override with `features_names=[...]` on at least one instance.
10. **Forgetting that window features run on the differenced series**: With `differentiation=1`, `roll_mean_7` is the mean of seven changes, not seven raw values. Adjust interpretation accordingly.

---
### Reference: calendar-features-reference

# Calendar Features — Reference

Two ways to add calendar/datetime features (month, day of week, hour, …) to a
forecaster:

- **Delegated** — pass a `CalendarFeatures` instance to the forecaster's
  `calendar_features` parameter. The forecaster generates the features
  automatically during **training and prediction**. No manual `exog`.
  **New in skforecast 0.23.0.**
- **Manual** — build the features yourself with `CalendarFeatures`
  (`fit_transform`) or `create_calendar_features` and pass the result as `exog`
  (or wire `CalendarFeatures` as `transformer_exog`). Required for forecasters
  that do not support the `calendar_features` parameter.

Both paths use the same `CalendarFeatures` class, so the constructor and
encoding options below apply to either workflow.

## Which workflow to use

| You want… | Use |
|-----------|-----|
| Simplest workflow, no manual exog at predict time | Delegated (`calendar_features=`) |
| One of the 4 supported forecasters (see table below) | Delegated |
| A forecaster without `calendar_features` support | Manual (`exog` / `transformer_exog`) |
| Calendar features inside a `Pipeline` / `ColumnTransformer` | Manual (`CalendarFeatures` as a transformer) |
| One-shot feature creation outside any forecaster | Manual (`create_calendar_features`) |

### Forecaster support for the `calendar_features` parameter

| Forecaster | `calendar_features` param | Calendar features path |
|------------|:-------------------------:|------------------------|
| `ForecasterRecursive` | ✓ | Delegated or manual |
| `ForecasterRecursiveMultiSeries` | ✓ | Delegated or manual |
| `ForecasterDirect` | ✓ | Delegated or manual |
| `ForecasterDirectMultiVariate` | ✓ | Delegated or manual |
| `ForecasterRecursiveClassifier` | ✗ | Manual only |
| `ForecasterRnn` | ✗ | Manual only |
| `ForecasterStats` | ✗ | Manual only |
| `ForecasterFoundation` | ✗ | Manual only |
| `ForecasterEquivalentDate` | ✗ | Manual only |

> All calendar tools are built into skforecast — no `feature_engine` (or any
> other extra) dependency is required.

## CalendarFeatures Constructor

```python
from skforecast.preprocessing import CalendarFeatures

calendar = CalendarFeatures(
    features=None,                   # list[str] | None — None extracts all supported features
    features_to_encode=None,         # list[str] | None — None encodes all encodable features
    encoding="cyclical",             # 'cyclical' | 'onehot' | 'spline' | None
    max_values=None,                 # dict[str, int] | None — override per-feature period
    spline_kwargs=None,              # dict | None — args for SplineTransformer (encoding='spline')
    keep_original_columns=True,      # bool — merge with the input X's columns
)
```

> The forecaster stores a **clone** of the instance you pass to
> `calendar_features`, so the same object can be reused safely.

### Supported features

`'year'`, `'month'`, `'week'`, `'day_of_week'`, `'day_of_month'`,
`'day_of_year'`, `'weekend'`, `'hour'`, `'minute'`, `'second'`, `'quarter'`.

By default all are extracted; pass `features=[...]` to subset.

### Supported encodings

| `encoding` | Output | Notes |
|-----------|--------|-------|
| `'cyclical'` (default) | `{feature}_sin`, `{feature}_cos` | sin/cos pair per cyclical feature |
| `'onehot'` | One column per known category (e.g. `month_1` … `month_12`) | Stable schema across train / predict |
| `'spline'` | `≈ max_val` columns per feature, periodic B-splines | Smooth alternative to onehot |
| `None` | Raw integer columns | No transformation |

`'year'` and `'weekend'` are **never** encoded (they are not cyclical) and are
always kept as raw integers regardless of `encoding`.

### Choosing an encoding

| Encoding | Columns per feature | Best for |
|----------|---------------------|----------|
| `'cyclical'` | 2 (sin + cos) | Compact, smooth — good default for tree and linear models |
| `'onehot'` | `max_val` | Stable schema; tree models can split on individual categories |
| `'spline'` | `≈ max_val` (dense) | Smooth + flexible; high-cardinality features (`day_of_year`) when memory allows |
| `None` | 1 | Tree models that benefit from raw ordinal values |

For high-cardinality features (`day_of_year` → 366, `day_of_month` → 31),
`'cyclical'` or `'spline'` are typically more memory-efficient than `'onehot'`.

### `max_values` defaults

```python
{'month': 12, 'week': 53, 'day_of_week': 7, 'day_of_month': 31,
 'day_of_year': 366, 'hour': 24, 'minute': 60, 'second': 60, 'quarter': 4}
```

These handle leap years and ISO week 53 correctly. Override only the keys you
need:

```python
calendar = CalendarFeatures(
    features=['month', 'hour'],
    encoding='cyclical',
    max_values={'month': 6},   # Custom semester period; hour keeps default 24
)
```

### `spline_kwargs`

`spline_kwargs` accepts any argument of
`sklearn.preprocessing.SplineTransformer` **except** `knots` (computed
internally from `max_values`) and `sparse_output` (incompatible with the
DataFrame output).

```python
calendar = CalendarFeatures(
    features=['day_of_year'],
    encoding='spline',
    spline_kwargs={'degree': 3, 'n_knots': 12},
)
```

### `features_to_encode` — encode only some features

Extract a feature but leave it as a raw integer:

```python
calendar = CalendarFeatures(
    features=['year', 'month', 'hour'],
    features_to_encode=['month', 'hour'],   # 'year' kept as raw int
    encoding='cyclical',
)
```

## Workflow A — Delegated (`calendar_features` parameter)

The forecaster generates calendar features from the datetime index during both
training and prediction — you never build or pass a calendar `exog`.

```python
from lightgbm import LGBMRegressor
from skforecast.preprocessing import CalendarFeatures
from skforecast.recursive import ForecasterRecursive

calendar = CalendarFeatures(
    features=['month', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    calendar_features=calendar,
)

# Index must be a pandas DatetimeIndex with frequency set
forecaster.fit(y=y_train)                       # calendar features built automatically
predictions = forecaster.predict(steps=24)      # calendar features built automatically
```

Key points:

- **Requires a `DatetimeIndex`.** A non-datetime index raises `TypeError`.
- **No horizon coverage needed.** The forecaster derives the prediction index
  itself, so unlike manual `exog` you do not pass future calendar values to
  `predict`.
- Calendar features are added as predictors **alongside** lags, window
  features, and any `exog`.
- Works in multi-series forecasters too: the same configuration is applied to
  the (shared) datetime index, producing identical calendar values per series.

### Related attributes (after `fit`)

| Attribute | Meaning |
|-----------|---------|
| `forecaster.calendar_features` | The cloned `CalendarFeatures` instance (or `None`) |
| `forecaster.calendar_features_names` | The `features` requested (e.g. `['month', 'hour']`) |
| `forecaster.X_train_calendar_features_names_out_` | Output column names added to `X_train` |

## Workflow B — Manual (`exog` / `transformer_exog` / function)

Use this for forecasters without `calendar_features` support, or when you need
the features inside a `Pipeline` / `ColumnTransformer`.

### As an exog DataFrame

```python
calendar = CalendarFeatures(
    features=['month', 'week', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)
exog_calendar = calendar.fit_transform(data)
# Columns: month_sin, month_cos, week_sin, week_cos,
#          day_of_week_sin, day_of_week_cos, hour_sin, hour_cos

calendar.get_feature_names_out()   # resolved output column names (sklearn API)
```

The `exog` passed to `predict()` must include future dates covering the **entire
forecast horizon**.

### Function form — `create_calendar_features`

For one-shot use without instantiating a transformer:

```python
from skforecast.preprocessing import create_calendar_features

exog_calendar = create_calendar_features(
    X=data,
    features=['month', 'day_of_week', 'hour'],
    encoding='cyclical',
    keep_original_columns=False,
)
```

All parameters match `CalendarFeatures`. Prefer the transformer when you want to
fit / re-apply the same configuration, plug it into a `Pipeline`, or pass it as
`transformer_exog`.

## Gotchas

- **Datetime index required (delegated path).** A non-`DatetimeIndex` raises
  `TypeError`. Set the frequency first (`data.asfreq('h')`).
- **Do not mix paths for the same features.** Passing `calendar_features=` and
  also adding manually-built calendar columns to `exog` with the same names
  raises `ValueError: Duplicated feature names detected in X_train`. Pick one
  path per feature.
- **`keep_original_columns` default is `True`.** When used to build `exog`
  manually, set `keep_original_columns=False` if you only want the calendar
  columns (the constructor default keeps the input `X` columns too).
- **`max_values` for `'week'` / `'day_of_year'`.** Defaults (53, 366) handle ISO
  week 53 and leap years. Use 52 / 365 only after verifying your data never
  reaches the maximum.

---
### Reference: rolling-stats-reference

# Rolling Features — Statistics Reference

## RollingFeatures Constructor

```python
from skforecast.preprocessing import RollingFeatures

rolling = RollingFeatures(
    stats,                           # str | list[str] (required)
    window_sizes,                    # int | list[int] (required)
    min_periods=None,                # int | list[int] | None — defaults to window_sizes
    features_names=None,             # list[str] | None — auto-generated if None
    fillna=None,                     # str | float | None — fill NaN in transform_batch
    kwargs_stats=                    # dict | None — extra args per statistic
        {'ewm': {'alpha': 0.3}},     # default
)
```

## Available Statistics

| Statistic | String | Description | Requires `kwargs_stats` |
|-----------|--------|-------------|:-:|
| Mean | `'mean'` | Arithmetic mean of the window | — |
| Std deviation | `'std'` | Standard deviation of the window | — |
| Minimum | `'min'` | Minimum value in the window | — |
| Maximum | `'max'` | Maximum value in the window | — |
| Sum | `'sum'` | Sum of values in the window | — |
| Median | `'median'` | Median value in the window | — |
| Min/Max ratio | `'ratio_min_max'` | `min / max` of the window | — |
| Coef. variation | `'coef_variation'` | `std / mean` of the window | — |
| Exp. weighted mean | `'ewm'` | Exponentially weighted mean | ✓ `{'ewm': {'alpha': ...}}` |

## Feature Name Generation

Default names follow the pattern: `roll_{stat}_{window_size}`

| Stats | Window | Generated name |
|-------|--------|---------------|
| `'mean'` | `7` | `roll_mean_7` |
| `'std'` | `14` | `roll_std_14` |
| `'ewm'` | `7` (alpha=0.3) | `roll_ewm_7_alpha_0.3` |

Override with custom names:

```python
rolling = RollingFeatures(
    stats=['mean', 'std'],
    window_sizes=[7, 14],
    features_names=['weekly_avg', 'biweekly_std'],
)
```

## Window Behavior

Rolling windows use `closed='left'` and `center=False` to avoid data leakage.
The last point in the window is **excluded** from calculations:

```
Window size = 3, calculating for time t:
Uses values: [t-3, t-2, t-1]  (NOT t itself)
```

## Configuration Patterns

### Same window for all stats

```python
rolling = RollingFeatures(
    stats=['mean', 'std', 'min', 'max'],
    window_sizes=7,                   # int → applied to all 4 stats
)
# Features: roll_mean_7, roll_std_7, roll_min_7, roll_max_7
```

### Different windows per stat

```python
rolling = RollingFeatures(
    stats=['mean', 'std', 'min', 'max'],
    window_sizes=[7, 7, 14, 14],      # list → must match length of stats
)
# Features: roll_mean_7, roll_std_7, roll_min_14, roll_max_14
```

### Repeated stats with different windows

```python
rolling = RollingFeatures(
    stats=['mean', 'mean', 'std', 'std'],
    window_sizes=[7, 30, 7, 30],
)
# Features: roll_mean_7, roll_mean_30, roll_std_7, roll_std_30
```

### Multiple RollingFeatures objects

```python
rolling_short = RollingFeatures(stats=['mean', 'std'], window_sizes=7)
rolling_long = RollingFeatures(stats=['mean', 'std'], window_sizes=30)

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    window_features=[rolling_short, rolling_long],  # list of RollingFeatures
)
```

### Exponentially weighted mean

```python
rolling = RollingFeatures(
    stats=['ewm'],
    window_sizes=7,
    kwargs_stats={'ewm': {'alpha': 0.3}},
)
# Feature: roll_ewm_7_alpha_0.3

# Multiple ewm with different alphas (use separate objects)
ewm_fast = RollingFeatures(stats=['ewm'], window_sizes=7, kwargs_stats={'ewm': {'alpha': 0.5}})
ewm_slow = RollingFeatures(stats=['ewm'], window_sizes=7, kwargs_stats={'ewm': {'alpha': 0.1}})

forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    window_features=[ewm_fast, ewm_slow],
)
```

## min_periods Parameter

Controls the minimum number of observations required to compute a statistic.

```python
# Default: min_periods = window_sizes (requires full window)
rolling = RollingFeatures(stats=['mean'], window_sizes=7)
# Equivalent to: min_periods=7 → first 6 values are NaN

# Allow partial windows
rolling = RollingFeatures(stats=['mean'], window_sizes=7, min_periods=1)
# min_periods=1 → computes mean even with 1 observation
```

## fillna Parameter

Used only in `transform_batch()` method (not in recursive prediction):

| Value | Behavior |
|-------|----------|
| `None` | No filling (NaN remains) |
| `'mean'` | Fill with mean of the feature |
| `'median'` | Fill with median of the feature |
| `'ffill'` | Forward fill |
| `'bfill'` | Backward fill |
| `float` | Fill with specific value (e.g., `0.0`) |

## Forecaster Compatibility

| Forecaster | `window_features` supported | Notes |
|------------|:--:|-------|
| ForecasterRecursive | ✓ | |
| ForecasterDirect | ✓ | |
| ForecasterRecursiveMultiSeries | ✓ | Same instance applied independently per series |
| ForecasterDirectMultiVariate | ✓ | Same instance applied independently per input series |
| ForecasterRecursiveClassifier | ✓ | Use `RollingFeaturesClassification` (not `RollingFeatures`) |
| ForecasterRnn | — | |
| ForecasterStats | — | |
| ForecasterEquivalentDate | — | |
| ForecasterFoundation | — | |

## How the forecaster uses `window_features`

When the forecaster builds the training matrix, it calls `wf.transform_batch(y)` on each instance once. The returned DataFrame is appended to the lag matrix as additional predictors. At prediction time, the forecaster slices the last `max_size_window_features` values of the running history and calls `wf.transform(...)` (numpy fast path) for each step.

**Effective `window_size`:**

```
window_size = max(max_lag, max_size_window_features) [+ differentiation]
```

Consequences:
- The first `window_size` rows of `y_train` cannot be used as targets (no full history available to compute predictors).
- `last_window` passed to `predict(..., last_window=...)` must contain at least `window_size` observations.
- A large rolling window combined with short series may leave too few training samples — check `len(y) - window_size` before fitting.

**Differentiation interaction:** the forecaster differentiates the series first, then computes window features on the differenced values. `roll_mean_7` with `differentiation=1` therefore averages seven *changes*, not seven raw observations. Inverse-differentiation only applies to the predictions, not to the predictors themselves.

**Multiple instances:** `window_features=[wf1, wf2]` calls `transform_batch` on each instance and concatenates the columns. Feature names must be unique across the combined set; duplicates silently collide. Use `RollingFeatures(features_names=[...])` to disambiguate.

**Lags-free mode:** `lags=None, window_features=...` is supported. `window_size` falls back to `max_size_window_features`. At least one of the two must be non-`None` — passing both as `None` raises `ValueError` from `set_lags` / `set_window_features`.

**Runtime mutation:** `forecaster.set_window_features(new_wf)` (analogous to `set_lags`) replaces the configuration after init and recomputes `window_size`, `window_features_names`, `window_features_class_names`, and the differentiator's window size. Useful in manual hyperparameter loops where the forecaster instance is reused across configurations.

## Custom window feature classes (protocol)

`window_features` accepts any object implementing the following contract — `RollingFeatures` is one such implementation, but user-defined classes are accepted as long as they conform.

**Required methods:**

| Method | Signature | Purpose |
|--------|-----------|---------|
| `transform_batch` | `(y: pd.Series) -> pd.DataFrame` | Computes features for the full training series. Used during `fit`. |
| `transform` | `(X: np.ndarray) -> np.ndarray` | Computes features for a single (or batched) prediction window. Used during `predict`. |

**Required attributes:**

| Attribute | Type | Purpose |
|-----------|------|---------|
| `max_window_size` | `int` | Largest window the class needs — drives the forecaster's `window_size`. |
| `features_names` | `list[str]` | Output column / feature names (must match `transform_batch` columns and the order produced by `transform`). |

**Validation enforced by the forecaster** (see `_create_window_features`):
- `transform_batch` must return a `pd.DataFrame` (`TypeError` otherwise).
- The DataFrame length must equal `len(y) - max_window_size` (`ValueError` otherwise).
- The DataFrame index must equal the forecaster's `train_index` (`ValueError` otherwise).

**Minimal custom example:**

```python
import numpy as np
import pandas as pd

class LastDeltaFeature:
    """Feature: difference between the last and the first value of the window."""

    def __init__(self, window_size: int):
        self.max_window_size = window_size
        self.features_names = [f'last_minus_first_{window_size}']

    def transform_batch(self, y: pd.Series) -> pd.DataFrame:
        # closed='left' to avoid leakage — last point excluded
        last = y.shift(1).rolling(self.max_window_size).apply(lambda w: w[-1], raw=True)
        first = y.shift(1).rolling(self.max_window_size).apply(lambda w: w[0], raw=True)
        out = (last - first).iloc[self.max_window_size:]
        return out.to_frame(name=self.features_names[0])

    def transform(self, X: np.ndarray) -> np.ndarray:
        # X shape: (window_length,) or (window_length, n_samples)
        if X.ndim == 1:
            return np.array([X[-1] - X[0]])
        return (X[-1, :] - X[0, :]).reshape(-1, 1)

# Use it like a RollingFeatures instance
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=7,
    window_features=LastDeltaFeature(window_size=14),
)
```

Combine custom and built-in instances freely:

```python
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=7,
    window_features=[
        RollingFeatures(stats=['mean', 'std'], window_sizes=14),
        LastDeltaFeature(window_size=14),
    ],
)
```

## Feature Selection Interaction

When using `select_features()`, the returned `selected_window_features` is a
**list of feature name strings** (e.g., `['roll_mean_7', 'roll_std_14']`),
NOT the `RollingFeatures` object. The original `RollingFeatures` instance
should still be passed to the forecaster.

```python
selected_lags, selected_wf, selected_exog = select_features(
    forecaster=forecaster,
    selector=RFECV(...),
    y=y_train,
)
# selected_wf = ['roll_mean_7', 'roll_std_14']  ← names, not objects
# The forecaster still uses the original RollingFeatures object
```

================================================================================
# SKILL: feature-selection
================================================================================

# Feature Selection

## When to Use

Use feature selection when:
- You have many lags or exogenous variables and want to reduce overfitting
- You want to identify which features matter most
- You need to speed up training by removing irrelevant features

### Related skills

- **Before**: `autocorrelation-and-lag-selection` (generate an informed candidate set of lags before running the selector)
- **Before**: `feature-engineering` (create the rolling, calendar, and exogenous features that the selector will rank)
- **After**: `hyperparameter-optimization` (tune the estimator on the reduced feature set)

## Single Series

`select_features` works with `ForecasterRecursive` and `ForecasterDirect`.

```python
from sklearn.feature_selection import RFECV
from sklearn.ensemble import RandomForestRegressor
from skforecast.recursive import ForecasterRecursive
from skforecast.preprocessing import RollingFeatures
from skforecast.feature_selection import select_features

# Create forecaster with many candidate features
rolling = RollingFeatures(stats=['mean', 'std', 'min', 'max'], window_sizes=[7, 14])

forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=48,  # Many lags — feature selection will reduce
    window_features=rolling,
)

# Run feature selection
selected_lags, selected_window_features, selected_exog = select_features(
    forecaster=forecaster,
    selector=RFECV(
        estimator=RandomForestRegressor(n_estimators=50, random_state=123),
        step=1,
        cv=3,
    ),
    y=y_train,
    exog=exog_train,
    select_only=None,          # 'autoreg' (lags+window), 'exog', or None (all)
    force_inclusion=None,      # Features to always keep (list or regex string)
    subsample=0.5,             # Use 50% of data for faster selection
    random_state=123,
    verbose=True,
)

# Apply selected lags to the same forecaster (simplest approach)
forecaster.set_lags(selected_lags)

# selected_window_features is a list of names (strings), not the RollingFeatures
# object. Use the names to verify which window features were selected.
print(f'Selected window features: {selected_window_features}')
print(f'Selected exog variables: {selected_exog}')
```

## Multi-Series

`select_features_multiseries` works with `ForecasterRecursiveMultiSeries` and `ForecasterDirectMultiVariate`.

> **Note:** When used with `ForecasterDirectMultiVariate`, `selected_lags` is returned as a `dict` (one entry per series) instead of a `list`.

```python
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.feature_selection import select_features_multiseries

forecaster = ForecasterRecursiveMultiSeries(
    estimator=RandomForestRegressor(n_estimators=100, random_state=123),
    lags=48,
    encoding='ordinal',
)

selected_lags, selected_window_features, selected_exog = select_features_multiseries(
    forecaster=forecaster,
    selector=RFECV(
        estimator=RandomForestRegressor(n_estimators=50, random_state=123),
        step=1,
        cv=3,
    ),
    series=series_df,
    exog=exog_df,
    select_only=None,
    force_inclusion=None,
    subsample=0.5,
    random_state=123,
    verbose=True,
)
```

## Force Inclusion

```python
# Always keep specific features regardless of selection
selected_lags, selected_wf, selected_exog = select_features(
    forecaster=forecaster,
    selector=selector,
    y=y_train,
    exog=exog_train,
    force_inclusion=['temperature', 'holiday'],  # Always keep these exog columns
)

# Regex pattern to force include
selected_lags, selected_wf, selected_exog = select_features(
    forecaster=forecaster,
    selector=selector,
    y=y_train,
    exog=exog_train,
    force_inclusion='^lag_',  # Keep all lag features
)
```

## Select Only Specific Feature Types

```python
# Only select among exogenous variables (keep all lags)
selected_lags, selected_wf, selected_exog = select_features(
    forecaster=forecaster,
    selector=selector,
    y=y_train,
    exog=exog_train,
    select_only='exog',  # Only select exog, keep all autoregressive features
)

# Only select among autoregressive features (keep all exog)
selected_lags, selected_wf, selected_exog = select_features(
    forecaster=forecaster,
    selector=selector,
    y=y_train,
    exog=exog_train,
    select_only='autoreg',  # Only select lags+window features, keep all exog
)
```

## Common Mistakes

1. **Using the wrong selector**: RFECV works best for recursive feature elimination. For faster selection, use `SelectFromModel`.
2. **Too small subsample**: If `subsample` is too small, selection may be unreliable. Use at least 0.3.
3. **Not updating forecaster**: After selection, update the forecaster with `forecaster.set_lags(selected_lags)` — the original is not modified in place by `select_features`.
4. **Running on full dataset**: Always run on training data only (`y_train`, `exog_train`).
5. **Confusing `selected_window_features` with `RollingFeatures`**: The returned `selected_window_features` is a list of **feature name strings** (e.g. `['mean_7', 'std_14']`), not the `RollingFeatures` object itself. Use these names to verify which window features were kept, but pass the original `RollingFeatures` instance to the forecaster.

================================================================================
# SKILL: drift-detection
================================================================================

# Drift Detection

## When to Use

Use drift detection to monitor whether new data falls outside the patterns seen during training. This helps decide when to retrain a model.

| Detector | Speed | Use Case |
|----------|-------|----------|
| `RangeDriftDetector` | Very fast | Real-time inference — checks if values are in training range |
| `PopulationDriftDetector` | Moderate | Batch monitoring — statistical tests for distribution shifts |

### Related skills

- **Before**: `forecasting-single-series` / `forecasting-multiple-series` (the detector is fitted on the training data of an existing forecaster)
- **Before**: `prediction-intervals` (intervals quantify uncertainty under the training distribution; drift detection flags when that distribution changes)

## RangeDriftDetector

Checks whether new observations fall within the ranges seen during training. Lightweight and suitable for real-time scoring.

> `fit()` accepts `series` and `exog` as a pandas Series, DataFrame, or dict
> (useful for multi-series pipelines with `ForecasterRecursiveMultiSeries`).

```python
from skforecast.drift_detection import RangeDriftDetector
from skforecast.recursive import ForecasterRecursive

# 1. Train the forecaster
forecaster = ForecasterRecursive(estimator=estimator, lags=24)
forecaster.fit(y=y_train, exog=exog_train)

# 2. Fit the drift detector on training data
detector = RangeDriftDetector()
detector.fit(series=y_train, exog=exog_train)

# 3. Check new data before making predictions
flag_drift, out_of_range_series, out_of_range_exog = detector.predict(
    last_window=new_data,
    exog=new_exog,
    verbose=True,           # Print drift details
    suppress_warnings=False,
)

if flag_drift:
    print("WARNING: New data contains values outside training range!")
    print(f"Out-of-range series features: {out_of_range_series}")
    print(f"Out-of-range exog features: {out_of_range_exog}")
```

## PopulationDriftDetector

Uses statistical tests to detect distribution shifts between reference (training) and new data.

> `fit(X)` and `predict(X)` expect a pandas DataFrame. For multi-series data,
> use a MultiIndex DataFrame with `(series_id, date)` index.

```python
from skforecast.drift_detection import PopulationDriftDetector

# 1. Create detector
detector = PopulationDriftDetector(
    chunk_size=100,                     # Split data into chunks of 100 obs
    threshold=3,                        # Multiplier for std deviation
    threshold_method='std',             # 'std' or 'quantile'
    max_out_of_range_proportion=0.1,    # Max 10% out-of-range allowed
)

# 2. Fit on reference (training) data
detector.fit(X=X_train)

# 3. Detect drift in new data
results, summary = detector.predict(X=X_new)

# results: DataFrame with per-chunk drift statistics
# summary: DataFrame with per-feature drift summary
print(summary)
```

### Chunk Size Options

```python
# Fixed number of observations per chunk
detector = PopulationDriftDetector(chunk_size=100)

# Time-based chunks
detector = PopulationDriftDetector(chunk_size='W')   # Weekly
detector = PopulationDriftDetector(chunk_size='M')   # Monthly
detector = PopulationDriftDetector(chunk_size='D')   # Daily

# No chunking — compare entire datasets
detector = PopulationDriftDetector(chunk_size=None)
```

### Threshold Methods

```python
# Standard deviation method: flag if statistic > mean + threshold * std
detector = PopulationDriftDetector(
    threshold=3,
    threshold_method='std',
)

# Quantile method: flag if statistic > empirical quantile
detector = PopulationDriftDetector(
    threshold=0.95,             # 95th percentile
    threshold_method='quantile',
)
```

## Integration with Forecasting Pipeline

```python
from skforecast.recursive import ForecasterRecursive
from skforecast.drift_detection import RangeDriftDetector

# Train
forecaster = ForecasterRecursive(estimator=estimator, lags=24)
forecaster.fit(y=y_train, exog=exog_train)

# Set up monitoring
detector = RangeDriftDetector()
detector.fit(series=y_train, exog=exog_train)

# Production loop
def predict_with_monitoring(new_window, new_exog):
    flag, _, _ = detector.predict(
        last_window=new_window, exog=new_exog, verbose=False
    )
    if flag:
        print("Drift detected — consider retraining the model")
    return forecaster.predict(steps=10, exog=new_exog)
```

## Common Mistakes

1. **Fitting detector on test data**: Always fit on training data — the reference distribution.
2. **Ignoring drift signals**: Drift doesn't mean the model is wrong, but it signals degradation risk.
3. **Over-sensitive thresholds**: Start with `threshold=3` (3 sigma) and adjust based on false positive rate.

================================================================================
# SKILL: deep-learning-forecasting
================================================================================

# Deep Learning Forecasting (RNN/LSTM)

## References

See [references/architecture-options.md](references/architecture-options.md) for
the complete `create_and_compile_model` signature, recurrent layer types,
output shape rules, exog architecture, custom Keras model requirements,
and `fit_kwargs` options.

## When to Use

Use `ForecasterRnn` when:
- You have large datasets (thousands of observations)
- Complex nonlinear patterns that tree-based models struggle with
- Multi-series problems where series share deep temporal patterns

**Requirements**: `pip install skforecast[deeplearning]` (installs keras)

### Related skills

- **Before**: `choosing-a-forecaster` (confirm `ForecasterRnn` is the right choice for the data size and pattern)
- **Before**: `feature-engineering` (RNN models still benefit from cyclical / calendar exogenous features)
- **After**: `hyperparameter-optimization` (tune RNN architecture and training hyperparameters)
- **After**: `prediction-intervals` (only conformal intervals are supported for `ForecasterRnn`)

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| `lags` in `ForecasterRnn` must match `create_and_compile_model(..., lags=...)` | Input shape mismatch error during fit | Use the same `lags` value in both calls |
| `ForecasterRnn` supports only `method='conformal'` for intervals | Error when calling `predict_interval(method='bootstrapping')` | Use `method='conformal'` |
| Pass `exog` to `create_and_compile_model` if exog is used in `fit()` / `predict()` | Architecture mismatch or failure when predicting with exog | Build the model with the same `exog` you train on |
| Scale inputs with `transformer_series=MinMaxScaler()` | Poor convergence; RNNs are scale-sensitive | Always set a scaler on `transformer_series` |

## Quick Start

```python
import pandas as pd
from skforecast.deep_learning import ForecasterRnn, create_and_compile_model
from sklearn.preprocessing import MinMaxScaler

# 1. Prepare data (DataFrame with DatetimeIndex, columns = series)
series = pd.read_csv('data.csv', index_col='date', parse_dates=True)
series = series.asfreq('h')

# 2. Create and compile a Keras model
model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels=series.columns.tolist(),  # All series
    recurrent_layer='LSTM',          # 'LSTM', 'GRU', or 'RNN'
    recurrent_units=[64, 32],        # Units per recurrent layer
    dense_units=[32],                # Units per dense layer
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
)

# 3. Create forecaster
forecaster = ForecasterRnn(
    levels=series.columns.tolist(),
    lags=48,
    estimator=model,
    transformer_series=MinMaxScaler(feature_range=(0, 1)),
    fit_kwargs={'epochs': 50, 'batch_size': 32, 'verbose': 0},
)

# 4. Train
forecaster.fit(series=series)

# 5. Predict
predictions = forecaster.predict(steps=24)
```

## Model Architecture with create_and_compile_model

```python
from skforecast.deep_learning import create_and_compile_model

# Simple LSTM
model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels='target',
    recurrent_layer='LSTM',
    recurrent_units=[64],
    dense_units=[32],
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
)

# Stacked LSTM (multiple recurrent layers)
model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels=series.columns.tolist(),
    recurrent_layer='LSTM',
    recurrent_units=[128, 64, 32],  # 3 stacked LSTM layers
    dense_units=[64, 32],           # 2 dense layers
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
)

# GRU variant (faster training)
model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels='target',
    recurrent_layer='GRU',
    recurrent_units=[64],
    dense_units=[32],
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
)

# Advanced: customize layer kwargs
model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels=series.columns.tolist(),
    recurrent_layer='LSTM',
    recurrent_units=[128, 64],
    recurrent_layers_kwargs={'activation': 'tanh'},   # default
    dense_units=[64],
    dense_layers_kwargs={'activation': 'relu'},        # default
    output_dense_layer_kwargs={'activation': 'linear'}, # default
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
    model_name='my_lstm_model',
)
```

## With Exogenous Variables

When using exogenous variables, pass `exog` to `create_and_compile_model` so it
builds the correct architecture (uses `TimeDistributed` layers internally).

```python
# exog must be a DataFrame covering the training period
exog = pd.DataFrame({'temperature': [...], 'holiday': [...]}, index=series.index)

model = create_and_compile_model(
    series=series,
    lags=48,
    steps=24,
    levels=series.columns.tolist(),
    exog=exog,                        # Passes exog info to build architecture
    recurrent_layer='LSTM',
    recurrent_units=[64, 32],
    dense_units=[32],
    compile_kwargs={'optimizer': 'adam', 'loss': 'mse'},
)

forecaster = ForecasterRnn(
    levels=series.columns.tolist(),
    lags=48,
    estimator=model,
    fit_kwargs={'epochs': 50, 'batch_size': 32, 'verbose': 0},
)

forecaster.fit(series=series, exog=exog)
predictions = forecaster.predict(steps=24, exog=exog_test)  # exog_test covers forecast horizon
```

## Custom Keras Model

```python
import keras

# Build your own model for full control (single level, no exog)
# Output units = steps * n_levels. For 1 level: steps. For N levels: steps * N + Reshape.
inputs = keras.layers.Input(shape=(48, 1))  # (lags, n_features)
x = keras.layers.LSTM(64, return_sequences=True)(inputs)
x = keras.layers.LSTM(32)(x)
x = keras.layers.Dense(32, activation='relu')(x)
outputs = keras.layers.Dense(24)(x)  # steps * n_levels (here 24 * 1)

model = keras.Model(inputs=inputs, outputs=outputs)
model.compile(optimizer='adam', loss='mse')

forecaster = ForecasterRnn(
    levels='target',
    lags=48,
    estimator=model,
    transformer_series=MinMaxScaler(feature_range=(0, 1)),
    fit_kwargs={'epochs': 100, 'batch_size': 32},
)
```

> **Multi-series custom model**: For N levels, the output layer should be
> `Dense(steps * n_levels)` followed by `Reshape((steps, n_levels))`.

## Prediction Intervals

```python
# ForecasterRnn supports conformal prediction only
forecaster.fit(series=series, store_in_sample_residuals=True)

predictions = forecaster.predict_interval(
    steps=24,
    method='conformal',           # Only 'conformal' supported
    interval=[0.1, 0.9],
    use_in_sample_residuals=True,
    use_binned_residuals=True,    # Better calibration with binned residuals
)
```

## Backtesting

```python
from skforecast.model_selection import backtesting_forecaster_multiseries, TimeSeriesFold

cv = TimeSeriesFold(
    steps=24,
    initial_train_size=len(series) - 200,
    refit=False,  # Retraining RNNs is expensive; set True only if needed
)

metric, predictions = backtesting_forecaster_multiseries(
    forecaster=forecaster,
    series=series,
    cv=cv,
    metric='mean_absolute_error',
)
```

## Common Mistakes

1. **Not scaling data**: RNNs are sensitive to scale. Always use `transformer_series=MinMaxScaler()`.
2. **Too few epochs**: Deep learning needs more training iterations. Start with 50-100 epochs.
3. **Wrong input shape**: The `lags` parameter in `ForecasterRnn` and `create_and_compile_model` must match.
4. **Refit=True in backtesting**: Retraining RNNs at every fold is very slow — use `refit=False` or `refit=5`.
5. **No GPU**: Training is slow on CPU. Use GPU if available.
6. **Using `predict_interval(method='bootstrapping')`**: ForecasterRnn only supports `method='conformal'`.
7. **Forgetting `exog` in `create_and_compile_model`**: If you use exog in `fit()`/`predict()`, you must also pass `exog` when building the model so the architecture accounts for the extra input features.

---
### Reference: architecture-options

# Deep Learning — Architecture Options Reference

## create_and_compile_model Signature

```python
from skforecast.deep_learning import create_and_compile_model

model = create_and_compile_model(
    series,                          # pd.DataFrame (required), input time series
    lags,                            # int | list[int] | np.ndarray | range (required)
    steps,                           # int (required), forecast horizon
    levels=None,                     # str | list[str] | None, output series (None = all)
    exog=None,                       # pd.Series | pd.DataFrame | None
    recurrent_layer='LSTM',          # 'LSTM' | 'GRU' | 'RNN'
    recurrent_units=100,             # int | list[int], units per recurrent layer
    recurrent_layers_kwargs={'activation': 'tanh'},   # dict | list[dict] | None
    dense_units=64,                  # int | list[int] | None
    dense_layers_kwargs={'activation': 'relu'},       # dict | list[dict] | None
    output_dense_layer_kwargs={'activation': 'linear'},  # dict | None
    compile_kwargs={'optimizer': Adam(), 'loss': MeanSquaredError()},  # dict
    model_name=None,                 # str | None
)
```

## Recurrent Layer Types

| Layer | Class | Speed | Memory | Best for |
|-------|-------|-------|--------|----------|
| `'LSTM'` | Long Short-Term Memory | Slowest | Highest | Long-range dependencies, default choice |
| `'GRU'` | Gated Recurrent Unit | Medium | Medium | Faster training, comparable performance |
| `'RNN'` | Simple RNN | Fastest | Lowest | Short sequences only, rarely used |

## Architecture Building Blocks

### Single recurrent layer

```python
model = create_and_compile_model(
    series=series, lags=48, steps=24, levels='target',
    recurrent_layer='LSTM',
    recurrent_units=64,           # single int → 1 recurrent layer
    dense_units=32,               # single int → 1 dense layer
)
# Architecture: Input → LSTM(64) → Dense(32) → Dense(24)
```

### Stacked recurrent layers

```python
model = create_and_compile_model(
    series=series, lags=48, steps=24, levels='target',
    recurrent_layer='LSTM',
    recurrent_units=[128, 64, 32],  # list → 3 stacked LSTM layers
    dense_units=[64, 32],            # list → 2 dense layers
)
# Architecture: Input → LSTM(128) → LSTM(64) → LSTM(32) → Dense(64) → Dense(32) → Dense(24)
```

### No dense layers

```python
model = create_and_compile_model(
    series=series, lags=48, steps=24, levels='target',
    recurrent_layer='GRU',
    recurrent_units=64,
    dense_units=None,               # None → no intermediate dense layers
)
# Architecture: Input → GRU(64) → Dense(24)  (output layer only)
```

## Output Layer Shape

The output layer is automatically sized based on `steps` and `levels`:

| Configuration | Output Dense units | Output shape |
|---------------|-------------------|--------------|
| 1 level, N steps | `steps` | `(batch, steps)` |
| M levels, N steps | `steps × M` + Reshape | `(batch, steps, M)` |

```python
# Single level
model = create_and_compile_model(
    series=series, lags=48, steps=24,
    levels='target',                 # 1 level → output: Dense(24)
)

# Multiple levels
model = create_and_compile_model(
    series=series, lags=48, steps=24,
    levels=['series_1', 'series_2', 'series_3'],  # 3 levels → output: Dense(72) + Reshape(24,3)
)
```

## Exogenous Variables Architecture

When `exog` is provided, the model creates a separate input branch with
`TimeDistributed` layers:

```python
# Without exog: single input branch
model = create_and_compile_model(
    series=series, lags=48, steps=24, levels='target',
)
# Input shape: (batch, 48, n_series)

# With exog: two input branches merged
model = create_and_compile_model(
    series=series, lags=48, steps=24, levels='target',
    exog=exog_df,                    # Must be provided at model creation
)
# Input shapes: series_input (batch, 48, n_series) + exog_input (batch, 48, n_exog)
```

> **Critical:** If you plan to use exog in `fit()` and `predict()`, you MUST
> pass `exog` to `create_and_compile_model()` so the architecture includes
> the exogenous input branch.

## Layer Kwargs Customization

### Same kwargs for all layers

```python
model = create_and_compile_model(
    ...,
    recurrent_units=[128, 64],
    recurrent_layers_kwargs={'activation': 'tanh', 'dropout': 0.2},  # applied to both
    dense_units=[64, 32],
    dense_layers_kwargs={'activation': 'relu'},  # applied to both
)
```

### Different kwargs per layer

```python
model = create_and_compile_model(
    ...,
    recurrent_units=[128, 64],
    recurrent_layers_kwargs=[
        {'activation': 'tanh', 'dropout': 0.3},  # first LSTM
        {'activation': 'tanh', 'dropout': 0.1},  # second LSTM
    ],
    dense_units=[64, 32],
    dense_layers_kwargs=[
        {'activation': 'relu'},
        {'activation': 'relu'},
    ],
)
```

## ForecasterRnn Constructor

```python
ForecasterRnn(
    levels,                          # str | list[str] (required)
    lags,                            # int | list[int] | np.ndarray | range (required)
    estimator=None,                  # Keras model (from create_and_compile_model or custom)
    transformer_series=MinMaxScaler(feature_range=(0, 1)),  # default: MinMaxScaler
    transformer_exog=MinMaxScaler(feature_range=(0, 1)),    # default: MinMaxScaler
    fit_kwargs=None,                 # dict, kwargs for model.fit()
    forecaster_id=None,              # str | int
)
```

### fit_kwargs

```python
forecaster = ForecasterRnn(
    levels=series.columns.tolist(),
    lags=48,
    estimator=model,
    fit_kwargs={
        'epochs': 50,                # number of training epochs
        'batch_size': 32,            # training batch size
        'verbose': 0,                # 0=silent, 1=progress bar, 2=one line per epoch
        'validation_split': 0.1,     # fraction for validation
        'callbacks': [EarlyStopping(patience=5)],  # Keras callbacks
    },
)
```

## Custom Keras Model Requirements

When building a custom model instead of using `create_and_compile_model`:

### Single level (no exog)

```python
import keras

inputs = keras.layers.Input(shape=(lags, n_series))  # (lags, number of input series)
x = keras.layers.LSTM(64, return_sequences=True)(inputs)
x = keras.layers.LSTM(32)(x)
x = keras.layers.Dense(32, activation='relu')(x)
outputs = keras.layers.Dense(steps)(x)               # units = steps

model = keras.Model(inputs=inputs, outputs=outputs)
model.compile(optimizer='adam', loss='mse')
```

### Multiple levels (no exog)

```python
inputs = keras.layers.Input(shape=(lags, n_series))
x = keras.layers.LSTM(64, return_sequences=True)(inputs)
x = keras.layers.LSTM(32)(x)
x = keras.layers.Dense(64, activation='relu')(x)
x = keras.layers.Dense(steps * n_levels)(x)           # units = steps × n_levels
outputs = keras.layers.Reshape((steps, n_levels))(x)   # reshape required

model = keras.Model(inputs=inputs, outputs=outputs)
model.compile(optimizer='adam', loss='mse')
```

### Input shape reference

| Scenario | Input shape | Output shape |
|----------|-------------|-------------|
| 1 level, no exog | `(lags, 1)` | `(steps,)` |
| M levels, no exog | `(lags, M)` | `(steps, M)` after Reshape |
| 1 level, K exog features | Two inputs: `(lags, 1)` + `(lags, K)` | `(steps,)` |
| M levels, K exog features | Two inputs: `(lags, M)` + `(lags, K)` | `(steps, M)` after Reshape |

## Prediction Intervals

ForecasterRnn only supports conformal prediction:

```python
forecaster.fit(series=series, store_in_sample_residuals=True)

predictions = forecaster.predict_interval(
    steps=24,
    method='conformal',           # only valid method
    interval=[0.1, 0.9],
    use_in_sample_residuals=True,
    use_binned_residuals=True,    # Better calibration with binned residuals
)
# NOTE: 'bootstrapping' is NOT supported for ForecasterRnn
# NOTE: predict_quantiles() and predict_dist() are NOT available
```

================================================================================
# SKILL: foundation-forecasting
================================================================================

# Foundation Model Forecasting (Zero-Shot)

## References

See [references/adapter-parameters.md](references/adapter-parameters.md) for the per-adapter constructor parameters of `ChronosAdapter`, `TimesFMAdapter`, `MoiraiAdapter`, `TabICLAdapter`, `TabPFNAdapter`, and `T0Adapter`.

## When to Use

Use `ForecasterFoundation` when:
- You want a **zero-shot baseline** before investing in model training.
- You have **very short histories** where ML models struggle.
- You need to forecast **cold-start** series (new product, new sensor).
- You want to compare against pre-trained generalist models.

Foundation models are **pre-trained on massive corpora** — `fit()` does not train them; it only stores the recent context and metadata.

## Stop Conditions

Scan before writing code. Each row lists a rule, the symptom when it is broken, and the recovery. Full pitfall catalog: the `troubleshooting-common-errors` skill.

| Rule | Symptom | Recovery |
|------|---------|----------|
| `fit()` stores context only; it never trains the model | Expecting training to happen or weights to update | Treat the model as pre-trained; evaluate with `backtesting_foundation` |
| Only Chronos-2, TabICL, TabPFN-TS, and T0 use `exog`; TimesFM 2.5 and Moirai-2 ignore it | `exog` silently dropped, no error raised | Pick an exog-capable adapter when covariates matter |
| TimesFM 2.5 and Moirai-2 restrict quantiles to `[0.1, 0.2, ..., 0.9]` | Requested quantile rejected or unsupported | Request only supported quantiles, or use an adapter allowing any quantile in (0, 1) |
| Each backend library must be installed separately | `ModuleNotFoundError` / `ImportError` on first use | `pip install` the matching backend (see Installation) |

## Installation

Foundation model backends are **not** bundled with skforecast. Install only the backend(s) you need:

```bash
pip install chronos-forecasting                                 # For Chronos-2
pip install git+https://github.com/google-research/timesfm.git  # For TimesFM 2.5
pip install uni2ts                                              # For Moirai-2
pip install tabicl[forecast]                                    # For TabICL
pip install tabpfn-time-series                                  # For TabPFN-TSpip install tfc-t0                                             # For T0```

Models are downloaded from HuggingFace on first use.

## Quick Start (single series)

```python
import pandas as pd
from skforecast.foundation import FoundationModel, ForecasterFoundation

# Data must have a DatetimeIndex with a frequency
data = pd.read_csv('data.csv', index_col='date', parse_dates=True).asfreq('h')

# 1. Configure a foundation model (adapter is resolved from model_id)
model = FoundationModel(
    model_id='autogluon/chronos-2-small',
    context_length=2048,      # Adapter-specific default: see reference
    device_map='auto',        # 'auto' picks CUDA > MPS > CPU
)

# 2. Wrap it in ForecasterFoundation for the skforecast API
forecaster = ForecasterFoundation(estimator=model)

# 3. "Fit" only stores the last context_length observations (no training)
forecaster.fit(series=data['target'])

# 4. Point forecast — returns long-format DataFrame: columns ['level', 'pred']
predictions = forecaster.predict(steps=24)
```

## Multi-Series (Global Zero-Shot Model)

Pass a wide `DataFrame`, a long-format `DataFrame` (MultiIndex), or a
`dict[str, pd.Series]` to `fit`.

```python
# series: wide DataFrame — each column is one series
forecaster.fit(series=series)

# Forecast all series
predictions = forecaster.predict(steps=24)

# Forecast a subset
predictions = forecaster.predict(steps=24, levels=['series_1', 'series_2'])
```

Chronos-2 supports `cross_learning=True` to share information across series
in the batch (ignored in single-series mode):

```python
model = FoundationModel(
    model_id='autogluon/chronos-2-small',
    cross_learning=True,
)
```

## With Exogenous Variables (Chronos-2, TabICL, TabPFN-TS and T0)

Chronos-2, TabICL, TabPFN-TS and T0 (`allow_exog=True`) accept exogenous variables. TimesFM 2.5 and Moirai-2 ignore them.

```python
# Historical + future exog (must cover the forecast horizon)
forecaster.fit(series=data['target'], exog=exog_train)

predictions = forecaster.predict(steps=24, exog=exog_test)
```

## Prediction Intervals and Quantiles

Foundation models output native quantile forecasts — no bootstrapping or conformal calibration is required.

```python
# Interval (lower/upper bounds from the model's quantiles)
predictions = forecaster.predict_interval(
    steps=24,
    interval=[0.1, 0.9],   # 80% prediction interval (quantiles, 0-1 scale)
)
# Columns: ['level', 'pred', 'lower_bound', 'upper_bound']

# Explicit quantiles
predictions = forecaster.predict_quantiles(
    steps=24,
    quantiles=[0.1, 0.5, 0.9],
)
# Columns: ['level', 'q_0.1', 'q_0.5', 'q_0.9']
```

For TimesFM 2.5 and Moirai-2, requested quantiles must be a subset of `[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]`. Chronos-2, TabICL, TabPFN-TS and TFC-T0 support any quantile in `(0, 1)`.

## Choosing a Model

| Model (`model_id` prefix)              | Exog | Default context | Best for                                          |
|----------------------------------------|:----:|----------------:|---------------------------------------------------|
| `autogluon/chronos-2-*` (Amazon)       | Yes  | 8192            | General-purpose, exog-friendly, cross-series info |
| `google/timesfm-2.5-*` (Google)        | No   | 512             | Long-horizon point/quantile forecasts             |
| `Salesforce/moirai-2.0-*` (Salesforce) | No   | 2048            | Multivariate pretraining, probabilistic forecasts |
| `soda-inria/tabicl` (Soda-INRIA)       | Yes  | 4096            | Tabular in-context learning, exog-aware           |
| `priorlabs/tabpfn-ts` (Prior Labs)     | Yes  | 32768           | Tabular foundation model, exog-aware, long context |
| `theforecastingcompany/t0` (TFC)       | Yes  | 8192            | Probabilistic forecasts, exog-aware (future covariates) |

The adapter is resolved automatically from the `model_id` prefix — no need to import adapter classes directly.

## Backtesting

Use the dedicated `backtesting_foundation` function — it is the only backtester that accepts a `ForecasterFoundation`. Refit is always disabled internally (the loaded model weights are preserved across folds) and probabilistic output is requested via `quantiles`, not `interval`.

```python
from skforecast.model_selection import backtesting_foundation, TimeSeriesFold

cv = TimeSeriesFold(
    steps=24,
    initial_train_size=len(series) - 200,
    refit=False,      # Refit is always disabled for foundation forecasters
)

metric, predictions = backtesting_foundation(
    forecaster=forecaster,
    series=series,
    cv=cv,
    metric='mean_absolute_error',
    quantiles=[0.1, 0.5, 0.9],   # Native model quantiles; no bootstrapping
)
```

## Override the Stored Context

Pass `context` at predict time to forecast from a different window without
refitting — useful for one-off predictions or custom backtesting loops:

```python
predictions = forecaster.predict(
    steps=24,
    context=new_window,        # pandas Series / DataFrame / dict
    context_exog=new_exog,     # Only with exog-aware adapters
    exog=future_exog,
)
```

If `context` is longer than the adapter's `context_length`, it is trimmed
automatically to the last `context_length` observations.

## Common Pitfalls

1. **Expecting `fit()` to train the model**: it only stores context. The weights come from HuggingFace.
2. **Index without frequency**: call `series.asfreq('h')` (or similar) before `fit` — skforecast requires a frequency.
3. **Passing `exog` to TimesFM 2.5 / Moirai-2**: ignored. Only Chronos-2, TabICL, TabPFN-TS and TFC-T0 support exogenous variables.
4. **Requesting unsupported quantiles**: TimesFM 2.5 and Moirai-2 are restricted to the nine deciles `0.1 … 0.9`.
5. **Large model downloads**: first call can be slow; consider using smaller variants (`*-small`) for experimentation.
6. **Forgetting to install the backend**: each foundation model requires its own library (`chronos-forecasting`, `timesfm`, `uni2ts`, `tabicl`, `tabpfn-time-series`, `tfc-t0`). Install only the one(s) you need.

---
### Reference: adapter-parameters

# Foundation Adapter Parameters

`FoundationModel` resolves the adapter automatically from `model_id`. All keyword arguments passed to `FoundationModel(...)` beyond `model_id` are forwarded to the chosen adapter's `__init__`.

```python
from skforecast.foundation import FoundationModel

model = FoundationModel(
    model_id='autogluon/chronos-2-small',
    context_length=2048,
    device_map='auto',
)
```

## ChronosAdapter — Amazon Chronos-2

- **`model_id` prefix**: `autogluon/chronos`
- **`allow_exog`**: `True` (past and future covariates)
- **Quantiles**: any value in `(0, 1)`

| Parameter        | Type    | Default  | Description                                                                    |
|------------------|---------|----------|--------------------------------------------------------------------------------|
| `model_id`       | str     | —        | HuggingFace model ID (e.g. `autogluon/chronos-2-small`).                       |
| `pipeline`       | object  | `None`   | Pre-loaded `BaseChronosPipeline`. If `None`, loaded lazily on first `predict`. |
| `context_length` | int     | `8192`   | Max historical observations kept as context.                                   |
| `predict_kwargs` | dict    | `None`   | Extra kwargs forwarded to the pipeline's `predict_quantiles`.                  |
| `device_map`     | str     | `'auto'` | Device placement: `'auto'` (CUDA > MPS > CPU), `'cuda'`, `'mps'`, `'cpu'`.     |
| `torch_dtype`    | object  | `None`   | Torch dtype for `from_pretrained` (e.g. `torch.bfloat16`).                     |
| `cross_learning` | bool    | `False`  | If `True`, shares information across series in multi-series batches.           |

## TimesFMAdapter — Google TimesFM 2.5

- **`model_id` prefix**: `google/timesfm`
- **`allow_exog`**: `False`
- **Supported quantiles**: `[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]`

| Parameter                | Type | Default | Description                                                         |
|--------------------------|------|---------|---------------------------------------------------------------------|
| `model_id`               | str  | —       | HuggingFace model ID (e.g. `google/timesfm-2.5-200m-pytorch`).      |
| `model`                  | obj  | `None`  | Pre-loaded & compiled TimesFM model. If `None`, loaded lazily.      |
| `context_length`         | int  | `512`   | Max historical observations kept as context.                        |
| `max_horizon`            | int  | `512`   | Max forecast horizon. `predict(steps=...)` must be ≤ this.          |
| `forecast_config_kwargs` | dict | `None`  | Extra kwargs forwarded to `timesfm.ForecastConfig` at compile time. |

The model is compiled lazily for the exact requested `steps` (up to `max_horizon`) to avoid unnecessary decode iterations.

## MoiraiAdapter — Salesforce Moirai-2

- **`model_id` prefix**: `Salesforce/moirai`
- **`allow_exog`**: `False`
- **Supported quantiles**: `[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]`

| Parameter        | Type | Default  | Description                                                              |
|------------------|------|----------|--------------------------------------------------------------------------|
| `model_id`       | str  | —        | HuggingFace model ID (e.g. `Salesforce/moirai-2.0-R-small`).             |
| `module`         | obj  | `None`   | Pre-loaded `Moirai2Module`. If `None`, loaded lazily.                    |
| `context_length` | int  | `2048`   | Max historical observations kept as context.                             |
| `device`         | str  | `'auto'` | Device placement: `'auto'` (CUDA > MPS > CPU), `'cuda'`, `'mps'`, `'cpu'`. |

## TabICLAdapter — Soda-INRIA TabICL

- **`model_id` prefix**: `soda-inria/tabicl`
- **`allow_exog`**: `True` (past and future covariates)
- **Quantiles**: any value in `(0, 1)`

| Parameter            | Type  | Default  | Description                                                                      |
|----------------------|-------|----------|----------------------------------------------------------------------------------|
| `model_id`           | str   | —        | HuggingFace model ID (e.g. `soda-inria/tabicl`).                                 |
| `model`              | obj   | `None`   | Pre-instantiated `TabICLForecaster`. If `None`, created lazily on first predict. |
| `context_length`     | int   | `4096`   | Max historical observations kept as context.                                     |
| `point_estimate`     | str   | `'mean'` | Point forecast method: `'mean'` or `'median'`.                                   |
| `tabicl_config`      | dict  | `None`   | Extra kwargs forwarded to `TabICLRegressor` at inference time.                   |
| `temporal_features`  | list  | `None`   | `TimeTransform` instances applied before inference. `None` = TabICL defaults; `[]` = disable all. |

## TabPFNAdapter — Prior Labs TabPFN-TS

- **`model_id` prefix**: `priorlabs/tabpfn`
- **`allow_exog`**: `True` (known-future covariates; covariates without future values are discarded by the library)
- **Quantiles**: any value in `(0, 1)`

| Parameter             | Type  | Default    | Description                                                                      |
|-----------------------|-------|------------|----------------------------------------------------------------------------------|
| `model_id`            | str   | —          | Model ID (e.g. `priorlabs/tabpfn-ts`). Used only for adapter resolution.         |
| `model`               | obj   | `None`     | Pre-instantiated `TabPFNTSPipeline`. If `None`, created lazily on first predict. |
| `context_length`      | int   | `32768`    | Max historical observations kept as context. Lower (e.g. 4096) for faster inference. |
| `mode`                | str   | `'local'`  | `'local'` (on-device inference, CUDA > MPS > CPU) or `'client'` (Prior Labs cloud API, no GPU needed). |
| `point_estimate`      | str   | `'median'` | Ensemble aggregation for the point forecast: `'mean'`, `'median'` or `'mode'`.   |
| `tabpfn_model_config` | dict  | `None`     | Extra config forwarded to the underlying TabPFN regressor (e.g. `model_path`, `device`). |
| `temporal_features`   | list  | `None`     | `FeatureGenerator` instances applied before inference. `None` = TabPFN-TS defaults; `[]` = disable all. |

## T0Adapter — The Forecasting Company T0

- **`model_id` prefix**: `theforecastingcompany/t0`
- **`allow_exog`**: `True` (future-known covariates; historical and known-future values are concatenated into the `[context + horizon]` covariate stream)
- **Quantiles**: any value in `(0, 1)` (native levels `0.1, 0.25, 0.5, 0.75, 0.9`; other levels are produced by inference-time interpolation)

| Parameter        | Type   | Default  | Description                                                                |
|------------------|--------|----------|----------------------------------------------------------------------------|
| `model_id`       | str    | —        | HuggingFace model ID (e.g. `theforecastingcompany/t0-alpha`).             |
| `model`          | obj    | `None`   | Pre-loaded `T0Forecaster`. If `None`, loaded lazily on first `predict`.    |
| `context_length` | int    | `8192`   | Max historical observations kept as context.                               |
| `device_map`     | str    | `'auto'` | Device placement: `'auto'` (CUDA > MPS > CPU), `'cuda'`, `'mps'`, `'cpu'`. |
| `torch_dtype`    | object | `None`   | Torch dtype the loaded model is cast to (e.g. `torch.bfloat16`).           |

Point forecasts use the median (quantile `0.5`). Covariates must be numeric; encode categoricals as numbers before passing them. A series with no future exog is forecast without covariates.

## Common Behavior

All adapters implement the same minimal interface:

- `fit(series, exog=None)` — stores context and metadata; no training.
- `predict(steps, context, context_exog, exog, quantiles)` — returns a   `dict[str, np.ndarray]` of shape `(steps, n_quantiles)` keyed by series name.
- `get_params()` / `set_params(**kwargs)` — sklearn-style parameter access.

Backend libraries (`chronos-forecasting`, `timesfm`, `uni2ts`, `tabicl`, `tabpfn-time-series`, `tfc-t0`) are imported **lazily** inside the adapter method that needs them, so only the backend for the adapter you actually use needs to be installed.

================================================================================
# SKILL: choosing-a-forecaster
================================================================================

# Choosing a Forecaster

## When to Use This Skill

Use this skill when the user needs help choosing a forecaster, comparing forecaster types (recursive vs direct, single vs multi-series), or understanding which skforecast class fits their problem.

### Related skills

- **After**: `forecasting-single-series` (apply the chosen forecaster to one target series)
- **After**: `forecasting-multiple-series` (apply the chosen forecaster to several series jointly)
- **After**: `autocorrelation-and-lag-selection` (analyse the series dynamics before configuring `lags`)
- **After**: `feature-engineering` (build the input feature set: calendar, rolling, exogenous)

## Overview

Skforecast is a **machine learning-first** library. The ML forecasters are the primary tools; statistical models (`ForecasterStats`) and naive baselines (`ForecasterEquivalentDate`) serve as comparison benchmarks.

## Step 1 — How Many Series?

| Scenario | Go to |
|----------|-------|
| **1 target series** (with or without exogenous variables) | → Step 2a |
| **Multiple series** to forecast simultaneously | → Step 2b |
| **Multiple series as drivers** to predict one target | → `ForecasterDirectMultiVariate` |
| **Categorical target** (e.g., low/medium/high) | → `ForecasterRecursiveClassifier` |

## Step 2a — Single Series

| Scenario | Recommended Forecaster | Why |
|----------|----------------------|-----|
| **General purpose** (start here) | `ForecasterRecursive` | Default choice. One model, recursive multi-step. Works with any sklearn-compatible estimator (LightGBM, XGBoost, CatBoost, RandomForest, etc.). Supports lags, window features, exog, differentiation, transformers, weight functions, and all probabilistic prediction methods (bootstrapping, conformal, quantiles, distributions) |
| **Horizon-dependent patterns** (e.g., predicting at 1h vs 24h requires different relationships) | `ForecasterDirect` | Trains one independent model per step — no error propagation. Better when the predictive relationship changes significantly across the forecast horizon. Requires `steps` at init; parallelizable with `n_jobs` |
| **Statistical baseline** | `ForecasterStats` | Wraps ARIMA, SARIMAX, ETS, ARAR. Use as a benchmark to compare against ML models, or when the series is very short (< 200 obs) and ML overfits |
| **Zero-shot / cold-start / no training data** | `ForecasterFoundation` | Wraps pre-trained foundation models (Chronos-2, TimesFM 2.5, Moirai-2, TabICL, TabPFN-TS, TFC-T0). `fit()` only stores context — no training. Good baseline and cold-start option. See the `foundation-forecasting` skill |
| **Naive baseline** | `ForecasterEquivalentDate` | Predicts using equivalent past dates (e.g., same weekday last week). Use as a sanity-check baseline |

## Step 2b — Multiple Series

| Scenario | Recommended Forecaster | Why |
|----------|----------------------|-----|
| **Forecast many series with a shared model** (start here) | `ForecasterRecursiveMultiSeries` | One global model learns cross-series patterns. Supports DataFrame or dict input (dict allows series with different date ranges). Encoding options: `'ordinal'` (default), `'ordinal_category'`, `'onehot'`, `None`. Supports per-series transformers, per-series differentiation, series_weights |
| **Other series are features for one target** | `ForecasterDirectMultiVariate` | All series become input features to predict a single `level`. Per-series lags via dict (`{'sales': [1,7], 'price': [1]}`). One model per step — no error propagation |
| **Deep learning / complex nonlinear patterns** | `ForecasterRnn` | Keras-based RNN/LSTM/GRU. Single model outputs all steps and levels simultaneously via 3D tensors. Only conformal intervals (no bootstrapping). Requires keras |
| **Zero-shot / pre-trained generalist** | `ForecasterFoundation` | Global zero-shot forecasts via Chronos-2 / TimesFM 2.5 / Moirai-2 / TabICL / TabPFN-TS / TFC-T0. `fit()` only stores context. Native quantile intervals. Chronos-2, TabICL, TabPFN-TS, and TFC-T0 support exog; TimesFM 2.5 & Moirai-2 do not. See the `foundation-forecasting` skill |

## Decision Flowchart

```
How many series?
│
├─► 1 series
│   │
│   ├─► Is it a classification problem? ──► Yes ──► ForecasterRecursiveClassifier
│   │
│   └─► Regression (continuous target)
│       │
│       ├─► Does the forecast relationship change across the horizon?
│       │   ├─► No / Unsure ──► ForecasterRecursive  ← START HERE
│       │   └─► Yes (step-specific patterns) ──► ForecasterDirect
│       │
│       └─► Compare with baselines:
│           ├─► ForecasterStats (Auto-ARIMA, ETS) as statistical benchmark
│           └─► ForecasterEquivalentDate as naive benchmark
│
└─► Multiple series
    │
    ├─► Want to forecast ALL series (global model)?
    │   └─► ForecasterRecursiveMultiSeries  ← START HERE
    │
    ├─► Want to use other series AS FEATURES for one target?
    │   └─► ForecasterDirectMultiVariate
    │
    └─► Need deep learning for very complex patterns?
        └─► ForecasterRnn

Zero-shot / no training data / cold-start (single or multi-series)?
    └─► ForecasterFoundation (Chronos-2 / TimesFM 2.5 / Moirai-2 / TabICL / TabPFN-TS / TFC-T0)
```

## Key Comparisons

### Recursive vs Direct (Single Series)

| Aspect | ForecasterRecursive | ForecasterDirect |
|--------|-------------------|-----------------|
| Models trained | 1 | N (one per step) |
| Error propagation | Yes (predictions feed into next step) | No (each step uses only observed data) |
| Forecast horizon | Flexible (any `steps` at predict time) | Fixed at init (`steps` required) |
| Training speed | Fast (one model) | Slower (N models, parallelizable via `n_jobs`) |
| Memory | Lower (1 estimator) | Higher (N estimators stored in `estimators_` dict) |
| Best for | Most cases, especially short-to-medium horizons | Long horizons where patterns change per step |
| Prediction intervals | Bootstrapping + conformal | Bootstrapping + conformal |

### RecursiveMultiSeries vs DirectMultiVariate (Multiple Series)

| Aspect | ForecasterRecursiveMultiSeries | ForecasterDirectMultiVariate |
|--------|-------------------------------|------------------------------|
| Goal | Forecast **all series** with one model | Use all series as **features** for one target |
| Input | DataFrame or dict of Series | DataFrame (all series as columns) |
| Target | All series (selected via `levels`) | One series (specified via `level`) |
| Strategy | Recursive (1 model, predictions feed back) | Direct (1 model per step, no error propagation) |
| Series identification | Encoding: `'ordinal'`, `'ordinal_category'`, `'onehot'`, `None` | All series create separate lag columns |
| Series with different ranges | Yes (via dict input) | No (all must share same range) |
| Per-series lags | No (same lags for all) | Yes (dict: `{'sales': [1,7], 'price': [1]}`) |
| Series weights | Yes (`series_weights` param) | No |
| Per-series transformers | Yes (`transformer_series` dict) | Yes (`transformer_series` dict) |
| Per-series differentiation | Yes (`differentiation` dict) | Yes (via `differentiator_` per series) |

### ML Forecasters vs Statistical Models

| Aspect | ML Forecasters | ForecasterStats |
|--------|---------------|-----------------|
| **Primary role** | Main forecasting tools | Comparison baseline |
| Best for | Medium-to-large datasets, complex patterns, many exogenous features | Short series, interpretability, parametric intervals |
| Estimators | Any sklearn-compatible (LightGBM, XGBoost, RF, etc.) | Arima, Sarimax, Ets, Arar |
| Lags / window features | Full support (`lags`, `RollingFeatures`) | No (model handles its own structure) |
| Differentiation | Built-in (`differentiation` param) | Handled within model (e.g., ARIMA `d` parameter) |
| Exogenous variables | Full support | Only SARIMAX |
| Prediction intervals | Bootstrapping (binned residuals) + conformal | Built-in parametric |
| Tuning | `grid_search_forecaster`, `random_search_forecaster`, `bayesian_search_forecaster` | `grid_search_stats`, `random_search_stats` |
| Backtesting | `backtesting_forecaster` / `backtesting_forecaster_multiseries` | `backtesting_stats` |
| Feature selection | `select_features` / `select_features_multiseries` | Not applicable |

## Feature Support Matrix

| Feature | Recursive | Direct | RecursiveMultiSeries | DirectMultiVariate | Rnn | Stats | EquivalentDate | Classifier |
|---------|:---------:|:------:|:-------------------:|:-----------------:|:---:|:-----:|:--------------:|:----------:|
| Lags | ✓ | ✓ | ✓ | ✓ (per-series dict) | ✓ | — | — | ✓ |
| Window features (`RollingFeatures`) | ✓ | ✓ | ✓ | ✓ | — | — | — | ✓ |
| Exogenous variables | ✓ | ✓ | ✓ | ✓ | ✓ | SARIMAX only | — | ✓ |
| Differentiation | ✓ | ✓ | ✓ (per-series) | ✓ (per-series) | — | Within model | — | — |
| Transformer y/series | ✓ | ✓ | ✓ (per-series) | ✓ (per-series) | ✓ | ✓ | — | — |
| Transformer exog | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | — | ✓ |
| Weight function | ✓ | ✓ | ✓ (per-series) | ✓ | — | — | — | ✓ |
| Bootstrapping intervals | ✓ | ✓ | ✓ | ✓ | — | — | — | — |
| Conformal intervals | ✓ | ✓ | ✓ | ✓ | ✓ | — | ✓ | — |
| Binned residuals | ✓ | ✓ | ✓ | ✓ | ✓ | — | ✓ | — |
| Quantile predictions | ✓ | ✓ | ✓ | ✓ | — | — | — | — |
| Distribution fitting | ✓ | ✓ | ✓ | ✓ | — | — | — | — |
| Class probabilities | — | — | — | — | — | — | — | ✓ |
| Feature importances | ✓ | ✓ | ✓ | ✓ | — | — | — | — |

> **Legend:** ✓ = supported, — = not supported/not applicable.

## Quick Start Recommendations

Once you have chosen a forecaster, follow these steps to get started:

1. **Define your problem**: 1 series → `ForecasterRecursive`; multiple series → `ForecasterRecursiveMultiSeries`
2. **Choose an estimator**: LightGBM (`LGBMRegressor`) is the best starting point — fast, handles categoricals, good defaults
3. **Add features**: Use `RollingFeatures` (rolling mean, std, min, max) and `CalendarFeatures` or `create_calendar_features` as exogenous variables
4. **Handle non-stationarity**: Use the `differentiation` parameter instead of manual differencing
5. **Evaluate with backtesting**: `backtesting_forecaster` + `TimeSeriesFold` for realistic multi-step evaluation
6. **Tune hyperparameters**: `bayesian_search_forecaster` (Optuna-based) — can include `lags` in the search space
7. **Add prediction intervals**: `predict_interval(method='bootstrapping', use_binned_residuals=True)` for uncertainty quantification
8. **Compare with baselines**: Use `ForecasterStats` (Auto-ARIMA: `Arima(order=None)`) and `ForecasterEquivalentDate` to verify the ML model adds value

================================================================================
# SKILL: troubleshooting-common-errors
================================================================================

# Troubleshooting Common Errors

This skill is the full pitfall catalog. The highest-risk skills also carry a focused `## Stop Conditions` table at the top (forecasting-single-series, forecasting-multiple-series, prediction-intervals, statistical-models, foundation-forecasting, feature-engineering, deep-learning-forecasting, hyperparameter-optimization); this skill remains the canonical source they point back to.

## Deprecated Import Paths

The most frequent LLM error. Old import paths no longer exist.

| Wrong (Deprecated) | Correct (v0.14.0+) |
|-------|---------|
| `from skforecast.ForecasterAutoreg import ForecasterAutoreg` | `from skforecast.recursive import ForecasterRecursive` |
| `from skforecast.ForecasterAutoregMultiSeries import ForecasterAutoregMultiSeries` | `from skforecast.recursive import ForecasterRecursiveMultiSeries` |
| `from skforecast.ForecasterAutoregDirect import ForecasterAutoregDirect` | `from skforecast.direct import ForecasterDirect` |
| `from skforecast.ForecasterAutoregMultiVariate import ForecasterAutoregMultiVariate` | `from skforecast.direct import ForecasterDirectMultiVariate` |
| `from skforecast.model_selection_multiseries import backtesting_forecaster_multiseries` | `from skforecast.model_selection import backtesting_forecaster_multiseries` |

## Wrong Class/Function Names

| Wrong | Correct |
|-------|---------|
| `ForecasterAutoreg` | `ForecasterRecursive` |
| `ForecasterAutoregMultiSeries` | `ForecasterRecursiveMultiSeries` |
| `ForecasterAutoregDirect` | `ForecasterDirect` |
| `ForecasterAutoregMultiVariate` | `ForecasterDirectMultiVariate` |
| `ForecasterSarimax` | `ForecasterStats(estimator=Sarimax(...))` |

## Removed Arguments

| Removed (v0.22.0+) | Replacement |
|---------------------|-------------|
| `regressor=...` | `estimator=...` (in all Forecasters) |

## Categorical Exogenous Variables

```python
# ❌ WRONG: setting categorical features directly on the estimator
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(categorical_feature=[0, 1]),
    lags=24,
)

# ✅ CORRECT: use categorical_features parameter on the forecaster
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    categorical_features='auto',  # or ['col_name_1', 'col_name_2']
)
```

## Data Issues

### "ValueError: The index of the series must be a DatetimeIndex with frequency"

```python
# Fix: set the frequency
data = data.asfreq('h')       # Hourly
data = data.asfreq('D')       # Daily
data = data.asfreq('MS')      # Monthly start
data = data.asfreq('QS')      # Quarterly start
```

### "ValueError: y contains NaN values"

```python
# Fix 1 (recommended for NaN-tolerant estimators): keep NaN rows
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(verbose=-1),  # LightGBM handles NaN natively
    lags=14,
    dropna_from_series=False,  # Default — NaN rows kept in training matrices
)
forecaster.fit(y=data['target'], suppress_warnings=True)

# Fix 2: drop rows with NaN from training matrices
forecaster = ForecasterRecursive(
    estimator=RandomForestRegressor(),
    lags=14,
    dropna_from_series=True,  # Drop NaN rows before fitting
)

# Fix 3: impute missing values before fitting
data = data.ffill()                      # Forward fill
data = data.interpolate(method='linear') # Linear interpolation
```

### "ValueError: exog must have the same index as y" / "exog does not cover forecast horizon"

```python
# Fix: exog for prediction must cover ALL future steps
# If predicting 10 steps ahead, exog_test must have at least 10 rows
# with dates matching the expected forecast dates
exog_test = exog.loc[forecast_start:forecast_end]
predictions = forecaster.predict(steps=10, exog=exog_test)
```

## Wrong Backtesting Function

```python
# ❌ WRONG: using backtesting_forecaster with ForecasterStats
from skforecast.model_selection import backtesting_forecaster
backtesting_forecaster(forecaster=forecaster_stats, y=y, cv=cv, metric=metric)  # Error!

# ✅ CORRECT: use backtesting_stats for statistical models
from skforecast.model_selection import backtesting_stats
backtesting_stats(forecaster=forecaster_stats, y=y, cv=cv, metric=metric)

# ❌ WRONG: using backtesting_forecaster with ForecasterRecursiveMultiSeries
backtesting_forecaster(forecaster=forecaster_multi, y=y, cv=cv, metric=metric)  # Error!

# ✅ CORRECT: use backtesting_forecaster_multiseries
from skforecast.model_selection import backtesting_forecaster_multiseries
backtesting_forecaster_multiseries(
    forecaster=forecaster_multi, series=series, cv=cv, metric=metric
)
```

## Wrong Search Function

```python
# ❌ WRONG: grid_search_forecaster with ForecasterStats
grid_search_forecaster(forecaster=forecaster_stats, y=y, cv=cv, param_grid=param_grid)

# ✅ CORRECT: grid_search_stats for statistical models
from skforecast.model_selection import grid_search_stats
grid_search_stats(forecaster=forecaster_stats, y=y, cv=cv, param_grid=param_grid)

# ❌ WRONG: grid_search_forecaster with ForecasterRecursiveMultiSeries
grid_search_forecaster(forecaster=forecaster_multi, y=y, cv=cv, param_grid=param_grid)

# ✅ CORRECT: grid_search_forecaster_multiseries
from skforecast.model_selection import grid_search_forecaster_multiseries
grid_search_forecaster_multiseries(
    forecaster=forecaster_multi, series=series, cv=cv, param_grid=param_grid
)
```

## Prediction Interval Errors

### "No in-sample residuals stored"

```python
# ❌ WRONG: fit without residuals, then call predict_interval
forecaster.fit(y=y_train)
forecaster.predict_interval(steps=10, method='bootstrapping')

# ✅ CORRECT: store residuals during fit
forecaster.fit(y=y_train, store_in_sample_residuals=True)
forecaster.predict_interval(steps=10, method='bootstrapping')
```

### Wrong interval method for a forecaster

| Forecaster | Supported Methods |
|------------|-------------------|
| `ForecasterRecursive` | `'bootstrapping'`, `'conformal'` |
| `ForecasterDirect` | `'bootstrapping'`, `'conformal'` |
| `ForecasterRecursiveMultiSeries` | `'bootstrapping'`, `'conformal'` (default: `'conformal'`) |
| `ForecasterDirectMultiVariate` | `'bootstrapping'`, `'conformal'` (default: `'conformal'`) |
| `ForecasterEquivalentDate` | `'conformal'` only |
| `ForecasterRnn` | `'conformal'` only |
| `ForecasterStats` | Built-in (uses `alpha` or `interval` parameter, no `method`) |
| `ForecasterRecursiveClassifier` | Not available — use `predict_proba()` |

## ETS Model API Confusion

```python
# ❌ WRONG (deprecated Ets API)
ets_model = Ets(error='add', trend='add', seasonal='add', seasonal_periods=12)

# ✅ CORRECT (current API)
ets_model = Ets(model='AAA', m=12)
# Model string: 1st char=Error, 2nd=Trend, 3rd=Seasonal
# A=Additive, M=Multiplicative, N=None, Z=Auto-select
```

## Function Mapping Reference

| Task | Single Series | Multi-Series | Statistical |
|------|--------------|-------------|-------------|
| **Backtesting** | `backtesting_forecaster` | `backtesting_forecaster_multiseries` | `backtesting_stats` |
| **Grid Search** | `grid_search_forecaster` | `grid_search_forecaster_multiseries` | `grid_search_stats` |
| **Random Search** | `random_search_forecaster` | `random_search_forecaster_multiseries` | `random_search_stats` |
| **Bayesian Search** | `bayesian_search_forecaster` | `bayesian_search_forecaster_multiseries` | N/A |
| **Feature Selection** | `select_features` | `select_features_multiseries` | N/A |

## Loading Serialized Forecasters from Older Versions

Forecasters saved (pickled/joblib) with older skforecast versions may fail to load or behave unexpectedly after upgrading. Internal attributes, class structures, and default values change between releases.

```python
# ❌ Common error when loading a forecaster saved with an older version
import joblib
forecaster = joblib.load('forecaster_v0.13.pkl')
# AttributeError: 'ForecasterRecursive' object has no attribute 'new_attribute'
# or: ModuleNotFoundError: No module named 'skforecast.ForecasterAutoreg'

# ✅ CORRECT: retrain the forecaster with the current version
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
)
forecaster.fit(y=y_train)
joblib.dump(forecaster, 'forecaster_v0.22.pkl')
```

**Best practices:**
- Always retrain and re-save forecasters after upgrading skforecast.
- Store training code (not just the serialized object) so models can be reproduced.
- Pin skforecast version in `requirements.txt` for production deployments.

================================================================================
# SKILL: complete-api-reference
================================================================================

# Complete API Reference

## When to Use This Skill

Use this when you need exact parameter names, types, defaults, or method signatures for any skforecast class or function.

## Overview

This skill contains the full constructor and method signatures for all
public skforecast classes and functions. See
[references/method-signatures.md](references/method-signatures.md) for
the complete reference, including:

- All forecaster constructors
- `fit()`, `predict()`, `predict_interval()`, `predict_quantiles()`, `predict_dist()` signatures
- `set_params()`, `set_lags()`, `set_out_sample_residuals()` signatures
- Method availability matrix (which forecaster supports which method)
- Backtesting, search, cross-validation, feature selection, and drift detection signatures

## Quick Index

### Forecaster Constructors
- `ForecasterRecursive` — single series, recursive strategy
- `ForecasterRecursiveMultiSeries` — multiple series, global model
- `ForecasterDirect` — single series, one model per step
- `ForecasterDirectMultiVariate` — multiple input series, one target
- `ForecasterRecursiveClassifier` — classification-based
- `ForecasterStats` — statistical models (ARIMA, ETS, SARIMAX, ARAR)
- `ForecasterEquivalentDate` — baseline using past offsets
- `ForecasterRnn` — deep learning (RNN/LSTM/GRU)
- `ForecasterFoundation` — zero-shot with foundation models (Chronos-2, TimesFM 2.5, Moirai-2, TabICL, TabPFN-TS, TFC-T0)
- `FoundationModel` — low-level foundation model wrapper used by `ForecasterFoundation`

### Forecaster Methods
- `fit()` — train the model
- `predict()` — generate point forecasts
- `predict_interval()` — generate prediction intervals

### Model Selection
- `backtesting_forecaster` — backtest single-series forecasters
- `backtesting_forecaster_multiseries` — backtest multi-series forecasters
- `backtesting_stats` — backtest statistical models
- `grid_search_forecaster` / `grid_search_forecaster_multiseries` / `grid_search_stats`
- `random_search_forecaster` / `random_search_forecaster_multiseries` / `random_search_stats`
- `bayesian_search_forecaster` / `bayesian_search_forecaster_multiseries`
- `TimeSeriesFold` — multi-step cross-validation
- `OneStepAheadFold` — fast one-step cross-validation

### Feature Selection
- `select_features` — single series
- `select_features_multiseries` — multi-series

### Drift Detection
- `RangeDriftDetector` — lightweight range check
- `PopulationDriftDetector` — statistical tests

### Preprocessing
- `RollingFeatures` — rolling window statistics
- `TimeSeriesDifferentiator` — differencing
- `CalendarFeatures` — calendar features

---
### Reference: method-signatures

# Method Signatures Reference

Complete constructor and method signatures for all skforecast public API.

## Forecaster Constructors

### ForecasterRecursive

```python
ForecasterRecursive(
    estimator=None,               # sklearn-compatible regressor
    lags=None,                    # int | list[int] | np.ndarray | range | None
    window_features=None,         # RollingFeatures | list[RollingFeatures] | None
    transformer_y=None,           # sklearn transformer for target variable
    transformer_exog=None,        # sklearn transformer | ColumnTransformer for exog
    categorical_features='auto',  # 'auto' | list[str] | None, categorical exog handling
    weight_func=None,             # Callable to weight training samples by index position
    differentiation=None,         # int, differencing order applied before training
    dropna_from_series=False,     # bool, drop NaN rows from training matrices
    fit_kwargs=None,              # dict, extra kwargs passed to estimator.fit()
    binner_kwargs=None,           # dict, kwargs for KBinsDiscretizer (binned residuals)
    forecaster_id=None,           # str | int, optional identifier
)
```

### ForecasterRecursiveMultiSeries

```python
ForecasterRecursiveMultiSeries(
    estimator=None,                # sklearn-compatible regressor
    lags=None,                     # int | list[int] | np.ndarray | range | None
    window_features=None,          # RollingFeatures | list[RollingFeatures] | None
    encoding='ordinal',            # 'ordinal' | 'ordinal_category' | 'onehot' | None
    transformer_series=None,       # sklearn transformer | dict[str, transformer] | None
    transformer_exog=None,         # sklearn transformer | ColumnTransformer | None
    categorical_features='auto',   # 'auto' | list[str] | None, categorical exog handling
    weight_func=None,              # Callable | dict[str, Callable] | None
    series_weights=None,           # dict[str, float] | None, relative weight of each series
    differentiation=None,          # int | dict[str, int | None] | None
    dropna_from_series=False,      # bool, allow NaN in individual series
    fit_kwargs=None,               # dict, extra kwargs passed to estimator.fit()
    binner_kwargs=None,            # dict, kwargs for KBinsDiscretizer (binned residuals)
    forecaster_id=None,            # str | int, optional identifier
)
```

### ForecasterDirect

```python
ForecasterDirect(
    steps,                        # int (required), number of steps to forecast
    estimator=None,               # sklearn-compatible regressor
    lags=None,                    # int | list[int] | np.ndarray | range | None
    window_features=None,         # RollingFeatures | list[RollingFeatures] | None
    transformer_y=None,           # sklearn transformer for target variable
    transformer_exog=None,        # sklearn transformer | ColumnTransformer for exog
    categorical_features='auto',  # 'auto' | list[str] | None, categorical exog handling
    weight_func=None,             # Callable to weight training samples by index position
    differentiation=None,         # int, differencing order applied before training
    dropna_from_series=False,     # bool, drop NaN rows from training matrices
    fit_kwargs=None,              # dict, extra kwargs passed to estimator.fit()
    binner_kwargs=None,           # dict, kwargs for KBinsDiscretizer (binned residuals)
    n_jobs='auto',                # int | str, parallel jobs for training one model per step
    forecaster_id=None,           # str | int, optional identifier
)
```

### ForecasterDirectMultiVariate

```python
ForecasterDirectMultiVariate(
    level,                         # str (required), name of the target series to predict
    steps,                         # int (required), number of steps to forecast
    estimator=None,                # sklearn-compatible regressor
    lags=None,                     # int | list | np.ndarray | range | dict[str, int|list] | None
    window_features=None,          # RollingFeatures | list[RollingFeatures] | None
    transformer_series=StandardScaler(),  # sklearn transformer | dict[str, transformer] | None
    transformer_exog=None,         # sklearn transformer | ColumnTransformer | None
    categorical_features='auto',   # 'auto' | list[str] | None, categorical exog handling
    weight_func=None,              # Callable to weight training samples by index position
    differentiation=None,          # int, differencing order applied before training
    dropna_from_series=False,      # bool, drop NaN rows from training matrices
    fit_kwargs=None,               # dict, extra kwargs passed to estimator.fit()
    binner_kwargs=None,            # dict, kwargs for KBinsDiscretizer (binned residuals)
    n_jobs='auto',                 # int | str, parallel jobs for training one model per step
    forecaster_id=None,            # str | int, optional identifier
)
```

### ForecasterRecursiveClassifier

```python
ForecasterRecursiveClassifier(
    estimator,                    # sklearn-compatible classifier (required, not optional)
    lags=None,                    # int | list[int] | np.ndarray | range | None
    window_features=None,         # RollingFeatures | list[RollingFeatures] | None
    features_encoding='auto',     # str, encoding for categorical exog features
    transformer_exog=None,        # sklearn transformer | ColumnTransformer for exog
    categorical_features='auto',  # 'auto' | list[str] | None, categorical exog handling
    weight_func=None,             # Callable to weight training samples by index position
    dropna_from_series=False,     # bool, drop NaN rows from training matrices
    fit_kwargs=None,              # dict, extra kwargs passed to estimator.fit()
    forecaster_id=None,           # str | int, optional identifier
)
# NOTE: No transformer_y, differentiation, or binner_kwargs.
# NOTE: Uses predict_proba() instead of predict_interval().
```

### ForecasterStats

```python
ForecasterStats(
    estimator=None,           # Arima | Sarimax | Ets | Arar | list of these
    transformer_y=None,       # sklearn transformer for target variable
    transformer_exog=None,    # sklearn transformer | ColumnTransformer for exog
    forecaster_id=None,       # str | int, optional identifier
)
```

### ForecasterEquivalentDate

```python
ForecasterEquivalentDate(
    offset,                   # int | pd.tseries.offsets.DateOffset (required)
    n_offsets=1,              # int, number of past offsets to aggregate
    agg_func=np.mean,         # Callable, function to aggregate multiple offsets
    binner_kwargs=None,       # dict, kwargs for KBinsDiscretizer (binned residuals)
    forecaster_id=None,       # str | int, optional identifier
)
```

### ForecasterRnn

```python
ForecasterRnn(
    estimator=None,                    # Keras model (use create_and_compile_model)
    levels,                            # str | list[str] (required), target series names
    lags,                              # int | list[int] | np.ndarray | range (required)
    transformer_series=MinMaxScaler(feature_range=(0, 1)),  # transformer | dict | None
    transformer_exog=MinMaxScaler(feature_range=(0, 1)),    # transformer | None
    fit_kwargs=None,                   # dict, extra kwargs passed to model.fit()
    binner_kwargs=None,                # dict, kwargs for KBinsDiscretizer (binned residuals)
    forecaster_id=None,                # str | int, optional identifier
)
```

### ForecasterFoundation

```python
FoundationModel(
    model_id,                  # str (required), e.g. 'autogluon/chronos-2-small'
    **kwargs,                  # Forwarded to the resolved adapter. Common keys:
                               #   context_length : int
                               #   device_map / device : 'auto' | 'cuda' | 'mps' | 'cpu'
                               #   torch_dtype : object (Chronos-2, T0)
                               #   cross_learning : bool (Chronos-2 only)
                               #   max_horizon, forecast_config_kwargs (TimesFM 2.5)
                               #   point_estimate, tabicl_config, temporal_features (TabICL)
                               #   mode, point_estimate, tabpfn_model_config, temporal_features (TabPFN-TS)
                               #   (T0 uses only context_length, device_map, torch_dtype)
)

ForecasterFoundation(
    estimator,                 # FoundationModel (required)
    forecaster_id=None,        # str | int, optional identifier
)
```

## Forecaster Methods: fit()

```python
# ForecasterRecursive, ForecasterDirect
forecaster.fit(
    y,                                # pd.Series with DatetimeIndex (required)
    exog=None,                        # pd.Series | pd.DataFrame | None
    store_last_window=True,           # bool
    store_in_sample_residuals=False,  # bool, set True before using predict_interval()
    random_state=123,                 # int, seed for residual sampling
    suppress_warnings=False           # bool
)

# ForecasterRecursiveMultiSeries
forecaster.fit(
    series,                           # pd.DataFrame | dict[str, pd.Series|pd.DataFrame] (required)
    exog=None,                        # pd.Series | pd.DataFrame | dict[str, pd.Series|pd.DataFrame] | None
    store_last_window=True,           # bool | list[str], True stores all, list stores specific series
    store_in_sample_residuals=False,  # bool, set True before using predict_interval()
    random_state=123,                 # int
    suppress_warnings=False           # bool
)

# ForecasterDirectMultiVariate, ForecasterRnn
forecaster.fit(
    series,                           # pd.DataFrame with multiple columns (required)
    exog=None,                        # pd.Series | pd.DataFrame | None
    store_last_window=True,           # bool
    store_in_sample_residuals=False,  # bool
    random_state=123,                 # int
    suppress_warnings=False           # bool
)

# ForecasterRecursiveClassifier
forecaster.fit(
    y,                                # pd.Series with DatetimeIndex (required)
    exog=None,                        # pd.Series | pd.DataFrame | None
    store_last_window=True,           # bool
    suppress_warnings=False           # bool
)
# NOTE: No store_in_sample_residuals or random_state.

# ForecasterStats
forecaster.fit(
    y,                                # pd.Series with DatetimeIndex (required)
    exog=None,                        # pd.Series | pd.DataFrame | None
    store_last_window=True,           # bool
    suppress_warnings=False           # bool
)

# ForecasterEquivalentDate
forecaster.fit(
    y,                                # pd.Series with DatetimeIndex (required)
    store_in_sample_residuals=False,  # bool
    random_state=123,                 # int
    suppress_warnings=False           # bool
)
# NOTE: No exog parameter (uses date offsets, not exogenous variables).

# ForecasterFoundation
forecaster.fit(
    series,                           # pd.Series | pd.DataFrame | dict[str, pd.Series] (required)
    exog=None,                        # pd.Series | pd.DataFrame | dict | None (Chronos-2 only)
)
# NOTE: "fit" does not train the model — it only stores the last
# context_length observations and metadata. Foundation models are
# pre-trained; training happens upstream on HuggingFace.
```

## Forecaster Methods: predict()

```python
# ForecasterRecursive
forecaster.predict(
    steps,                    # int | str | pd.Timestamp (required)
    last_window=None,         # pd.Series | pd.DataFrame | None
    exog=None,                # pd.Series | pd.DataFrame | None
    check_inputs=True,        # bool
    suppress_warnings=False   # bool
) -> pd.Series

# ForecasterRecursiveMultiSeries
forecaster.predict(
    steps,                    # int (required)
    levels=None,              # str | list[str] | None, which series to predict
    last_window=None,         # pd.DataFrame | None
    exog=None,                # pd.Series | pd.DataFrame | dict | None
    suppress_warnings=False,  # bool
    check_inputs=True         # bool
) -> pd.DataFrame

# ForecasterDirect
forecaster.predict(
    steps=None,               # int | list[int] | None, subset of trained steps
    last_window=None,         # pd.Series | pd.DataFrame | None
    exog=None,                # pd.Series | pd.DataFrame | None
    check_inputs=True,        # bool
    suppress_warnings=False   # bool
) -> pd.Series

# ForecasterDirectMultiVariate
forecaster.predict(
    steps=None,               # int | list[int] | None, subset of trained steps
    last_window=None,         # pd.DataFrame | None
    exog=None,                # pd.Series | pd.DataFrame | None
    suppress_warnings=False,  # bool
    check_inputs=True         # bool
) -> pd.DataFrame

# ForecasterRecursiveClassifier
forecaster.predict(
    steps,                    # int | str | pd.Timestamp (required)
    last_window=None,         # pd.Series | pd.DataFrame | None
    exog=None                 # pd.Series | pd.DataFrame | None
) -> pd.Series
# Also: predict_proba(steps, last_window=None, exog=None) -> pd.DataFrame

# ForecasterStats
forecaster.predict(
    steps,                    # int (required)
    last_window=None,         # pd.Series | None
    last_window_exog=None,    # pd.Series | pd.DataFrame | None, exog for last_window period
    exog=None,                # pd.Series | pd.DataFrame | None, exog for forecast period
    suppress_warnings=False   # bool
) -> pd.Series | pd.DataFrame

# ForecasterEquivalentDate
forecaster.predict(
    steps,                    # int (required)
    last_window=None,         # pd.Series | None
    check_inputs=True,        # bool
    suppress_warnings=False   # bool
) -> pd.Series

# ForecasterRnn
forecaster.predict(
    steps=None,               # int | list[int] | None, subset of trained steps
    levels=None,              # str | list[str] | None, which series to predict
    last_window=None,         # pd.DataFrame | None
    exog=None,                # pd.Series | pd.DataFrame | None
    suppress_warnings=False,  # bool
    check_inputs=True         # bool
) -> pd.DataFrame

# ForecasterFoundation
forecaster.predict(
    steps,                    # int (required)
    levels=None,              # str | list[str] | None, subset of series
    context=None,             # pd.Series | pd.DataFrame | dict | None, override stored context
    context_exog=None,        # pd.Series | pd.DataFrame | dict | None, historical exog
    exog=None,                # pd.Series | pd.DataFrame | dict | None, future exog (Chronos-2 only)
    check_inputs=True         # bool
) -> pd.DataFrame             # Long-format: columns ['level', 'pred']
# Also:
#   predict_interval(steps, ..., interval=[0.1, 0.9]) -> ['level','pred','lower_bound','upper_bound']
#   predict_quantiles(steps, ..., quantiles=[0.1, 0.5, 0.9]) -> ['level','q_0.1','q_0.5','q_0.9']
```

## Forecaster Methods: predict_interval()

```python
# ForecasterRecursive
forecaster.predict_interval(
    steps,                              # int | str | pd.Timestamp (required)
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    method='bootstrapping',             # 'bootstrapping' | 'conformal'
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    n_boot=250,                         # int, number of bootstrap samples
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterRecursiveMultiSeries (NOTE: default method='conformal')
forecaster.predict_interval(
    steps,                              # int (required)
    levels=None,                        # str | list[str] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | dict | None
    method='conformal',                 # 'bootstrapping' | 'conformal'
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirect
forecaster.predict_interval(
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    method='bootstrapping',             # 'bootstrapping' | 'conformal'
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirectMultiVariate (NOTE: default method='conformal')
forecaster.predict_interval(
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    method='conformal',                 # 'bootstrapping' | 'conformal'
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterStats (NOTE: different interface — uses alpha, no method/n_boot)
forecaster.predict_interval(
    steps,                              # int (required)
    last_window=None,                   # pd.Series | None
    last_window_exog=None,              # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    alpha=0.05,                         # float, significance level
    interval=None,                      # list[float] | tuple[float] | None, quantiles 0-1
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterEquivalentDate (NOTE: only 'conformal' method supported)
forecaster.predict_interval(
    steps,                              # int (required)
    last_window=None,                   # pd.Series | None
    method='conformal',                 # only 'conformal' supported
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=None,                  # Any, ignored (API compatibility)
    exog=None,                          # Any, ignored (API compatibility)
    n_boot=None,                        # Any, ignored (API compatibility)
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterRnn (NOTE: only 'conformal' method supported)
forecaster.predict_interval(
    steps=None,                         # int | list[int] | None
    levels=None,                        # str | list[str] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    method='conformal',                 # only 'conformal' supported
    interval=[0.05, 0.95],              # float (coverage) | list[float] | tuple[float], quantiles 0-1
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    suppress_warnings=False,            # bool
    n_boot=None,                        # Any, ignored (API compatibility)
    random_state=None,                  # Any, ignored (API compatibility)
) -> pd.DataFrame

# ForecasterRecursiveClassifier: No predict_interval(). Use predict_proba() instead.
```

## Backtesting Functions

```python
backtesting_forecaster(
    forecaster,                         # ForecasterRecursive | ForecasterDirect |
                                        # ForecasterEquivalentDate | ForecasterRecursiveClassifier
    y,                                  # pd.Series with DatetimeIndex
    cv,                                 # TimeSeriesFold
    metric,                             # str | Callable | list[str | Callable]
    exog=None,                          # pd.Series | pd.DataFrame | None
    interval=None,                      # float | list[float] | tuple[float] | str | distribution | None
    interval_method='bootstrapping',    # 'bootstrapping' | 'conformal'
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    return_predictors=False,            # bool
    n_jobs='auto',                      # int | str
    verbose=False,                      # bool
    show_progress=True,                 # bool
    suppress_warnings=False             # bool
) -> tuple[pd.DataFrame, pd.DataFrame]

backtesting_forecaster_multiseries(
    forecaster,                         # ForecasterRecursiveMultiSeries |
                                        # ForecasterDirectMultiVariate | ForecasterRnn
    series,                             # pd.DataFrame | dict[str, pd.Series | pd.DataFrame]
    cv,                                 # TimeSeriesFold
    metric,                             # str | Callable | list[str | Callable]
    levels=None,                        # str | list[str] | None
    add_aggregated_metric=True,         # bool
    exog=None,                          # pd.Series | pd.DataFrame | dict | None
    interval=None,                      # float | list[float] | tuple[float] | str | distribution | None
    interval_method='conformal',        # 'bootstrapping' | 'conformal' (NOTE: default 'conformal')
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    return_predictors=False,            # bool
    n_jobs='auto',                      # int | str
    verbose=False,                      # bool
    show_progress=True,                 # bool
    suppress_warnings=False             # bool
) -> tuple[pd.DataFrame, pd.DataFrame]

backtesting_stats(
    forecaster,                         # ForecasterStats
    y,                                  # pd.Series with DatetimeIndex
    cv,                                 # TimeSeriesFold
    metric,                             # str | Callable | list[str | Callable]
    exog=None,                          # pd.Series | pd.DataFrame | None
    alpha=None,                         # float | None, significance level
    interval=None,                      # list[float] | tuple[float] | None
    freeze_params=True,                 # bool, if True only first fold fits the model
    n_jobs='auto',                      # int | str
    verbose=False,                      # bool
    show_progress=True,                 # bool
    suppress_warnings=False             # bool
) -> tuple[pd.DataFrame, pd.DataFrame]
```

## Hyperparameter Search Functions

### Single Series

```python
grid_search_forecaster(
    forecaster,              # ForecasterRecursive | ForecasterDirect
    y,                       # pd.Series with DatetimeIndex
    cv,                      # TimeSeriesFold | OneStepAheadFold
    param_grid,              # dict, sklearn-style parameter grid
    metric,                  # str | Callable | list[str | Callable]
    exog=None,               # pd.Series | pd.DataFrame | None
    lags_grid=None,          # list[int | list | np.ndarray | range] | dict | None
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None, path to save results incrementally
) -> pd.DataFrame

random_search_forecaster(
    forecaster,              # ForecasterRecursive | ForecasterDirect
    y,                       # pd.Series with DatetimeIndex
    cv,                      # TimeSeriesFold | OneStepAheadFold
    param_distributions,     # dict, parameter distributions for sampling
    metric,                  # str | Callable | list[str | Callable]
    exog=None,               # pd.Series | pd.DataFrame | None
    lags_grid=None,          # list[int | list | np.ndarray | range] | dict | None
    n_iter=10,               # int, number of random parameter combinations
    random_state=123,        # int
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None
) -> pd.DataFrame

bayesian_search_forecaster(
    forecaster,              # ForecasterRecursive | ForecasterDirect
    y,                       # pd.Series with DatetimeIndex
    cv,                      # TimeSeriesFold | OneStepAheadFold
    search_space,            # Callable, Optuna trial search space function
    metric,                  # str | Callable | list[str | Callable]
    exog=None,               # pd.Series | pd.DataFrame | None
    n_trials=20,             # int, number of Optuna trials
    random_state=123,        # int
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None,        # str | None
    kwargs_create_study=None,     # dict | None, kwargs for optuna.create_study()
    kwargs_study_optimize=None    # dict | None, kwargs for study.optimize()
) -> tuple[pd.DataFrame, object]
```

### Multi-Series

```python
grid_search_forecaster_multiseries(
    forecaster,              # ForecasterRecursiveMultiSeries | ForecasterDirectMultiVariate | ForecasterRnn
    series,                  # pd.DataFrame | dict[str, pd.Series | pd.DataFrame]
    cv,                      # TimeSeriesFold | OneStepAheadFold
    param_grid,              # dict
    metric,                  # str | Callable | list[str | Callable]
    aggregate_metric=['weighted_average', 'average', 'pooling'],  # str | list[str]
    levels=None,             # str | list[str] | None
    exog=None,               # pd.Series | pd.DataFrame | dict | None
    lags_grid=None,          # list | dict | None
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None
) -> pd.DataFrame

random_search_forecaster_multiseries(
    forecaster,              # ForecasterRecursiveMultiSeries | ForecasterDirectMultiVariate | ForecasterRnn
    series,                  # pd.DataFrame | dict[str, pd.Series | pd.DataFrame]
    cv,                      # TimeSeriesFold | OneStepAheadFold
    param_distributions,     # dict
    metric,                  # str | Callable | list[str | Callable]
    aggregate_metric=['weighted_average', 'average', 'pooling'],  # str | list[str]
    levels=None,             # str | list[str] | None
    exog=None,               # pd.Series | pd.DataFrame | dict | None
    lags_grid=None,          # list | dict | None
    n_iter=10,               # int
    random_state=123,        # int
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None
) -> pd.DataFrame

bayesian_search_forecaster_multiseries(
    forecaster,              # ForecasterRecursiveMultiSeries | ForecasterDirectMultiVariate | ForecasterRnn
    series,                  # pd.DataFrame | dict[str, pd.Series | pd.DataFrame]
    cv,                      # TimeSeriesFold | OneStepAheadFold
    search_space,            # Callable, Optuna trial search space function
    metric,                  # str | Callable | list[str | Callable]
    aggregate_metric=['weighted_average', 'average', 'pooling'],  # str | list[str]
    levels=None,             # str | list[str] | None
    exog=None,               # pd.Series | pd.DataFrame | dict | None
    n_trials=20,             # int
    random_state=123,        # int
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None,        # str | None
    kwargs_create_study=None,     # dict | None
    kwargs_study_optimize=None    # dict | None
) -> tuple[pd.DataFrame, object]
```

### Statistical Models

```python
grid_search_stats(
    forecaster,              # ForecasterStats
    y,                       # pd.Series with DatetimeIndex
    cv,                      # TimeSeriesFold
    param_grid,              # dict
    metric,                  # str | Callable | list[str | Callable]
    exog=None,               # pd.Series | pd.DataFrame | None
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None
) -> pd.DataFrame

random_search_stats(
    forecaster,              # ForecasterStats
    y,                       # pd.Series with DatetimeIndex
    cv,                      # TimeSeriesFold
    param_distributions,     # dict
    metric,                  # str | Callable | list[str | Callable]
    exog=None,               # pd.Series | pd.DataFrame | None
    n_iter=10,               # int
    random_state=123,        # int
    return_best=True,        # bool
    n_jobs='auto',           # int | str
    verbose=False,           # bool
    show_progress=True,      # bool
    suppress_warnings=False, # bool
    output_file=None         # str | None
) -> pd.DataFrame
```

## Forecaster Methods: predict_quantiles()

```python
# ForecasterRecursive
forecaster.predict_quantiles(
    steps,                              # int | str | pd.Timestamp (required)
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    quantiles=[0.05, 0.5, 0.95],       # list[float] | tuple[float]
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterRecursiveMultiSeries
forecaster.predict_quantiles(
    steps,                              # int (required)
    levels=None,                        # str | list[str] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | dict | None
    quantiles=[0.05, 0.5, 0.95],       # list[float] | tuple[float]
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirect
forecaster.predict_quantiles(
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    quantiles=[0.05, 0.5, 0.95],       # list[float] | tuple[float]
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirectMultiVariate
forecaster.predict_quantiles(
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    quantiles=[0.05, 0.5, 0.95],       # list[float] | tuple[float]
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False,            # bool
    levels=None,                        # Any, ignored (API compatibility)
) -> pd.DataFrame

# NOT available in: ForecasterRecursiveClassifier, ForecasterStats,
#                    ForecasterEquivalentDate, ForecasterRnn
```

## Forecaster Methods: predict_dist()

```python
# ForecasterRecursive
forecaster.predict_dist(
    steps,                              # int | str | pd.Timestamp (required)
    distribution,                       # scipy.stats distribution object (required)
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterRecursiveMultiSeries
forecaster.predict_dist(
    steps,                              # int (required)
    distribution,                       # scipy.stats distribution object (required)
    levels=None,                        # str | list[str] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | dict | None
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirect
forecaster.predict_dist(
    distribution,                       # scipy.stats distribution object (required)
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.Series | pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False             # bool
) -> pd.DataFrame

# ForecasterDirectMultiVariate
forecaster.predict_dist(
    distribution,                       # scipy.stats distribution object (required)
    steps=None,                         # int | list[int] | None
    last_window=None,                   # pd.DataFrame | None
    exog=None,                          # pd.Series | pd.DataFrame | None
    n_boot=250,                         # int
    use_in_sample_residuals=True,       # bool
    use_binned_residuals=True,          # bool
    random_state=123,                   # int
    suppress_warnings=False,            # bool
    levels=None,                        # Any, ignored (API compatibility)
) -> pd.DataFrame

# NOT available in: ForecasterRecursiveClassifier, ForecasterStats,
#                    ForecasterEquivalentDate, ForecasterRnn
```

## Forecaster Methods: set_out_sample_residuals()

```python
# ForecasterRecursive, ForecasterDirect
forecaster.set_out_sample_residuals(
    y_true,                  # np.ndarray | pd.Series (required)
    y_pred,                  # np.ndarray | pd.Series (required)
    append=False,            # bool, append to existing residuals
    random_state=123         # int
) -> None

# ForecasterRecursiveMultiSeries, ForecasterDirectMultiVariate, ForecasterRnn
forecaster.set_out_sample_residuals(
    y_true,                  # dict[str, np.ndarray | pd.Series] (required)
    y_pred,                  # dict[str, np.ndarray | pd.Series] (required)
    append=False,            # bool
    random_state=123         # int
) -> None

# ForecasterEquivalentDate (same as single series)
forecaster.set_out_sample_residuals(
    y_true,                  # np.ndarray | pd.Series (required)
    y_pred,                  # np.ndarray | pd.Series (required)
    append=False,            # bool
    random_state=123         # int
) -> None

# NOT available in: ForecasterRecursiveClassifier, ForecasterStats
```

## Forecaster Methods: set_params() and set_lags()

```python
# set_params — available in all forecasters except ForecasterEquivalentDate
forecaster.set_params(
    params                   # dict[str, object] (required)
) -> None
# ForecasterStats also accepts dict[str, dict] for multiple models

# set_lags — available in all forecasters except ForecasterStats and ForecasterEquivalentDate
forecaster.set_lags(
    lags=None                # int | list[int] | np.ndarray | range | None
) -> None
# ForecasterDirectMultiVariate also accepts dict[str, int | list]
# ForecasterRnn: set_lags() exists but is a no-op for API consistency
```

## Method Availability Matrix

| Method | Recursive | Direct | RecursiveMultiSeries | DirectMultiVariate | Rnn | Stats | EquivalentDate | Classifier |
|--------|:---------:|:------:|:-------------------:|:-----------------:|:---:|:-----:|:--------------:|:----------:|
| `predict()` | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |
| `predict_interval()` | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | — |
| `predict_quantiles()` | ✓ | ✓ | ✓ | ✓ | — | — | — | — |
| `predict_dist()` | ✓ | ✓ | ✓ | ✓ | — | — | — | — |
| `predict_proba()` | — | — | — | — | — | — | — | ✓ |
| `set_params()` | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | — | ✓ |
| `set_lags()` | ✓ | ✓ | ✓ | ✓ | ✓* | — | — | ✓ |
| `set_out_sample_residuals()` | ✓ | ✓ | ✓ | ✓ | ✓ | — | ✓ | — |

> ✓ = supported, — = not available, ✓* = exists but is a no-op

## Cross-Validation Classes

```python
TimeSeriesFold(
    steps,                          # int (required), forecast horizon
    initial_train_size=None,        # int | str | pd.Timestamp | None
    fold_stride=None,               # int | None, if None equals steps
    window_size=None,               # int | None, set automatically by forecaster
    differentiation=None,           # int | None, set automatically by forecaster
    refit=False,                    # bool | int, refit model each fold or every n folds
    fixed_train_size=True,          # bool, fixed vs expanding window
    gap=0,                          # int, observations between train end and test start
    skip_folds=None,                # int | list[int] | None
    allow_incomplete_fold=True,     # bool
    return_all_indexes=False,       # bool
    verbose=True                    # bool
)

OneStepAheadFold(
    initial_train_size,             # int | str | pd.Timestamp (required)
    window_size=None,               # int | None, set automatically by forecaster
    differentiation=None,           # int | None, set automatically by forecaster
    return_all_indexes=False,       # bool
    verbose=True                    # bool
)
```

## Feature Selection Functions

```python
select_features(
    forecaster,              # ForecasterRecursive | ForecasterDirect
    selector,                # sklearn feature selector (RFECV, SelectFromModel, etc.)
    y,                       # pd.Series | pd.DataFrame
    exog=None,               # pd.Series | pd.DataFrame | None
    select_only=None,        # 'autoreg' | 'exog' | None (select all)
    force_inclusion=None,    # list[str] | str (regex) | None
    subsample=0.5,           # int | float, proportion or number of samples
    random_state=123,        # int
    verbose=True             # bool
) -> tuple[list[int], list[str], list[str]]
# Returns: (selected_lags, selected_window_features, selected_exog)

select_features_multiseries(
    forecaster,              # ForecasterRecursiveMultiSeries
    selector,                # sklearn feature selector
    series,                  # pd.DataFrame | dict[str, pd.Series | pd.DataFrame]
    exog=None,               # pd.Series | pd.DataFrame | dict | None
    select_only=None,        # 'autoreg' | 'exog' | None
    force_inclusion=None,    # list[str] | str (regex) | None
    subsample=0.5,           # int | float
    random_state=123,        # int
    verbose=True             # bool
) -> tuple[list[int] | dict[str, int], list[str], list[str]]
```

## Drift Detection Classes

```python
# RangeDriftDetector — lightweight out-of-range detector
RangeDriftDetector()  # No constructor parameters

RangeDriftDetector.fit(
    series=None,             # pd.DataFrame | pd.Series | dict | None
    exog=None,               # pd.DataFrame | pd.Series | dict | None
)

RangeDriftDetector.predict(
    last_window=None,        # pd.Series | pd.DataFrame | dict | None
    exog=None,               # pd.Series | pd.DataFrame | dict | None
    verbose=True,            # bool
    suppress_warnings=False  # bool
) -> tuple[bool, list[str], list[str] | dict[str, list[str]]]

# PopulationDriftDetector — statistical tests for distribution drift
PopulationDriftDetector(
    chunk_size=None,                    # int | str | None
    threshold=3,                        # int | float
    threshold_method='std',             # 'std' | 'quantile'
    max_out_of_range_proportion=0.1     # float
)

PopulationDriftDetector.fit(X)          # Reference dataset
PopulationDriftDetector.predict(X) -> tuple[pd.DataFrame, pd.DataFrame]
```

## Preprocessing Classes

```python
RollingFeatures(
    stats,                   # str | list[str], e.g. ['mean', 'std', 'min', 'max']
    window_sizes,            # int | list[int], int applies to all stats
    min_periods=None,        # int | list[int] | None
    features_names=None,     # list[str] | None, custom names for features
    fillna=None,             # str | float | None
    kwargs_stats={'ewm': {'alpha': 0.3}}  # dict | None, kwargs for specific stats
)

TimeSeriesDifferentiator(
    order=1,                 # int, differencing order
    window_size=None         # int | None
)

CalendarFeatures(
    features=None,           # list[str] | None, e.g. ['year', 'month', 'day_of_week', 'hour']
    encoding='cyclical',     # 'cyclical' | 'onehot' | None
    max_values=None          # dict[str, int] | None, max values for cyclical encoding
)
```