<!-- 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 machine learning models

This document is for skforecast v0.21.0+. If you are using an older version, check the documentation at skforecast.org.

Skforecast is a Python library that simplifies time series forecasting using machine learning. It works with any estimator compatible with the scikit-learn API (LightGBM, XGBoost, CatBoost, Keras, etc.).

## Quick Info

- Version: 0.21.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
├── stats/                   # Arima, Sarimax, Ets, Arar (sklearn-compatible wrappers)
├── preprocessing/           # TimeSeriesDifferentiator, RollingFeatures, DateTimeFeatureTransformer,
│                            # 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
- Statistical models in `stats/` are wrapped by `ForecasterStats` (in `recursive/`)
- `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) |
| 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)
```

## 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
predictions = forecaster.predict_interval(
    steps=10,
    interval=[10, 90],  # 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=[10, 90],                          # Prediction intervals as percentiles (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, best_trial = 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
)
```

## Statistical Models (ARIMA, ETS, ARAR)

Statistical models are wrapped by `ForecasterStats` for a sklearn-compatible interface. Available models: `Arima`, `Sarimax`, `Ets`, `Arar`.

```python
from skforecast.recursive import ForecasterStats
from skforecast.stats import Arima, Ets, Sarimax, Arar

# 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
forecaster = ForecasterStats(estimator=Arima(order=None, seasonal=True, 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'))
```

## 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

# 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 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 DateTimeFeatureTransformer
from skforecast.preprocessing import create_datetime_features
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

# 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 calculate_lag_autocorrelation
from skforecast.plot import backtesting_gif_creator
from skforecast.plot import set_dark_theme

# Exceptions and Warnings
from skforecast.exceptions import set_warnings_style
from skforecast.exceptions import warn_skforecast_categories
```

## 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 × 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**: Forecasters don't accept NaN by default
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.21.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.

## 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
    differentiation=None,        # e.g., 1 for first-order differencing
)

# 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 [5, 95] (90%). Here [10, 90] creates an 80% interval.
forecaster.fit(y=y_train, store_in_sample_residuals=True)
predictions_interval = forecaster.predict_interval(
    steps=10,
    interval=[10, 90],          # 80% prediction interval
    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,
)
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. 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.

## 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},
    differentiation=None,
    dropna_from_series=False, # Set True if individual series may have NaN
)

# 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
)
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

## 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=[10, 90],  # or use alpha=0.2 for 80% interval
)
```

## Auto-ARIMA (Automatic Order Selection)

```python
# Set order=None to enable automatic order selection
auto_arima = Arima(order=None, seasonal=True, 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: 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 |

## 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, best_trial = 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)
```

## 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, best_trial = 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,
)
```

## 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, best_trial = 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,
)
```

## 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, best_trial = bayesian_search_forecaster(
    ...,
    kwargs_create_study={
        'sampler': optuna.samplers.TPESampler(seed=123),
        'direction': 'minimize',
    },
    kwargs_study_optimize={
        'timeout': 600,  # seconds
        'gc_after_trial': True,
    },
)
```

## Return Values

| Function | Returns | Best trial object |
|----------|---------|:-:|
| `grid_search_*` | `pd.DataFrame` sorted by metric | — |
| `random_search_*` | `pd.DataFrame` sorted by metric | — |
| `bayesian_search_*` | `tuple[pd.DataFrame, optuna.FrozenTrial]` | ✓ |
| `*_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).

================================================================================
# SKILL: prediction-intervals
================================================================================

# Prediction Intervals

## 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.

## When to Use

Use prediction intervals to quantify forecast uncertainty. Skforecast offers three methods:

| Method | Forecasters | Description |
|--------|-------------|-------------|
| **Bootstrapping** | Recursive, Direct | Resample from training residuals |
| **Conformal** | All ML forecasters | Distribution-free intervals via conformal prediction |
| **Built-in** | ForecasterStats (ARIMA, ETS) | Parametric intervals from the statistical model |

## 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=[10, 90],              # Percentiles → 80% interval (default is [5, 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=[10, 90],
    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=[10, 90],         # Or use alpha=0.05 for 95% interval
)
```

## 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=[10, 90],
    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=[10, 90],
    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=[10, 90],
    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 [10, 90] interval
```

## Common Mistakes

1. **Forgetting `store_in_sample_residuals=True`**: Required in `fit()` before using `predict_interval(method='bootstrapping')`.
2. **Wrong default method for multi-series**: `ForecasterRecursiveMultiSeries` and `ForecasterDirectMultiVariate` default to `method='conformal'`, not `'bootstrapping'`.
3. **Mixing `alpha` and `interval`**: `ForecasterStats` supports both `alpha` (e.g., `alpha=0.05` for 95% interval) and `interval=[lo, hi]`. ML forecasters only support `interval`.
4. **Not evaluating coverage**: Always check if actual coverage matches nominal interval width.

---
### Reference: interval-compatibility

# Prediction Intervals — Compatibility Reference

## 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=[5, 95],                # percentiles → 90% interval
    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 — 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=[5, 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=[5, 95],
    use_in_sample_residuals=True,
    n_boot=None,                     # Any, accepted but ignored
    use_binned_residuals=None,       # Any, accepted but ignored
    random_state=None,               # Any, accepted but ignored
)
# NOTE: `n_boot`, `use_binned_residuals`, 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 percentiles |
| `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` |
| `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: feature-engineering
================================================================================

