If you like Skforecast , please give us a star on GitHub! ⭐️
Skforecast: time series forecasting with Python and Scikit-learn
More about forecasting
- ARIMA and SARIMAX models with python
- Forecasting time series with gradient boosting: Skforecast, XGBoost, LightGBM and CatBoost
- Modelling time series trend with tree based models
- Forecasting energy demand with machine learning
- Prediction intervals in forecasting models
- Intermittent demand forecasting with skforecast
- Global Forecasting Models: Comparative Analysis of Single and Multi-Series Forecasting Modeling
- Forecasting web traffic with machine learning and Python
- Bitcoin price prediction with Python, when the past does not repeat itself
- Stacking ensemble of machine learning models to improve forecasting
Introduction¶
A time series is a succession of chronologically ordered data spaced at equal or unequal intervals. The forecasting process consists of predicting the future value of a time series, either by modeling the series solely based on its past behavior (autoregressive) or by using other external variables.
This guide explores the use of scikit-learn regression models for time series forecasting. Specifically, it introduces skforecast, an intuitive library equipped with essential classes and functions to customize any Scikit-learn regression model to effectively address forecasting challenges.
✎ Note
This document serves as an introductory guide to machine learning based forecasting using skforecast. For more advanced and detailed examples, please explore: These resources delve deeper into diverse applications, offering insights and practical demonstrations of advanced techniques in time series forecasting using machine learning methodologies.Machine learning for forecasting¶
In order to apply machine learning models to forecasting problems, the time series has to be transformed into a matrix in which each value is related to the time window (lags) that precedes it.
In a time series context, a lag with respect to a time step $t$ is defined as the values of the series at previous time steps. For example, lag 1 is the value at time step $t − 1$ and lag $m$ is the value at time step $t − m$.
This type of transformation also allows to include additional variables.
Once data have been rearranged into the new shape, any regression model can be trained to predict the next value (step) of the series. During model training, every row is considered a separate data instance, where values at lags 1, 2, ... $p$ are considered predictors for the target quantity of the time series at time step $t + 1$.
Multi-Step Time Series Forecasting¶
When working with time series, it is seldom needed to predict only the next element in the series ($t_{+1}$). Instead, the most common goal is to predict a whole future interval (($t_{+1}$), ..., ($t_{+n}$)) or a far point in time ($t_{+n}$). Several strategies allow generating this type of prediction.
Recursive multi-step forecasting¶
Since the value $t_{n-1}$ is required to predict $t_{n}$, and $t_{n-1}$ is unknown, a recursive process is applied in which, each new prediction, is based on the previous one. This process is known as recursive forecasting or recursive multi-step forecasting and can be easily generated with the ForecasterAutoreg and ForecasterAutoregCustom classes.
Direct multi-step forecasting¶
Direct multi-step forecasting consists of training a different model for each step of the forecast horizon. For example, to predict the next 5 values of a time series, 5 different models are trained, one for each step. As a result, the predictions are independent of each other.
The main complexity of this approach is to generate the correct training matrices for each model. The ForecasterAutoregDirect class of the skforecast library automates this process. It is also important to bear in mind that this strategy has a higher computational cost since it requires the train of multiple models. The following diagram shows the process for a case in which the response variable and two exogenous variables are available.
Multiple output forecasting¶
Some machine learning models, such as long short-term memory (LSTM) neural networks, can predict multiple values of a sequence simultaneously (one-shot). This strategy is not currently implemented in the skforecast library, but is expected to be included in future versions.
Libraries¶
The libraries used in this document are:
# Data manipulation
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset
# Plots
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 10
# Modeling and Forecasting
# ==============================================================================
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.preprocessing import StandardScaler
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregCustom import ForecasterAutoregCustom
from skforecast.ForecasterAutoregDirect import ForecasterAutoregDirect
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from skforecast.utils import save_forecaster
from skforecast.utils import load_forecaster
# Warnings configuration
# ==============================================================================
import warnings
warnings.filterwarnings('once')
Data¶
A time series is available with the monthly expenditure (millions of dollars) on corticosteroid drugs that the Australian health system had between 1991 and 2008. It is intended to create an autoregressive model capable of predicting future monthly expenditures. The data used in the examples of this document have been obtained from the magnificent book Forecasting: Principles and Practice by Rob J Hyndman and George Athanasopoulos.
# Data download
# ==============================================================================
data = fetch_dataset(name='h2o_exog', raw=True)
The column date has been stored as a string. To convert it to datetime the pd.to_datetime() function can be use. Once in datetime format, and to make use of Pandas functionalities, it is set as an index. Also, since the data is monthly, the frequency is set as Monthly Started 'MS'.
# Data preparation
# ==============================================================================
data = data.rename(columns={'fecha': 'date'})
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = data.set_index('date')
data = data.asfreq('MS')
data = data.sort_index()
data.head()
When using the asfreq() method in Pandas, any gaps in the time series will be filled with NaN values to match the specified frequency. Therefore, it is essential to check for any missing values that may occur after this transformation.
# Missing values
# ==============================================================================
print(f'Number of rows with missing values: {data.isnull().any(axis=1).mean()}')
Although it is unnecessary, since a frequency has been established, it is possible to verify that the time series is complete.
# Verify that a temporary index is complete
# ==============================================================================
(data.index == pd.date_range(start=data.index.min(),
end=data.index.max(),
freq=data.index.freq)).all()
# Fill gaps in a temporary index
# ==============================================================================
# data.asfreq(freq='30min', fill_value=np.nan)
The last 36 months are used as the test set to evaluate the predictive capacity of the model.
# Split data into train-test
# ==============================================================================
steps = 36
data_train = data[:-steps]
data_test = data[-steps:]
print(f"Train dates : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Test dates : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
ax.legend();
Recursive multi-step forecasting¶
ForecasterAutoreg¶
With the ForecasterAutoreg class, a forecasting model is created and trained using a RandomForestRegressor regressor with a time window of 6 lags. This means that the model uses the previous 6 months as predictors.
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = RandomForestRegressor(random_state=123),
lags = 6
)
forecaster.fit(y=data_train['y'])
forecaster
Prediction¶
Once the model is trained, the test data is predicted (36 months into the future).
# Predictions
# ==============================================================================
steps = 36
predictions = forecaster.predict(steps=steps)
predictions.head(5)
# Plot predictions versus test data
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
predictions.plot(ax=ax, label='predictions')
ax.legend();
Forecasting error in the test set¶
The error that the model makes in its predictions is quantified. In this case, the metric used is the mean squared error (mse).
# Test error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions
)
print(f"Test error (MSE): {error_mse}")
Hyperparameter tuning¶
The trained ForecasterAutoreg uses a 6 lag time window and a Random Forest model with the default hyperparameters. However, there is no reason why these values are the most suitable. Skforecast provide several search strategies to find the best combination of hyperparameters and lags. In this case, the grid_search_forecaster function is used. It compares the results obtained with each combinations of hyperparameters and lags, and identify the best one.
💡 Tip
The computational cost of hyperparameter tuning depends heavily on the backtesting approach chosen to evaluate each hyperparameter combination. In general, the duration of the tuning process increases with the number of re-trains involved in the backtesting.
To effectively speed up the prototyping phase, it is highly recommended to adopt a two-step strategy. First, use refit=False during the initial search to narrow down the range of values. Then, focus on the identified region of interest and apply a tailored backtesting strategy that meets the specific requirements of the use case. For more tips on backtesting strategies, see Hyperparameter tuning and lags selection.
# Hyperparameter grid search
# ==============================================================================
steps = 36
forecaster = ForecasterAutoreg(
regressor = RandomForestRegressor(random_state=123),
lags = 12 # This value will be replaced in the grid search
)
# Candidate values for lags
lags_grid = [10, 20]
# Candidate values for regressor's hyperparameters
param_grid = {
'n_estimators': [100, 500],
'max_depth': [3, 5, 10]
}
results_grid = grid_search_forecaster(
forecaster = forecaster,
y = data_train['y'],
param_grid = param_grid,
lags_grid = lags_grid,
steps = steps,
refit = False,
metric = 'mean_squared_error',
initial_train_size = int(len(data_train)*0.5),
fixed_train_size = False,
return_best = True,
n_jobs = 'auto',
verbose = False
)
# Search results
# ==============================================================================
results_grid
The best results are obtained using a time window of 20 lags and a Random Forest set up of {'max_depth': 3, 'n_estimators': 500}.
Final model¶
Finally, a ForecasterAutoreg is trained with the optimal configuration found. This step is not necessary if return_best = True is specified in the grid_search_forecaster function.
# Create and train forecaster with the best hyperparameters and lags found
# ==============================================================================
regressor = RandomForestRegressor(n_estimators=500, max_depth=3, random_state=123)
forecaster = ForecasterAutoreg(
regressor = regressor,
lags = 20
)
forecaster.fit(y=data_train['y'])
# Predictions
# ==============================================================================
predictions = forecaster.predict(steps=steps)
# Plot predictions versus test data
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
predictions.plot(ax=ax, label='predictions')
ax.legend();
# Test error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions
)
print(f"Test error (MSE): {error_mse}")
The optimal combination of hyperparameters significantly reduces test error.
Backtesting¶
To obtain a robust estimate of the model's predictive capacity, a backtesting process is carried out. The process of backtesting consists of evaluating the performance of a predictive model by applying it retrospectively to historical data. Therefore, it is a special type of cross-validation applied to the previous period(s).
✎ Note
To ensure an accurate evaluation of your model and gain confidence in its predictive performance on new data, it is critical to employ an appropriate backtesting strategy. Factors such as use case characteristics, available computing resources and time intervals between predictions need to be considered to determine which strategy to use.
In general, the more closely the backtesting process resembles the actual scenario in which the model is used, the more reliable the estimated metric will be. For more information about backtesting, visit Which strategy should I use?.
Backtesting with refit and increasing training size (fixed origin)
The model is trained each time before making predictions. With this configuration, the model uses all the data available so far. It is a variation of the standard cross-validation but, instead of making a random distribution of the observations, the training set increases sequentially, maintaining the temporal order of the data.
Backtesting with refit and fixed training size (rolling origin)
A technique similar to the previous one but, in this case, the forecast origin rolls forward, therefore, the size of training remains constant. This is also known as time series cross-validation or walk-forward validation.
Backtesting with intermittent refit
The model is retrained every $n$ iterations of predictions.
💡 Tip
This strategy usually achieves a good balance between the computational cost of retraining and avoiding model degradation.
Backtesting without refit
After an initial train, the model is used sequentially without updating it and following the temporal order of the data. This strategy has the advantage of being much faster since the model is trained only once. However, the model does not incorporate the latest data available, so it may lose predictive capacity over time.
Skforecast library has multiple backtesting strategies implemented. Regardless of which one is used, it is important not to include test data in the search process to avoid overfitting problems.
For this example, a backtesting with refit strategy is followed. Internally, the process that the function applies is:
In the first iteration, the model is trained with the observations selected for the initial training (in this case, 87). Then, the next 36 observations are used to validate the predictions of this first model.
In the second iteration, the model is retrained by adding, to the initial training set, the previous 36 validation observations (87 + 36). In the same way, the next 36 observations are established as the new validation set.
This process is repeated until all available observations are used. Following this strategy, the training set increases in each iteration with as many observations as steps are being predicted.
# Backtesting
# ==============================================================================
steps = 36
n_backtesting = 36*3 # The last 9 years are separated for the backtest
metric, predictions_backtest = backtesting_forecaster(
forecaster = forecaster,
y = data['y'],
initial_train_size = len(data) - n_backtesting,
fixed_train_size = False,
steps = steps,
metric = 'mean_squared_error',
refit = True,
verbose = True,
show_progress = True
)
print(f"Backtest metric (MSE): {metric}")
# Plot backtest predictions vs real values
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data.loc[predictions_backtest.index, 'y'].plot(ax=ax, label='test')
predictions_backtest.plot(ax=ax, label='predictions')
ax.legend();
Model explanaibility (Feature importance)¶
Due to the complex nature of many modern machine learning models, such as ensemble methods, they often function as black boxes, making it difficult to understand why a particular prediction was made. Explanability techniques aim to demystify these models, providing insight into their inner workings and helping to build trust, improve transparency, and meet regulatory requirements in various domains. Enhancing model explainability not only helps to understand model behavior, but also helps to identify biases, improve model performance, and enable stakeholders to make more informed decisions based on machine learning insights.
Skforecast is compatible with some of the most popular model explainability methods: model-specific feature importances, SHAP values, and partial dependence plots.
# Predictors importances
# ==============================================================================
forecaster.get_feature_importances()
⚠ Warning
Theget_feature_importances() method will only return values if the forecaster's regressor has either the coef_ or feature_importances_ attribute, which is the default in scikit-learn.
Forecasting with exogenous variables¶
In the previous example, only lags of the predicted variable itself were used as predictors. In certain scenarios, it is possible to have information about other variables, whose future value is known, so could serve as additional predictors in the model.
Continuing with the previous example, a new variable is simulated whose behavior is correlated with the modeled time series and which is to be included as a predictor.
# Data download
# ==============================================================================
data = fetch_dataset(name='h2o_exog', raw=True)
# Data preparation
# ==============================================================================
data = data.rename(columns={'fecha': 'date'})
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = data.set_index('date')
data = data.asfreq('MS')
data = data.sort_index()
fig, ax = plt.subplots(figsize=(6, 2.5))
data['y'].plot(ax=ax, label='y')
data['exog_1'].plot(ax=ax, label='exogenous variable')
ax.legend();
# Split data into train-test
# ==============================================================================
steps = 36
data_train = data[:-steps]
data_test = data[-steps:]
print(f"Train dates : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Test dates : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = RandomForestRegressor(random_state=123),
lags = 8
)
forecaster.fit(y=data_train['y'], exog=data_train['exog_1'])
forecaster
Since the ForecasterAutoreg has been trained with an exogenous variable, the value of this variable must be passed to predict(). The future information about the exogenous variable must be available when making predictions.
# Predictions
# ==============================================================================
predictions = forecaster.predict(steps=steps, exog=data_test['exog_1'])
# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
predictions.plot(ax=ax, label='predictions')
ax.legend();
# Test error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions
)
print(f"Test error (MSE): {error_mse}")
Forecasting with custom predictors¶
In addition to the lags, it may be interesting to incorporate other characteristics of the time series in some scenarios. For example, the moving average of the last n values could be used to capture the trend of the series.
The ForecasterAutoregCustom class behaves very similar to the ForecasterAutoreg class seen in the previous sections, but with the difference that it is the user who defines the function used to create the predictors.
The first example of the document is repeated. In this case, the predictors are the first 10 lags and the moving average of the values of the last 20 months.
# Data download
# ==============================================================================
data = fetch_dataset(name='h2o_exog', raw=True)
# Data preparation
# ==============================================================================
data = data.rename(columns={'fecha': 'date'})
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = data.set_index('date')
data = data.asfreq('MS')
data = data.sort_index()
# Split data into train-test
# ==============================================================================
steps = 36
data_train = data[:-steps]
data_test = data[-steps:]
print(f"Train dates : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Test dates : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
A ForecasterAutoregCustom is created and trained from a RandomForestRegressor regressor. The create_predictor() function, which calculates the first 10 lags and the moving average of the last 20 values, is used to create the predictors.
# Function to calculate predictors from time series
# ==============================================================================
def custom_predictors(y):
"""
Create first 10 lags of a time series.
Calculate moving average with window 20.
"""
lags = y[-1:-11:-1] # window size needed = 10
mean = np.mean(y[-20:]) # window size needed = 20
predictors = np.hstack([lags, mean])
return predictors
⚠ Warning
Thewindow_size parameter specifies the size of the data window that fun_predictors uses to generate each row of predictors. Choosing the appropriate value for this parameter is crucial to avoid losing data when constructing the training matrices.
In this case, the
window_size required by the mean is the largest (most restrictive), so window_size = 20.
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoregCustom(
regressor = RandomForestRegressor(random_state=123),
fun_predictors = custom_predictors,
window_size = 20
)
forecaster.fit(y=data_train['y'])
forecaster
It is possible to access the custom function used to create the predictors.
# Custom function source code
# ==============================================================================
print(forecaster.source_code_fun_predictors)
Using the method create_train_X_y, is it posible to acces the matrices that are created internally in the training process.
# Training matrix
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=data_train['y'])
display(X_train.head(5))
display(y_train.head(5))
# Predictions
# ==============================================================================
steps = 36
predictions = forecaster.predict(steps=steps)
# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
predictions.plot(ax=ax, label='predictions')
ax.legend();
# Test error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions
)
print(f"Test error (MSE): {error_mse}")
Direct multi-step forecasting¶
The ForecasterAutoreg and ForecasterAutoregCustom models follow a recursive prediction strategy in which each new prediction builds on the previous one. An alternative is to train a model for each of the steps to be predicted. This strategy, commonly known as direct multi-step forecasting, is computationally more expensive than recursive since it requires training several models. However, in some scenarios, it achieves better results. These kinds of models can be obtained with the ForecasterAutoregDirect class and can include one or multiple exogenous variables.
⚠ Warning
`ForecasterAutoregDirect` may require very long training times, as one model is fitted for each step.
ForecasterAutoregDirect¶
Unlike when using ForecasterAutoreg or ForecasterAutoregCustom, the number of steps to be predicted must be indicated in the ForecasterAutoregDirect. This means that it is not possible to predict steps beyond the value defined at their creation when executing the predict() method.
For this example, a linear model with Ridge penalty is used as a regressor. These models require the predictors to be standardized, so it is combined with a StandardScaler.
For more detailed documentation how to use transformers and pipelines, visit: skforecast with scikit-learn and transformers pipelines.
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoregDirect(
regressor = Ridge(random_state=123),
steps = 36,
lags = 8,
transformer_y = StandardScaler()
)
forecaster
# Hyperparameter Grid search
# ==============================================================================
from skforecast.exceptions import LongTrainingWarning
warnings.simplefilter('ignore', category=LongTrainingWarning)
forecaster = ForecasterAutoregDirect(
regressor = Ridge(random_state=123),
steps = 36,
lags = 8, # This value will be replaced in the grid search
transformer_y = StandardScaler()
)
param_grid = {'alpha': np.logspace(-5, 5, 10)}
lags_grid = [5, 12, 20]
results_grid = grid_search_forecaster(
forecaster = forecaster,
y = data_train['y'],
param_grid = param_grid,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_squared_error',
initial_train_size = int(len(data_train)*0.5),
fixed_train_size = False,
return_best = True,
n_jobs = 'auto',
verbose = False
)
# Grid Search results
# ==============================================================================
results_grid.head()
The best results are obtained using a time window of 12 lags and a Ridge setting {'alpha': 0.021544}.
# Predictions
# ==============================================================================
predictions = forecaster.predict()
# Plot predictions versus test data
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_train['y'].plot(ax=ax, label='train')
data_test['y'].plot(ax=ax, label='test')
predictions.plot(ax=ax, label='predictions')
ax.legend();
# Test error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions
)
print(f"Test error (MSE) {error_mse}")
Prediction intervals¶
A prediction interval defines the interval within which the true value of the target variable can be expected to be found with a given probability. Rob J Hyndman and George Athanasopoulos, in their book Forecasting: Principles and Practice, list multiple ways to estimate prediction intervals, most of which require that the residuals (errors) of the model to be normally distributed. If this cannot be assumed, one can resort to bootstrapping, which requires only that the residuals be uncorrelated. This is one of the methods available in skforecast. A more detailed explanation of prediction intervals can be found in the Probabilistic forecasting: prediction intervals and prediction distribution user guide.
# Data download
# ==============================================================================
data = fetch_dataset(name='h2o_exog', raw=True)
# Data preparation
# ==============================================================================
data = data.rename(columns={'fecha': 'date'})
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = data.set_index('date')
data = data.asfreq('MS')
data = data.sort_index()
# Split data into train-test
# ==============================================================================
steps = 36
data_train = data[:-steps]
data_test = data[-steps:]
print(f"Train dates : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Test dates : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = Ridge(alpha=0.1, random_state=765),
lags = 15
)
forecaster.fit(y=data_train['y'])
# Prediction intervals
# ==============================================================================
predictions = forecaster.predict_interval(
steps = steps,
interval = [1, 99],
n_boot = 500
)
predictions.head(5)
# Prediction error
# ==============================================================================
error_mse = mean_squared_error(
y_true = data_test['y'],
y_pred = predictions.iloc[:, 0]
)
print(f"Test error (MSE): {error_mse}")
# Plot forecasts with prediction intervals
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data_test['y'].plot(ax=ax, label='test')
predictions['pred'].plot(ax=ax, label='prediction')
ax.fill_between(
predictions.index,
predictions['lower_bound'],
predictions['upper_bound'],
color = 'red',
alpha = 0.2
)
ax.legend();
# Backtest with prediction intervals
# ==============================================================================
n_backtesting = 36*3 # The last 9 years are separated for backtesting
steps = 36
forecaster = ForecasterAutoreg(
regressor = Ridge(alpha=0.1, random_state=765),
lags = 15
)
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data['y'],
initial_train_size = len(data) - n_backtesting,
fixed_train_size = False,
steps = steps,
metric = 'mean_squared_error',
refit = True,
interval = [1, 99],
n_boot = 100,
n_jobs = 'auto',
verbose = True
)
print(f"Test error (MSE): {error_mse}")
# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
data.loc[predictions.index, 'y'].plot(ax=ax, label='test')
predictions['pred'].plot(ax=ax, label='predicciones')
ax.fill_between(
predictions.index,
predictions['lower_bound'],
predictions['upper_bound'],
color = 'red',
alpha = 0.2
)
ax.legend();
# Predicted interval coverage
# ==============================================================================
inside_interval = np.where(
(data.loc[predictions.index, 'y'] >= predictions['lower_bound']) & \
(data.loc[predictions.index, 'y'] <= predictions['upper_bound']),
True,
False
)
coverage = inside_interval.mean()
print(f"Predicted interval coverage: {round(100*coverage, 2)} %")
Custom metric¶
In the backtesting (backtesting_forecaster) and hyperparameter optimization (grid_search_forecaster) processes, besides the frequently used metrics such as mean_squared_error or mean_absolute_error, it is possible to use any custom function as long as:
It includes the arguments:
y_true: true values of the series.y_pred: predicted values.
It returns a numeric value (
floatorint).The metric is reduced as the model improves. Only applies in the
grid_search_forecasterfunction ifreturn_best=True(train the forecaster with the best model).
It allows evaluating the predictive capability of the model in a wide range of scenarios, for example:
Consider only certain months, days, hours...
Consider only dates that are holidays.
Consider only the last step of the predicted horizon.
The following example shows how to forecast a 12-month horizon but considering only the last 3 months of each year to calculate the interest metric.
# Custom metric
# ==============================================================================
def custom_metric(y_true, y_pred):
"""
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_absolute_error(y_true[mask], y_pred[mask])
return metric
# Backtesting
# ==============================================================================
steps = 36
n_backtesting = 36*3 # The last 9 years are separated for backtesting
metric, predictions_backtest = backtesting_forecaster(
forecaster = forecaster,
y = data['y'],
initial_train_size = len(data) - n_backtesting,
fixed_train_size = False,
steps = steps,
refit = True,
metric = custom_metric,
n_jobs = 'auto',
verbose = True
)
print(f"Backtest error: {metric}")
Save and load models¶
Skforecast models can be stored and loaded from disck using pickle or joblib library. To simply the process, two functions are available: save_forecaster and load_forecaster, simple example is shown below. For more detailed documentation, visit: skforecast save and load forecaster.
# Create forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(RandomForestRegressor(random_state=123), lags=3)
forecaster.fit(y=data['y'])
forecaster.predict(steps=3)
# Save forecaster
# ==============================================================================
save_forecaster(forecaster, file_name='forecaster.joblib', verbose=False)
# Load forecaster
# ==============================================================================
forecaster_loaded = load_forecaster('forecaster.joblib')
# Predict
# ==============================================================================
forecaster_loaded.predict(steps=3)
⚠ Warning
When usig forecasters with custom functions, such asForecasterAutoregCustom, the function used to create the predictors must be defined before loading the forecaster. Otherwise, an error will be raised. Therefore, it is recommended to save the function in a separate file and import it before loading the forecaster.
Use forecaster in production¶
In projects related to forecasting, it is common to generate a model after the experimentation and development phase. For this model to have a positive impact on the business, it must be able to be put into production and generate forecasts from time to time with which to decide. This need has widely guided the development of the skforecast library.
Suppose predictions have to be generated on a weekly basis, for example, every Monday. By default, when using the predict method on a trained forecaster object, predictions start right after the last training observation. Therefore, the model could be retrained weekly, just before the first prediction is needed, and call its predict method.
This strategy, although simple, may not be possible to use for several reasons:
Model training is very expensive and cannot be run as often.
The history with which the model was trained is no longer available.
The prediction frequency is so high that there is no time to train the model between predictions.
In these scenarios, the model must be able to predict at any time, even if it has not been recently trained.
Every model generated using skforecast has the last_window argument in its predict method. Using this argument, it is possible to provide only the past values needs to create the autoregressive predictors (lags) and thus, generate the predictions without the need to retrain the model.
For more detailed documentation, visit: skforecast forecaster in production.
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = RandomForestRegressor(random_state=123),
lags = 6
)
forecaster.fit(y=data_train['y'])
Since the model uses the last 6 lags as predictors, last_window must contain at least the 6 values previous to the moment where the prediction starts.
# Predict with last window
# ==============================================================================
last_window = data_test['y'][-6:]
forecaster.predict(last_window=last_window, steps=4)
If the forecaster uses exogenous variables, besides last_window, the argument exog must contain the future values of the exogenous variables.
Other functionalities¶
There are too many functionalities offered by skforecast to be shown in a single document. The reader is encouraged to read the following list of links to learn more about the library.
Examples and tutorials¶
We also recommend to explore our collection of examples and tutorials where you will find numerous use cases that will give you a practical insight into how to use this powerful library. Examples and tutorials.
Session information¶
import session_info
session_info.show(html=False)
Bibliography¶
Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia.
Time Series Analysis and Forecasting with ADAM Ivan Svetunkov
Joseph, M. (2022). Modern time series forecasting with Python: Explore industry-ready time series forecasting using modern machine learning and Deep Learning. Packt Publishing.
Python for Finance: Mastering Data-Driven Finance
Citation¶
How to cite this document
If you use this document or any part of it, please acknowledge the source, thank you!
Skforecast: time series forecasting with Python and Scikit-learn by Joaquín Amat Rodrigo and Javier Escobar Ortiz, available under Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) at https://cienciadedatos.net/documentos/py27-time-series-forecasting-python-scikitlearn.html
How to cite skforecast
If you use skforecast for a scientific publication, we would appreciate it if you cite the published software.
Zenodo:
Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2023). skforecast (v0.11.0). Zenodo. https://doi.org/10.5281/zenodo.8382788
APA:
Amat Rodrigo, J., & Escobar Ortiz, J. (2023). skforecast (Version 0.11.0) [Computer software]. https://doi.org/10.5281/zenodo.8382788
BibTeX:
@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.11.0}, month = {11}, year = {2023}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }
Did you like the article? Your support is important
Website maintenance has high cost, your contribution will help me to continue generating free educational content. Many thanks! 😊
This work by Joaquín Amat Rodrigo is licensed under a Attribution-NonCommercial-ShareAlike 4.0 International.