Hyperparameter tuning and lags selection¶
Hyperparameter tuning is a key step in building accurate and robust machine learning models. Hyperparameters are configuration values that cannot be learned directly from data and must be defined by the user before training. These values can significantly affect model performance, and carefully tuning them helps improve both accuracy and generalization.
In forecasting models, the selection of lags (past time steps used as predictors) is considered an additional hyperparameter, as it directly influences the model's input structure and learning capacity.
Hyperparameter tuning consists of systematically evaluating combinations of hyperparameters (including lags) to find the configuration that yields the best predictive performance. The skforecast library supports several tuning strategies: grid search, random search, and Bayesian search. These strategies can be used with either backtesting or one-step-ahead validation to determine the optimal parameter set for a given forecasting task.
✏️ Note
All backtesting and hyperparameter search functions in the model_selection module include the n_jobs argument, enabling multi-process parallelization to improve computational performance.
Its effectiveness depends on factors like the estimator type, the number of model fits to perform, and the volume of data. When n_jobs is set to 'auto', the level of parallelization is automatically determined using heuristic rules designed to select the most efficient configuration for each scenario.
For more information, see the guide Parallelization in skforecast.
Validation strategies¶
Hyperparameter and lag tuning involves systematically testing different values or combinations of hyperparameters (and/or lags) to find the optimal configuration that gives the best performance. The skforecast library provides two different methods to evaluate each candidate configuration:
Backtesting: Simulates a real deployment scenario by generating multi-step forecasts in repeated iterations, using the defined forecast horizon and retraining frequency. This approach provides a realistic estimate of performance over time. Use the
TimeSeriesFoldclass for this validation strategy. More information.One-Step-Ahead: Evaluates model performance using only one-step-ahead forecasts (). This method is faster because the model is trained only once and each prediction uses the observed values of the lags (no recursive forecasting is needed), but it only tests the model's performance in the immediate next time step. Use the
OneStepAheadFoldclass for the one-step-ahead strategy. More information.
The two methods often, but not always, select similar configurations. A common workflow is to use one-step-ahead validation to explore a large search space quickly, and then backtest the final model (or the best few candidates) to obtain a reliable estimate of its multi-step performance.
Libraries and data¶
# Libraries
# ==============================================================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from lightgbm import LGBMRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from skforecast.datasets import fetch_dataset
from skforecast.recursive import ForecasterRecursive
from skforecast.plot import set_dark_theme
from skforecast.model_selection import (
TimeSeriesFold,
OneStepAheadFold,
backtesting_forecaster,
grid_search_forecaster,
random_search_forecaster,
bayesian_search_forecaster
)
The data used in this guide is the h2o dataset: the monthly expenditure, in Australian dollars, on corticosteroid drugs in the Australian health system between July 1991 and June 2008 (204 observations).
# Download data
# ==============================================================================
data = fetch_dataset(
name="h2o", raw=True, kwargs_read={"names": ["y", "datetime"], "header": 0}
)
# Data preprocessing
# ==============================================================================
data['datetime'] = pd.to_datetime(data['datetime'], format='%Y-%m-%d')
data = data.set_index('datetime')
data = data.sort_index()
data = data.asfreq('MS')
data = data[['y']]
data.head(3)
╭────────────────────────────────────── h2o ───────────────────────────────────────╮ │ Description: │ │ Monthly expenditure ($AUD) on corticosteroid drugs that the Australian health │ │ system had between 1991 and 2008. │ │ │ │ Source: │ │ Hyndman R (2023). fpp3: Data for Forecasting: Principles and Practice(3rd │ │ Edition). http://pkg.robjhyndman.com/fpp3package/,https://github.com/robjhyndman │ │ /fpp3package, http://OTexts.com/fpp3. │ │ │ │ URL: │ │ https://raw.githubusercontent.com/skforecast/skforecast- │ │ datasets/main/data/h2o.csv │ │ │ │ Shape: 204 rows x 2 columns │ ╰──────────────────────────────────────────────────────────────────────────────────╯
| y | |
|---|---|
| datetime | |
| 1991-07-01 | 0.429795 |
| 1991-08-01 | 0.400906 |
| 1991-09-01 | 0.432159 |
The series is split chronologically into three consecutive sets. A random split must never be used with time series, since it would let the model learn from observations that come after the ones it predicts.
Training (until
end_train): used to fit the models evaluated during the search.Validation (from
end_traintoend_val): used to compare the candidate configurations.Test (after
end_val): held out during the search and used only at the end, to estimate the performance of the selected model.
⚠ Warning
The search functions receive only the training and validation data (data.loc[:end_val]), and initial_train_size is set to the end of the training set. This way, each candidate configuration is trained with the training set and evaluated on the validation set.
The test set must not be passed to the search. If it were, the configuration would be selected using the same data employed to evaluate it, and the reported error would be optimistic (data leakage).
# Train-val-test dates
# ==============================================================================
# The time 23:59:00 keeps 2001-01-01 in the training set, so that
# data.loc[end_train:] starts in the next month and the sets do not overlap
end_train = '2001-01-01 23:59:00'
end_val = '2006-01-01 23:59:00'
print(
f"Train dates : {data.index.min()} --- {data.loc[:end_train].index.max()}"
f" (n={len(data.loc[:end_train])})"
)
print(
f"Validation dates : {data.loc[end_train:].index.min()} --- {data.loc[:end_val].index.max()}"
f" (n={len(data.loc[end_train:end_val])})"
)
print(
f"Test dates : {data.loc[end_val:].index.min()} --- {data.index.max()}"
f" (n={len(data.loc[end_val:])})"
)
print()
# Plot
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(7, 3))
data.loc[:end_train, 'y'].plot(ax=ax, label='train')
data.loc[end_train:end_val, 'y'].plot(ax=ax, label='validation')
data.loc[end_val:, 'y'].plot(ax=ax, label='test')
ax.legend()
plt.show()
Train dates : 1991-07-01 00:00:00 --- 2001-01-01 00:00:00 (n=115) Validation dates : 2001-02-01 00:00:00 --- 2006-01-01 00:00:00 (n=60) Test dates : 2006-02-01 00:00:00 --- 2008-06-01 00:00:00 (n=29)
Grid search¶
Grid search is a popular hyperparameter tuning technique that evaluates an exhaustive list of combinations of hyperparameters and lags to find the optimal configuration for a forecasting model. To perform a grid search with the skforecast library, two grids are needed: one with different lags (lags_grid) and another with the hyperparameters (param_grid). The lags_grid can be a list of lag configurations (an int means all lags from 1 to that value, a list means specific lags) or a dict whose keys are custom labels for each configuration. These labels are shown in the lags_label column of the results.
The grid search process involves the following steps:
grid_search_forecasterreplaces thelagsargument with the first option appearing inlags_grid.The function validates all combinations of hyperparameters presented in
param_gridusing backtesting or one-step-ahead validation.The function repeats these two steps until it has evaluated all possible combinations of lags and hyperparameters.
If
return_best = True, the original forecaster is trained with the best lags and hyperparameters configuration found during the grid search process, using all the data passed to the search.
In the examples of this guide, the lags and hyperparameters of a ForecasterRecursive with a LightGBM estimator are tuned.
💡 Tip
When using backtesting as the validation strategy, the computational cost of the tuning largely depends on the strategy used to evaluate each hyperparameter combination. In general, the more re-trainings required, the longer the tuning process will take.
To speed up the prototyping phase, a two-step approach is recommended. First, run the search with refit=False to explore a broad range of values quickly. Then, refine the search within the most promising region using a tailored backtesting strategy aligned with the specific needs of the use case.
For more guidance, refer to the following resource: Which backtesting strategy should I use?.
# Grid search hyperparameters and lags
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Lags used as predictors
lags_grid = {
'lags_1': 3,
'lags_2': 10,
'lags_3': [1, 2, 3, 20]
}
# Estimator hyperparameters
param_grid = {
'n_estimators': [50, 100],
'max_depth': [5, 10, 15]
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = end_train, # A date can be used instead of a number of observations
refit = False
)
results = grid_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
param_grid = param_grid,
lags_grid = lags_grid,
cv = cv,
metric = 'mean_squared_error',
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
results
| lags | lags_label | params | mean_squared_error | max_depth | n_estimators | |
|---|---|---|---|---|---|---|
| 0 | [1, 2, 3] | lags_1 | {'max_depth': 5, 'n_estimators': 100} | 0.043875 | 5 | 100 |
| 1 | [1, 2, 3] | lags_1 | {'max_depth': 10, 'n_estimators': 100} | 0.043875 | 10 | 100 |
| 2 | [1, 2, 3] | lags_1 | {'max_depth': 15, 'n_estimators': 100} | 0.043875 | 15 | 100 |
| 3 | [1, 2, 3, 20] | lags_3 | {'max_depth': 15, 'n_estimators': 100} | 0.044074 | 15 | 100 |
| 4 | [1, 2, 3, 20] | lags_3 | {'max_depth': 10, 'n_estimators': 100} | 0.044074 | 10 | 100 |
| 5 | [1, 2, 3, 20] | lags_3 | {'max_depth': 5, 'n_estimators': 100} | 0.044074 | 5 | 100 |
| 6 | [1, 2, 3] | lags_1 | {'max_depth': 5, 'n_estimators': 50} | 0.045423 | 5 | 50 |
| 7 | [1, 2, 3] | lags_1 | {'max_depth': 15, 'n_estimators': 50} | 0.045423 | 15 | 50 |
| 8 | [1, 2, 3] | lags_1 | {'max_depth': 10, 'n_estimators': 50} | 0.045423 | 10 | 50 |
| 9 | [1, 2, 3, 20] | lags_3 | {'max_depth': 15, 'n_estimators': 50} | 0.046221 | 15 | 50 |
| 10 | [1, 2, 3, 20] | lags_3 | {'max_depth': 5, 'n_estimators': 50} | 0.046221 | 5 | 50 |
| 11 | [1, 2, 3, 20] | lags_3 | {'max_depth': 10, 'n_estimators': 50} | 0.046221 | 10 | 50 |
| 12 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 5, 'n_estimators': 100} | 0.047896 | 5 | 100 |
| 13 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 10, 'n_estimators': 100} | 0.047896 | 10 | 100 |
| 14 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 15, 'n_estimators': 100} | 0.047896 | 15 | 100 |
| 15 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 15, 'n_estimators': 50} | 0.051399 | 15 | 50 |
| 16 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 5, 'n_estimators': 50} | 0.051399 | 5 | 50 |
| 17 | [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] | lags_2 | {'max_depth': 10, 'n_estimators': 50} | 0.051399 | 10 | 50 |
Since return_best = True, the forecaster object is updated with the best configuration found and refitted using all the data passed to the search (in this example, the training and validation sets, data.loc[:end_val]). The test set remains unseen, so it can be used to obtain an unbiased estimate of the performance of the final model. This final model can then be used for future predictions on new data.
The best configuration found uses lags [1, 2, 3] and n_estimators = 100, with a validation mean squared error of 0.043875. The three values of max_depth obtain the same error (see the note below), and the forecaster keeps the first row of the results (max_depth = 5), as shown in its parameters.
forecaster
ForecasterRecursive
General Information
- Estimator: LGBMRegressor
- Lags: [1 2 3]
- Window features: None
- Calendar features: None
- Window size: 3
- Series name: y
- Exogenous included: False
- Categorical features: auto
- Weight function included: False
- Differentiation order: None
- Drop NaN from series: False
- Creation date: 2026-10-08 19:49:23
- Last fit date: 2026-10-08 19:49:24
- Skforecast version: 0.26.0
- Python version: 3.14.3
- Forecaster id: None
Exogenous Variables
None
Data Transformations
- Transformer for y: None
- Transformer for exog: None
Training Information
- Training range: [Timestamp('1991-07-01 00:00:00'), Timestamp('2006-01-01 00:00:00')]
- Training index type: DatetimeIndex
- Training index frequency: MS
Estimator Parameters
-
{'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0, 'importance_type': 'split', 'learning_rate': 0.1, 'max_depth': 5, 'min_child_samples': 20, 'min_child_weight': 0.001, 'min_split_gain': 0.0, 'n_estimators': 100, 'n_jobs': None, 'num_leaves': 31, 'objective': None, 'random_state': 123, 'reg_alpha': 0.0, 'reg_lambda': 0.0, 'subsample': 1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0, 'verbose': -1}
Fit Kwargs
-
{}
The search only used the training and validation sets. To estimate the performance of the selected model on unseen data, the refitted forecaster is evaluated on the test set with backtesting_forecaster. With initial_train_size equal to the number of observations up to the end of the validation set and refit=False, the forecaster is trained with the same data used by return_best and then predicts the test set in folds of 12 months.
# Backtesting of the best model on the test set
# ==============================================================================
cv_test = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_val]),
refit = False
)
metric_test, predictions_test = backtesting_forecaster(
forecaster = forecaster,
y = data['y'],
cv = cv_test,
metric = 'mean_squared_error',
show_progress = False
)
metric_test
| mean_squared_error | |
|---|---|
| 0 | 0.058782 |
Random search¶
Random search (random_search_forecaster) is another hyperparameter tuning strategy available in the skforecast library. In contrast to grid search, which tries out all possible combinations of hyperparameters and lags, random search evaluates only a fixed number of combinations of hyperparameter values, sampled at random from the specified possibilities. The number of combinations that are evaluated is given by n_iter.
It is important to note that random sampling is only applied to the model hyperparameters, but not to the lags. All lags specified by the user are evaluated. In the following example, n_iter = 5 combinations of hyperparameters are evaluated for each of the 2 lag configurations (the same 5 combinations for both), 10 models in total.
# Random search hyperparameters and lags
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Lags used as predictors
lags_grid = [3, 5]
# Estimator hyperparameters
param_distributions = {
'n_estimators': np.arange(start=10, stop=100, step=1, dtype=int),
'max_depth': np.arange(start=5, stop=30, step=1, dtype=int)
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]),
refit = False,
)
results = random_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
lags_grid = lags_grid,
param_distributions = param_distributions,
cv = cv,
n_iter = 5,
metric = 'mean_squared_error',
return_best = True,
random_state = 123,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
results.head(4)
| lags | lags_label | params | mean_squared_error | n_estimators | max_depth | |
|---|---|---|---|---|---|---|
| 0 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'n_estimators': 96, 'max_depth': 19} | 0.043131 | 96 | 19 |
| 1 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'n_estimators': 94, 'max_depth': 28} | 0.043171 | 94 | 28 |
| 2 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'n_estimators': 77, 'max_depth': 17} | 0.043663 | 77 | 17 |
| 3 | [1, 2, 3] | [1, 2, 3] | {'n_estimators': 96, 'max_depth': 19} | 0.043868 | 96 | 19 |
The best configuration uses lags 1 to 5, n_estimators = 96 and max_depth = 19, with a validation mean squared error of 0.043131. It is slightly lower than the best error of the grid search (0.043875), and it was obtained with 10 evaluations instead of 18. The improvement comes from the lag configuration (5 lags were not included in the first grid): with 3 lags, the error (0.043868) is practically the same as in the grid search.
Rows 0 and 3 have the same hyperparameters, which shows that the sampled combinations are reused for every lag configuration.
Bayesian search¶
Grid and random search can yield good results, especially when the search space is well-defined. However, these methods do not consider past results, which limits their ability to focus on the most promising regions and avoid uninformative ones.
A more efficient alternative is Bayesian optimization (bayesian_search_forecaster), which builds a probabilistic model of the objective function, typically the validation metric (e.g. MAE or MSE). Based on the results observed so far, the algorithm iteratively refines the search, concentrating on regions with the highest potential. This approach reduces the number of evaluations needed by prioritizing the most relevant hyperparameter combinations. It is especially useful when the search space is large or model training is computationally expensive.
✏️ Note
Unlike the grid and random search functions, bayesian_search_forecaster has no lags_grid argument. Instead, the lags are included directly in the search_space, so that they are optimized jointly with the other estimator hyperparameters during the search.
In skforecast, Bayesian optimization is implemented using Optuna (by default, with the Tree-structured Parzen Estimator, TPE, sampler) and its Study object. The goal of the optimization is to minimize the metric returned by the validation strategy (either backtesting or one-step-ahead).
You can customize the optimization process by passing additional arguments through the kwargs_create_study and kwargs_study_optimize parameters. These are forwarded to Optuna’s create_study and optimize method, respectively.
To define the hyperparameter search space, the search_space argument must be a function that takes an Optuna Trial object and returns a dictionary of parameters to evaluate.
✏️ Note
The direction of the hyperparameter optimization can be changed by setting the direction argument in kwargs_create_study to 'minimize' or 'maximize'. By default, the direction is set to 'minimize' in regression and 'maximize' in classification tasks.
kwargs_create_study = {
'direction': 'minimize' # or 'maximize'
}
# Bayesian search hyperparameters and lags with Optuna
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Search space
def search_space(trial):
return {
'lags' : trial.suggest_categorical('lags', [3, 5]),
'n_estimators' : trial.suggest_int('n_estimators', 10, 20),
'min_child_samples': trial.suggest_int('min_child_samples', 1, 10),
'colsample_bytree' : trial.suggest_float('colsample_bytree', 0.5, 1.0)
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]), # A date can also be used: end_train
refit = False,
)
results, study = bayesian_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
search_space = search_space,
cv = cv,
metric = 'mean_absolute_error',
n_trials = 20,
random_state = 123,
return_best = False,
n_jobs = 'auto',
verbose = False,
show_progress = True,
suppress_warnings = False,
kwargs_create_study = None,
kwargs_study_optimize = None
)
results.head(4)
| trial_number | lags | params | mean_absolute_error | n_estimators | min_child_samples | colsample_bytree | |
|---|---|---|---|---|---|---|---|
| 0 | 17 | [1, 2, 3, 4, 5] | {'n_estimators': 14, 'min_child_samples': 1, '... | 0.112030 | 14.0 | 1.0 | 0.963702 |
| 1 | 2 | [1, 2, 3, 4, 5] | {'n_estimators': 14, 'min_child_samples': 1, '... | 0.142633 | 14.0 | 1.0 | 0.699022 |
| 2 | 15 | [1, 2, 3, 4, 5] | {'n_estimators': 13, 'min_child_samples': 1, '... | 0.146512 | 13.0 | 1.0 | 0.539784 |
| 3 | 10 | [1, 2, 3, 4, 5] | {'n_estimators': 12, 'min_child_samples': 1, '... | 0.151249 | 12.0 | 1.0 | 0.562491 |
The best trial (number 2) uses lags 1 to 5, n_estimators = 14 and min_child_samples = 1, with a mean absolute error of 0.142633. Note that this search uses the mean absolute error (MAE), whereas the previous ones used the mean squared error (MSE), so their values cannot be compared.
The search space includes min_child_samples because, as explained in the grid search section, it is the hyperparameter that limits the size of the trees in such a small dataset. The four best trials use its smallest possible value (1), which allows the trees to grow larger.
In addition to the results DataFrame, bayesian_search_forecaster also returns the Optuna study object. This gives the user full access to the optimization internals, enabling further analysis, visualization, or customization of the search process.
For more information, refer to the Optuna Study class.
# Trials history of the Optuna study
# ==============================================================================
study.trials_dataframe().head(4)
| number | value | datetime_start | datetime_complete | duration | params_colsample_bytree | params_lags | params_min_child_samples | params_n_estimators | user_attrs_mean_absolute_error | system_attrs_tpe:relative_params:0 | state | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0.176708 | 2026-10-08 19:49:25.640627 | 2026-10-08 19:49:25.664989 | 0 days 00:00:00.024362 | 0.859734 | 3 | 6 | 12 | 0.176708 | NaN | COMPLETE |
| 1 | 1 | 0.184199 | 2026-10-08 19:49:25.665252 | 2026-10-08 19:49:25.699322 | 0 days 00:00:00.034070 | 0.696059 | 5 | 5 | 17 | 0.184199 | NaN | COMPLETE |
| 2 | 2 | 0.142633 | 2026-10-08 19:49:25.699554 | 2026-10-08 19:49:25.747555 | 0 days 00:00:00.048001 | 0.699022 | 5 | 1 | 14 | 0.142633 | NaN | COMPLETE |
| 3 | 3 | 0.199588 | 2026-10-08 19:49:25.747879 | 2026-10-08 19:49:25.768530 | 0 days 00:00:00.020651 | 0.765914 | 3 | 6 | 11 | 0.199588 | NaN | COMPLETE |
# Optuna best trial in the study
# ==============================================================================
study.best_trial
FrozenTrial(number=17, state=<TrialState.COMPLETE: 1>, values=[0.1120303641092155], datetime_start=datetime.datetime(2026, 10, 8, 19, 49, 26, 196978), datetime_complete=datetime.datetime(2026, 10, 8, 19, 49, 26, 247308), params={'lags': 5, 'n_estimators': 14, 'min_child_samples': 1, 'colsample_bytree': 0.9637021314088274}, user_attrs={'mean_absolute_error': 0.1120303641092155}, system_attrs={'tpe:relative_params:0': '{"colsample_bytree": 0.9637021314088274, "lags": 5, "min_child_samples": 1, "n_estimators": 14}'}, intermediate_values={}, distributions={'lags': CategoricalDistribution(choices=(3, 5)), 'n_estimators': IntDistribution(high=20, log=False, low=10, step=1), 'min_child_samples': IntDistribution(high=10, log=False, low=1, step=1), 'colsample_bytree': FloatDistribution(high=1.0, log=False, low=0.5, step=None)}, trial_id=17, value=None)
One-step-ahead validation¶
As described in Validation strategies, the one-step-ahead strategy evaluates each configuration with one-step-ahead forecasts () only. It is faster than backtesting but does not measure multi-step performance. Use the OneStepAheadFold class for the one-step-ahead strategy.
With this strategy, the forecaster is trained only once, using the first initial_train_size observations. Then, each observation of the validation set is predicted using the actual (observed) values of its lags. Since no recursive predictions are needed, errors do not accumulate along the forecast horizon.
💡 Tip
For a more detailed comparison of the results (execution time and metric) obtained with each strategy, visit Hyperparameters and lags search: backtesting vs one-step-ahead.
# Bayesian search with OneStepAheadFold
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Search space: the search_space function defined in the previous section is reused
# Folds
cv = OneStepAheadFold(
initial_train_size = len(data.loc[:end_train]) # A date can also be used: end_train
)
results, study = bayesian_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
search_space = search_space,
cv = cv,
metric = 'mean_absolute_error',
n_trials = 20, # Increase the number of trials to get better results
random_state = 123,
return_best = False,
n_jobs = 'auto',
verbose = False,
show_progress = True,
suppress_warnings = False,
kwargs_create_study = None,
kwargs_study_optimize = None
)
results.head(4)
╭─────────────────────────── OneStepAheadValidationWarning ────────────────────────────╮ │ One-step-ahead predictions are used for faster model comparison, but they may not │ │ fully represent multi-step prediction performance. It is recommended to backtest the │ │ final model for a more accurate multi-step performance estimate. │ │ │ │ Category : skforecast.exceptions.OneStepAheadValidationWarning │ │ Location : │ │ /Users/javier.escobar/code/github/skforecast/skforecast/model_selection/_utils.py:65 │ │ 3 │ │ Suppress : warnings.simplefilter('ignore', category=OneStepAheadValidationWarning) │ ╰──────────────────────────────────────────────────────────────────────────────────────╯
| trial_number | lags | params | mean_absolute_error | n_estimators | min_child_samples | colsample_bytree | |
|---|---|---|---|---|---|---|---|
| 0 | 16 | [1, 2, 3, 4, 5] | {'n_estimators': 19, 'min_child_samples': 7, '... | 0.178049 | 19.0 | 7.0 | 0.960016 |
| 1 | 9 | [1, 2, 3, 4, 5] | {'n_estimators': 20, 'min_child_samples': 6, '... | 0.181041 | 20.0 | 6.0 | 0.806447 |
| 2 | 15 | [1, 2, 3, 4, 5] | {'n_estimators': 19, 'min_child_samples': 6, '... | 0.181092 | 19.0 | 6.0 | 0.879257 |
| 3 | 7 | [1, 2, 3, 4, 5] | {'n_estimators': 19, 'min_child_samples': 10, ... | 0.181476 | 19.0 | 10.0 | 0.750918 |
The best configurations found with the two validation strategies differ: backtesting selected n_estimators = 14 and min_child_samples = 1, while one-step-ahead validation selected n_estimators = 20 and min_child_samples = 6 (both with 5 lags). The three best one-step-ahead trials obtain exactly the same error because they only differ in colsample_bytree, and with 5 lags all of these values (0.831, 0.896 and 0.806) make LightGBM sample the same 4 of the 5 predictors, so the fitted models are identical. The metric values (0.143 and 0.181) should not be compared directly, since they measure the error of different forecasts: 12-step-ahead forecasts in the first case and one-step-ahead forecasts in the second. This is why the model selected with one-step-ahead validation should be evaluated with backtesting before being used.
💡 Tip
With one-step-ahead validation, the selected value of n_estimators (20) is the upper limit of its range (10 to 20). When the best value lies on the boundary of the search space, the optimum may be outside the explored range, so it is worth widening the range (for example, n_estimators between 10 and 200) and repeating the search. In the backtesting search, min_child_samples = 1 is also on the boundary, but in this case 1 is already the smallest meaningful value.
Hyperparameter tuning with custom metric¶
In addition to standard metrics such as mean_squared_error or mean_absolute_error, users can define custom metric functions, provided they accept the arguments y_true (true values), y_pred (predicted values) and optionally y_train (train values), and return a numeric value (float or int). The arguments y_true and y_pred are pandas Series indexed by the dates of the predictions, so the index can be used inside the function to select the observations to evaluate.
This flexibility allows evaluating model performance under specific conditions, for example, focusing only on certain months, days, non-holiday periods, or the last step of the forecast horizon.
To illustrate this, consider a scenario where a 12-month forecast is generated, but only the last three months of each year (October, November and December) are relevant for evaluation. This can be handled by defining a custom metric function that filters the desired months before computing the error, and then passing that function to the backtesting or hyperparameter tuning process.
The example below shows how to optimize model parameters using a custom metric focused on October, November and December.
💡 Tip
More information about time series forecasting metrics can be found in the Metrics guide.
# Custom metric
# ==============================================================================
def custom_metric(y_true, y_pred, y_train=None):
"""
Calculate the mean squared error using only the predicted values of the last
3 months of the year.
"""
mask = y_true.index.month.isin([10, 11, 12])
metric = mean_squared_error(y_true[mask], y_pred[mask])
return metric
# Grid search hyperparameters and lags with custom metric
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Lags used as predictors
lags_grid = [3, 10, [1, 2, 3, 20]]
# Estimator hyperparameters
param_grid = {
'n_estimators': [50, 100],
'max_depth': [5, 10, 15]
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]),
refit = False,
)
results = grid_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
cv = cv,
param_grid = param_grid,
lags_grid = lags_grid,
metric = custom_metric,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
results.head(4)
| lags | lags_label | params | custom_metric | max_depth | n_estimators | |
|---|---|---|---|---|---|---|
| 0 | [1, 2, 3, 20] | [1, 2, 3, 20] | {'max_depth': 15, 'n_estimators': 100} | 0.068182 | 15 | 100 |
| 1 | [1, 2, 3, 20] | [1, 2, 3, 20] | {'max_depth': 10, 'n_estimators': 100} | 0.068182 | 10 | 100 |
| 2 | [1, 2, 3, 20] | [1, 2, 3, 20] | {'max_depth': 5, 'n_estimators': 100} | 0.068182 | 5 | 100 |
| 3 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 5, 'n_estimators': 100} | 0.070472 | 5 | 100 |
With the custom metric, the best configuration uses lags [1, 2, 3, 20] (custom metric = 0.068182), whereas the grid search based on the mean squared error selected lags [1, 2, 3]. Adding lag 20 slightly worsens the overall error (0.044074 vs 0.043875 in the first grid search) but improves the forecasts of the last three months of the year, which are the ones that matter in this scenario. The metric used in the search should reflect the goal of the use case.
Compare multiple metrics¶
The functions grid_search_forecaster, random_search_forecaster, and bayesian_search_forecaster support the evaluation of multiple metrics for each forecaster configuration by passing a list of metric functions. This list can combine metrics passed by name (e.g. 'mean_absolute_error'), metric functions (e.g. mean_squared_error from scikit-learn) and custom-defined ones.
When multiple metrics are provided, the first metric in the list is used to select the best model.
💡 Tip
More information about time series forecasting metrics can be found in the Metrics guide.
# Grid search hyperparameters and lags with multiple metrics
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Lags used as predictors
lags_grid = [3, 10, [1, 2, 3, 20]]
# Estimator hyperparameters
param_grid = {
'n_estimators': [50, 100],
'max_depth': [5, 10, 15]
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]),
refit = False,
)
results = grid_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
param_grid = param_grid,
lags_grid = lags_grid,
cv = cv,
metric = ['mean_absolute_error', mean_squared_error, custom_metric],
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
results.head(4)
| lags | lags_label | params | mean_absolute_error | mean_squared_error | custom_metric | max_depth | n_estimators | |
|---|---|---|---|---|---|---|---|---|
| 0 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 5, 'n_estimators': 100} | 0.183594 | 0.043875 | 0.070472 | 5 | 100 |
| 1 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 10, 'n_estimators': 100} | 0.183594 | 0.043875 | 0.070472 | 10 | 100 |
| 2 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 15, 'n_estimators': 100} | 0.183594 | 0.043875 | 0.070472 | 15 | 100 |
| 3 | [1, 2, 3, 20] | [1, 2, 3, 20] | {'max_depth': 15, 'n_estimators': 100} | 0.184901 | 0.044074 | 0.068182 | 15 | 100 |
The results table includes one column per metric, and the configurations are ranked by the first metric in the list (mean_absolute_error). This metric selects lags [1, 2, 3], although custom_metric would have preferred lags [1, 2, 3, 20] (0.068182 vs 0.070472). When return_best = True, the order of the metrics therefore determines the final model.
Compare multiple estimators¶
The search process can be easily extended to compare several machine learning models. This can be achieved by using a simple for loop that iterates over each estimator and applies the desired function (for example, grid_search_forecaster). This approach allows for a more thorough exploration and can help you select the best model.
# Models to compare
# ==============================================================================
models = [
RandomForestRegressor(random_state=123),
LGBMRegressor(random_state=123, verbose=-1),
Ridge()
]
# Hyperparameters to search for each model
param_grids = {
'RandomForestRegressor': {'n_estimators': [50, 100], 'max_depth': [5, 15]},
'LGBMRegressor': {'n_estimators': [20, 50], 'max_depth': [5, 10]},
'Ridge': {'alpha': [0.01, 0.1, 1]}
}
# Lags used as predictors
lags_grid = [3, 5]
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]),
refit = False,
)
results_list = []
for model in models:
model_name = type(model).__name__
print(f"Grid search for estimator: {model_name}")
print("-------------------------")
forecaster = ForecasterRecursive(
estimator = model,
lags = 3 # Placeholder, the value will be overwritten
)
results = grid_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
param_grid = param_grids[model_name],
lags_grid = lags_grid,
cv = cv,
metric = 'mean_squared_error',
return_best = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Create a column with model name
results['model'] = model_name
results_list.append(results)
df_results = pd.concat(results_list, ignore_index=True)
df_results = df_results.sort_values(by='mean_squared_error', ignore_index=True)
df_results.head(10)
Grid search for estimator: RandomForestRegressor -------------------------
Grid search for estimator: LGBMRegressor -------------------------
Grid search for estimator: Ridge -------------------------
| lags | lags_label | params | mean_squared_error | max_depth | n_estimators | model | alpha | |
|---|---|---|---|---|---|---|---|---|
| 0 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 15, 'n_estimators': 50} | 0.015294 | 15.0 | 50.0 | RandomForestRegressor | NaN |
| 1 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 5, 'n_estimators': 50} | 0.018543 | 5.0 | 50.0 | RandomForestRegressor | NaN |
| 2 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'max_depth': 5, 'n_estimators': 50} | 0.033746 | 5.0 | 50.0 | RandomForestRegressor | NaN |
| 3 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 15, 'n_estimators': 100} | 0.035851 | 15.0 | 100.0 | RandomForestRegressor | NaN |
| 4 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'max_depth': 5, 'n_estimators': 100} | 0.037824 | 5.0 | 100.0 | RandomForestRegressor | NaN |
| 5 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'max_depth': 15, 'n_estimators': 100} | 0.038111 | 15.0 | 100.0 | RandomForestRegressor | NaN |
| 6 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 10, 'n_estimators': 50} | 0.045423 | 10.0 | 50.0 | LGBMRegressor | NaN |
| 7 | [1, 2, 3] | [1, 2, 3] | {'max_depth': 5, 'n_estimators': 50} | 0.045423 | 5.0 | 50.0 | LGBMRegressor | NaN |
| 8 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'max_depth': 10, 'n_estimators': 50} | 0.045519 | 10.0 | 50.0 | LGBMRegressor | NaN |
| 9 | [1, 2, 3, 4, 5] | [1, 2, 3, 4, 5] | {'max_depth': 5, 'n_estimators': 50} | 0.045519 | 5.0 | 50.0 | LGBMRegressor | NaN |
The best configuration is a RandomForestRegressor with lags 1 to 3, max_depth = 15 and n_estimators = 50, with a validation mean squared error of 0.015325. This is about a third of the best error obtained with LGBMRegressor (0.045423), and none of the Ridge configurations appear among the 10 best.
However, this result should be interpreted with caution. The errors of the random forest vary considerably between similar configurations: with the same lags and max_depth = 15, increasing the number of trees from 50 to 100 raises the error from 0.015325 to 0.028400, although adding trees to a random forest usually has little effect. The validation set contains only 60 observations (5 folds of 12 months), so part of these differences may be due to chance (for example, to the random seed used to build the trees), and the error of the best configuration is an optimistic estimate. Before adopting the selected model, evaluate it on the test set, as shown in the grid search section.
Saving results to file¶
The results of the grid and random search can be saved to a file by setting the output_file argument to the desired path. The results are saved in a tab-separated values (TSV) format containing the lags, hyperparameters and metrics of each configuration evaluated during the search. If a file with the same name already exists, it is overwritten.
The file is updated after each evaluation, which means that if the search is stopped in the middle of the process, the results of the configurations evaluated so far are already stored in the file. This can be useful for further analysis or to keep a record of the tuning process.
In bayesian_search_forecaster, the output_file argument works differently: it stores the log messages generated by Optuna during the optimization (one line per trial) instead of a table.
# Save results to file
# ==============================================================================
forecaster = ForecasterRecursive(
estimator = LGBMRegressor(random_state=123, verbose=-1),
lags = 10 # Placeholder, the value will be overwritten
)
# Lags used as predictors
lags_grid = [3, 10, [1, 2, 3, 20]]
# Estimator hyperparameters
param_grid = {
'n_estimators': [50, 100],
'max_depth': [5, 10, 15]
}
# Folds
cv = TimeSeriesFold(
steps = 12,
initial_train_size = len(data.loc[:end_train]),
refit = False
)
results = grid_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_val, 'y'],
param_grid = param_grid,
lags_grid = lags_grid,
cv = cv,
metric = 'mean_squared_error',
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True,
output_file = "results_grid_search.txt"
)
# Read results file
# ==============================================================================
results_file = pd.read_csv("results_grid_search.txt", sep="\t")
results_file
| lags | lags_label | params | mean_squared_error | max_depth | n_estimators | |
|---|---|---|---|---|---|---|
| 0 | [1 2 3] | [1 2 3] | {'max_depth': 5, 'n_estimators': 50} | 0.045423 | 5 | 50 |
| 1 | [1 2 3] | [1 2 3] | {'max_depth': 5, 'n_estimators': 100} | 0.043875 | 5 | 100 |
| 2 | [1 2 3] | [1 2 3] | {'max_depth': 10, 'n_estimators': 50} | 0.045423 | 10 | 50 |
| 3 | [1 2 3] | [1 2 3] | {'max_depth': 10, 'n_estimators': 100} | 0.043875 | 10 | 100 |
| 4 | [1 2 3] | [1 2 3] | {'max_depth': 15, 'n_estimators': 50} | 0.045423 | 15 | 50 |
| 5 | [1 2 3] | [1 2 3] | {'max_depth': 15, 'n_estimators': 100} | 0.043875 | 15 | 100 |
| 6 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 5, 'n_estimators': 50} | 0.051399 | 5 | 50 |
| 7 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 5, 'n_estimators': 100} | 0.047896 | 5 | 100 |
| 8 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 10, 'n_estimators': 50} | 0.051399 | 10 | 50 |
| 9 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 10, 'n_estimators': 100} | 0.047896 | 10 | 100 |
| 10 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 15, 'n_estimators': 50} | 0.051399 | 15 | 50 |
| 11 | [ 1 2 3 4 5 6 7 8 9 10] | [ 1 2 3 4 5 6 7 8 9 10] | {'max_depth': 15, 'n_estimators': 100} | 0.047896 | 15 | 100 |
| 12 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 5, 'n_estimators': 50} | 0.046221 | 5 | 50 |
| 13 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 5, 'n_estimators': 100} | 0.044074 | 5 | 100 |
| 14 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 10, 'n_estimators': 50} | 0.046221 | 10 | 50 |
| 15 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 10, 'n_estimators': 100} | 0.044074 | 10 | 100 |
| 16 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 15, 'n_estimators': 50} | 0.046221 | 15 | 50 |
| 17 | [ 1 2 3 20] | [ 1 2 3 20] | {'max_depth': 15, 'n_estimators': 100} | 0.044074 | 15 | 100 |
The lags are stored in the file as text (for example, [1 2 3]), so they must be converted back to a list of integers before being reused to create a forecaster, as shown below.
# Convert the lags stored as text back to lists of integers
# ==============================================================================
results_file['lags'] = results_file['lags'].apply(
lambda x: [int(lag) for lag in x.strip('[]').split()]
)
results_file['lags'].iloc[[0, 6, 12]] # One row of each lag configuration
0 [1, 2, 3] 6 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] 12 [1, 2, 3, 20] Name: lags, dtype: object