# Feature Engineering

## References

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, rolling statistics, cyclical encoding,
sunlight features, differencing, or data scaling.

## Overview

| Tool | Package | Purpose |
|------|---------|---------|
| `DatetimeFeatures` | feature_engine | Extract calendar features from datetime index |
| `CyclicalFeatures` | feature_engine | Encode cyclical features with sin/cos |
| `RollingFeatures` | skforecast | Rolling window statistics (mean, std, min, max, etc.) |
| `differentiation` param | skforecast | Make non-stationary series stationary |
| `astral` | astral | Sunrise, sunset, daylight hours |

## Calendar Features with feature_engine

### Manual extraction (pandas)

```python
import pandas as pd

# Data must have a DatetimeIndex with frequency set
data = data.asfreq('h')

data['year'] = data.index.year
data['month'] = data.index.month
data['day_of_week'] = data.index.dayofweek
data['hour'] = data.index.hour
```

### Automated extraction (DatetimeFeatures)

```python
from feature_engine.datetime import DatetimeFeatures

features_to_extract = ['month', 'week', 'day_of_week', 'hour']
calendar_transformer = DatetimeFeatures(
    variables           = 'index',
    features_to_extract = features_to_extract,
    drop_original       = True,
)

calendar_features = calendar_transformer.fit_transform(data)
```

> `DatetimeFeatures` is sklearn-compatible and can be passed directly as
> `transformer_exog` in skforecast forecasters.

## Cyclical Encoding

Cyclical features (hour, day_of_week, month) should NOT be treated as linear
integers — hour 23 is only 1 hour from hour 0. Use sin/cos encoding to
preserve the cyclical relationship.

```python
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures

# Step 1: Extract calendar features
features_to_extract = ['month', 'week', 'day_of_week', 'hour']
calendar_transformer = DatetimeFeatures(
    variables           = 'index',
    features_to_extract = features_to_extract,
    drop_original       = True,
)
calendar_features = calendar_transformer.fit_transform(data)

# Step 2: Encode as cyclical (sin/cos)
features_to_encode = ['month', 'week', 'day_of_week', 'hour']
max_values = {
    'month': 12,
    'week': 52,
    'day_of_week': 7,
    'hour': 24,
}
cyclical_encoder = CyclicalFeatures(
    variables     = features_to_encode,
    max_values    = max_values,
    drop_original = True,
)
exog_calendar = cyclical_encoder.fit_transform(calendar_features)
# Produces columns: month_sin, month_cos, week_sin, week_cos, ...
```

## Sunlight Features

Sunrise/sunset times can be powerful features for energy, transport, or
activity-related series.

```python
from astral.sun import sun
from astral import LocationInfo

location = LocationInfo('Washington, D.C.', 'USA')
sunrise_hour = [sun(location.observer, date=date)['sunrise'] for date in data.index]
sunset_hour = [sun(location.observer, date=date)['sunset'] for date in data.index]

# Round to the nearest hour
sunrise_hour = pd.Series(sunrise_hour, index=data.index).dt.round('h').dt.hour
sunset_hour = pd.Series(sunset_hour, index=data.index).dt.round('h').dt.hour

sun_light_features = pd.DataFrame({
    'sunrise_hour': sunrise_hour,
    'sunset_hour': sunset_hour,
})
sun_light_features['daylight_hours'] = (
    sun_light_features['sunset_hour'] - sun_light_features['sunrise_hour']
)
```

## 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}},
)
```

## 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
```

## Data Transformers (Scaling)

```python
from sklearn.preprocessing import StandardScaler, MinMaxScaler

# Scale target variable — transformer applied automatically during fit/predict
forecaster = ForecasterRecursive(
    estimator=LGBMRegressor(),
    lags=24,
    transformer_y=StandardScaler(),
    transformer_exog=StandardScaler(),
)

# For multi-series, different transformers per series
from skforecast.recursive import ForecasterRecursiveMultiSeries

forecaster = ForecasterRecursiveMultiSeries(
    estimator=LGBMRegressor(),
    lags=24,
    transformer_series={
        'series_1': StandardScaler(),
        'series_2': MinMaxScaler(),
    },
)
```

## Combining Features — Full Example

```python
import pandas as pd
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from skforecast.preprocessing import RollingFeatures
from skforecast.recursive import ForecasterRecursive
from sklearn.preprocessing import StandardScaler
from lightgbm import LGBMRegressor

# 1. Calendar features with cyclical encoding
calendar_transformer = DatetimeFeatures(
    variables='index',
    features_to_extract=['month', 'day_of_week', 'hour'],
    drop_original=True,
)
cyclical_encoder = CyclicalFeatures(
    variables=['month', 'day_of_week', 'hour'],
    max_values={'month': 12, 'day_of_week': 7, 'hour': 24},
    drop_original=True,
)
exog_calendar = cyclical_encoder.fit_transform(
    calendar_transformer.fit_transform(data)
)

# 2. Combine with other exogenous variables
exog = pd.concat([exog_external, exog_calendar], axis=1)

# 3. 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])
```

## 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). Always use sin/cos encoding.
2. **Forgetting frequency on index**: Calendar transformers require `DatetimeIndex` with frequency set (`data.asfreq('h')`).
3. **Not covering forecast horizon with exog**: Calendar features for `predict()` must include future dates covering the entire forecast horizon.
4. **Over-engineering features**: Start with lags only, then add rolling features and calendar features incrementally. Validate each addition with backtesting.

---
### 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 |
|------------|:--:|
| ForecasterRecursive | ✓ |
| ForecasterDirect | ✓ |
| ForecasterRecursiveMultiSeries | ✓ |
| ForecasterDirectMultiVariate | ✓ |
| ForecasterRecursiveClassifier | ✓ |
| ForecasterRnn | — |
| ForecasterStats | — |
| ForecasterEquivalentDate | — |

## 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

## 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 |

## 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)

## 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=[10, 90],
    use_in_sample_residuals=True,
)
```

## 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=[10, 90],
    use_in_sample_residuals=True,
)
# NOTE: 'bootstrapping' is NOT supported for ForecasterRnn
# NOTE: predict_quantiles() and predict_dist() are NOT available
```

================================================================================
# 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.

## 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 |
| **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 |

## 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
```

## 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 `DateTimeFeatureTransformer` or `create_datetime_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

## Deprecated Import Paths

The most frequent LLM error. Old import paths no longer exist.

| Wrong (Deprecated) | Correct (v0.21.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(...))` |

## 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: handle 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, ...)  # Error!

# ✅ CORRECT: use backtesting_stats for statistical models
from skforecast.model_selection import backtesting_stats
backtesting_stats(forecaster=forecaster_stats, ...)

# ❌ WRONG: using backtesting_forecaster with ForecasterRecursiveMultiSeries
backtesting_forecaster(forecaster=forecaster_multi, ...)  # Error!

# ✅ CORRECT: use backtesting_forecaster_multiseries
from skforecast.model_selection import backtesting_forecaster_multiseries
backtesting_forecaster_multiseries(forecaster=forecaster_multi, series=series, ...)
```

## Wrong Search Function

```python
# ❌ WRONG: grid_search_forecaster with ForecasterStats
grid_search_forecaster(forecaster=forecaster_stats, ...)

# ✅ CORRECT: grid_search_stats for statistical models
from skforecast.model_selection import grid_search_stats
grid_search_stats(forecaster=forecaster_stats, ...)

# ❌ WRONG: grid_search_forecaster with ForecasterRecursiveMultiSeries
grid_search_forecaster(forecaster=forecaster_multi, ...)

# ✅ CORRECT: grid_search_forecaster_multiseries
from skforecast.model_selection import grid_search_forecaster_multiseries
grid_search_forecaster_multiseries(forecaster=forecaster_multi, series=series, ...)
```

## 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 |

================================================================================
# 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)

### 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
- `DateTimeFeatureTransformer` — 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
    weight_func=None,         # Callable to weight training samples by index position
    differentiation=None,     # int, differencing order applied before training
    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
    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
    weight_func=None,         # Callable to weight training samples by index position
    differentiation=None,     # int, differencing order applied before training
    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
    weight_func=None,              # Callable to weight training samples by index position
    differentiation=None,          # int, differencing order applied before training
    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
    weight_func=None,         # Callable to weight training samples by index position
    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(
    levels,                            # str | list[str] (required), target series names
    lags,                              # int | list[int] | np.ndarray | range (required)
    estimator=None,                    # Keras model (use create_and_compile_model)
    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()
    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).
```

## 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
```

## 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=[5, 95],                   # float | list[float] | tuple[float]
    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=[5, 95],                   # float | 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_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=[5, 95],                   # float | 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 (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=[5, 95],                   # float | 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

# 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
    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=[5, 95],                   # float | list[float] | tuple[float]
    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=[5, 95],                   # float | list[float] | tuple[float]
    use_in_sample_residuals=True,       # bool
    suppress_warnings=False,            # bool
    n_boot=None,                        # Any, ignored (API compatibility)
    use_binned_residuals=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
)

DateTimeFeatureTransformer(
    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
)
```