Window and custom features¶

When forecasting time series data, it may be useful to consider additional characteristics beyond just the lagged values. For example, the moving average of the previous n values may help to capture the trend in the series. The window_features argument allows the inclusion of additional predictors created with the previous values of the series.

Libraries and data¶

In [1]:

# Libraries
# ==============================================================================
import numpy as np
import pandas as pd
from lightgbm import LGBMRegressor
from skforecast.datasets import fetch_dataset
from skforecast.preprocessing import RollingFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.recursive import ForecasterRecursiveMultiSeries

# Libraries # ============================================================================== import numpy as np import pandas as pd from lightgbm import LGBMRegressor from skforecast.datasets import fetch_dataset from skforecast.preprocessing import RollingFeatures from skforecast.recursive import ForecasterRecursive from skforecast.recursive import ForecasterRecursiveMultiSeries

In [2]:

# Data
# ==============================================================================
data = fetch_dataset(name="h2o", raw=False)
data.index.name = 'datetime'
data = data.rename(columns={'x': 'y'})
y = data['y']
y

# Data # ============================================================================== data = fetch_dataset(name="h2o", raw=False) data.index.name = 'datetime' data = data.rename(columns={'x': 'y'}) y = data['y'] y

h2o
---
Monthly expenditure ($AUD) on corticosteroid drugs that the Australian health
system had between 1991 and 2008.
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.
Shape of the dataset: (204, 1)

Out[2]:

datetime
1991-07-01    0.429795
1991-08-01    0.400906
1991-09-01    0.432159
1991-10-01    0.492543
1991-11-01    0.502369
                ...   
2008-02-01    0.761822
2008-03-01    0.649435
2008-04-01    0.827887
2008-05-01    0.816255
2008-06-01    0.762137
Freq: MS, Name: y, Length: 204, dtype: float64

RollingFeatures¶

The RollingFeatures class available is skforecast allows the creation of some of the most commonly used predictors:

• 'mean': the mean of the previous n values.
• 'std': the standard deviation of the previous n values.
• 'min': the minimum of the previous n values.
• 'max': the maximum of the previous n values.
• 'sum': the sum of the previous n values.
• 'median': the median of the previous n values.
• 'ratio_min_max': the ratio between the minimum and maximum of the previous n values.
• 'coef_variation': the coefficient of variation of the previous n values.

The user can specify these predictors by passing a list of strings to the stats argument. The user can also specify the window size for each of these predictors by passing a list of integers to the window_size argument, if the value is the same for all the predictors, the user can pass a single integer.

The following example demonstrates how to use the RollingFeatures class to include rolling statistics (mean, minimum, and maximum values). Here, the rolling mean is computed with a window size of 10 and 20, while the minimum and maximum values use a window size of 10.

In [3]:

# Window features
# ==============================================================================
window_features = RollingFeatures(
                      stats        = ['mean', 'mean', 'min', 'max'],
                      window_sizes = [20, 10, 10, 10]
                  )

# Window features # ============================================================================== window_features = RollingFeatures( stats = ['mean', 'mean', 'min', 'max'], window_sizes = [20, 10, 10, 10] )

In [4]:

# Create and fit forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor       = LGBMRegressor(random_state=123, verbose=-1),
                 lags            = 3,
                 window_features = window_features
             )

forecaster.fit(y=y)
forecaster

# Create and fit forecaster # ============================================================================== forecaster = ForecasterRecursive( regressor = LGBMRegressor(random_state=123, verbose=-1), lags = 3, window_features = window_features ) forecaster.fit(y=y) forecaster

Out[4]:

ForecasterRecursive

General Information

• Regressor: LGBMRegressor
• Lags: [1 2 3]
• Window features: ['roll_mean_20', 'roll_mean_10', 'roll_min_10', 'roll_max_10']
• Window size: 20
• Exogenous included: False
• Weight function included: False
• Differentiation order: None
• Creation date: 2024-11-12 15:20:38
• Last fit date: 2024-11-12 15:20:39
• Skforecast version: 0.14.0
• Python version: 3.11.10
• 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('2008-06-01 00:00:00')]
• Training index type: DatetimeIndex
• Training index frequency: MS

Regressor Parameters

{'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0, 'importance_type': 'split', 'learning_rate': 0.1, 'max_depth': -1, '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

{}

🛈 API Reference    🗎 User Guide

In [5]:

# Predict
# ==============================================================================
steps = 36
predictions = forecaster.predict(steps=steps)
predictions.head(3)

# Predict # ============================================================================== steps = 36 predictions = forecaster.predict(steps=steps) predictions.head(3)

Out[5]:

2008-07-01    0.784687
2008-08-01    0.907246
2008-09-01    1.078849
Freq: MS, Name: pred, dtype: float64

The window size needed by this Forecaster is the maximum of the window sizes between the lagged values and the rolling features. In this case, this value is 20.

In [6]:

# Predict
# ==============================================================================
print("Window size lags            : ", forecaster.max_lag)
print("Window size window features : ", forecaster.max_size_window_features)
print("Window size Forecaster      : ", forecaster.window_size)

# Predict # ============================================================================== print("Window size lags : ", forecaster.max_lag) print("Window size window features : ", forecaster.max_size_window_features) print("Window size Forecaster : ", forecaster.window_size)

Window size lags            :  3
Window size window features :  20
Window size Forecaster      :  20

Additional arguments can be passed to the RollingFeatures:

• features_names: By default, feature names follow the pattern roll_<stat>_<window_size>. For instance, a rolling mean with a window size of 20 is named roll_mean_20. Users can also assign custom names to each feature using a list of strings.

• min_periods: allows specifying the minimum number of observations required to compute the statistics during the training matrix generation (same as the min_periods argument of pandas rolling). It can be a single integer or a list of integers, one for each statistic.

• fill_na: define the strategy for handling missing values during the training matrix generation. Available methods are: 'mean', 'median', 'ffill', 'bfill', or a float value.

By inspecting the training matrices, it is possible to check that the rolling features have been correctly included.

In [7]:

# Training matrices used internally to fit the regressor
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=y)
X_train

# Training matrices used internally to fit the regressor # ============================================================================== X_train, y_train = forecaster.create_train_X_y(y=y) X_train

Out[7]:

	lag_1	lag_2	lag_3	roll_mean_20	roll_mean_10	roll_min_10	roll_max_10
datetime
1993-03-01	0.387554	0.751503	0.771258	0.496401	0.534009	0.361801	0.771258
1993-04-01	0.427283	0.387554	0.751503	0.496275	0.540557	0.387554	0.771258
1993-05-01	0.413890	0.427283	0.387554	0.496924	0.540893	0.387554	0.771258
1993-06-01	0.428859	0.413890	0.427283	0.496759	0.535440	0.387554	0.771258
1993-07-01	0.470126	0.428859	0.413890	0.495638	0.534906	0.387554	0.771258
...	...	...	...	...	...	...	...
2008-02-01	1.219941	1.176589	1.163534	0.980390	0.995834	0.561760	1.219941
2008-03-01	0.761822	1.219941	1.176589	0.978582	1.015840	0.745258	1.219941
2008-04-01	0.649435	0.761822	1.219941	0.966838	1.006258	0.649435	1.219941
2008-05-01	0.827887	0.649435	0.761822	0.955750	1.005253	0.649435	1.219941
2008-06-01	0.816255	0.827887	0.649435	0.946778	0.991464	0.649435	1.219941

184 rows × 7 columns

It is also possible to use the RollingFeatures class outside the forecaster to gain a deeper insight into its behaviour. The transform method computes the rolling features for a given numpy array, which is assumed to contain as many past observations as the maximum window size required to compute the features. The output is a numpy array with the rolling features.

In [8]:

# Create rolling features from a given array
# ==============================================================================
x = np.arange(20)
window_features.transform(X=x)

# Create rolling features from a given array # ============================================================================== x = np.arange(20) window_features.transform(X=x)

Out[8]:

array([ 9.5, 14.5, 10. , 19. ])

The transform_batch method is designed to transform a whole pandas Series from which multiple rolling windows can be extracted. The output is a pandas DataFrame with the rolling features.

In [9]:

# Create rolling features from a pandas series
# ==============================================================================
window_features.transform_batch(y).head(3)

# Create rolling features from a pandas series # ============================================================================== window_features.transform_batch(y).head(3)

Out[9]:

	roll_mean_20	roll_mean_10	roll_min_10	roll_max_10
datetime
1993-03-01	0.496401	0.534009	0.361801	0.771258
1993-04-01	0.496275	0.540557	0.387554	0.771258
1993-05-01	0.496924	0.540893	0.387554	0.771258

The reason for these two different data transformation methods is that the first is used during prediction, where the forecaster only has access to the last window of the series. In contrast, the second is used during training, where the forecaster has access to the entire series.

Create your custom window features¶

RollingFeatures is very useful for including some of the most commonly used predictors. However, users may need to include additional predictors that are not provided by this class. In such cases, users can create their own custom class to compute the desired features and include them in the forecaster.

The custom class must have these 2 methods:

• transform_batch: method to compute the features in batch from a pandas Series. This method will be used to compute the features during the training process. It must return a pandas DataFrame containing the rolling features.

• transform: method to compute the features from a numpy array. This method will be used to compute the features during the prediction process. It must return a numpy array containing the computed statistics.

and these 2 attributes:

• window_sizes: size of the rolling window required to compute the features. It must be a list of integers.

• features_names: list with the names of the output features. It must be a list of strings.

The follwing example shows how to create a custom class to include the rolling skewness and kurtosis with a window size of 20.

💡 Tip

If you have any doubt when creating your custom class, you can check the source code of the RollingFeatures class available in the API Reference.

In [10]:

# Custom class to create rolling skewness features
# ==============================================================================
from scipy.stats import skew


class RollingSkewness():
    """
    Custom class to create rolling skewness features.
    """

    def __init__(self, window_sizes, features_names: list = 'rolling_skewness'):
        
        if not isinstance(window_sizes, list):
            window_sizes = [window_sizes]
        self.window_sizes = window_sizes
        self.features_names = features_names

    def transform_batch(self, X: pd.Series) -> pd.DataFrame:
        
        rolling_obj = X.rolling(window=self.window_sizes[0], center=False, closed='left')
        rolling_skewness = rolling_obj.skew()
        rolling_skewness = pd.DataFrame({
                               self.features_names: rolling_skewness
                           }).dropna()

        return rolling_skewness

    def transform(self, X: np.ndarray) -> np.ndarray:
        
        X = X[~np.isnan(X)]
        if len(X) > 0:
            rolling_skewness = np.array([skew(X, bias=False)])
        else:
            rolling_skewness = np.array([np.nan])
        
        return rolling_skewness

# Custom class to create rolling skewness features # ============================================================================== from scipy.stats import skew class RollingSkewness(): """ Custom class to create rolling skewness features. """ def __init__(self, window_sizes, features_names: list = 'rolling_skewness'): if not isinstance(window_sizes, list): window_sizes = [window_sizes] self.window_sizes = window_sizes self.features_names = features_names def transform_batch(self, X: pd.Series) -> pd.DataFrame: rolling_obj = X.rolling(window=self.window_sizes[0], center=False, closed='left') rolling_skewness = rolling_obj.skew() rolling_skewness = pd.DataFrame({ self.features_names: rolling_skewness }).dropna() return rolling_skewness def transform(self, X: np.ndarray) -> np.ndarray: X = X[~np.isnan(X)] if len(X) > 0: rolling_skewness = np.array([skew(X, bias=False)]) else: rolling_skewness = np.array([np.nan]) return rolling_skewness

In [11]:

# Transform batch
# ==============================================================================
window_features = RollingSkewness(window_sizes=3)
window_features.transform_batch(y)

# Transform batch # ============================================================================== window_features = RollingSkewness(window_sizes=3) window_features.transform_batch(y)

Out[11]:

	rolling_skewness
datetime
1991-10-01	-1.696160
1991-11-01	0.897261
1991-12-01	-1.602797
1992-01-01	1.681518
1992-02-01	-0.778727
...	...
2008-02-01	1.359033
2008-03-01	-1.674974
2008-04-01	1.466482
2008-05-01	-0.747574
2008-06-01	-1.705640

201 rows × 1 columns

In [12]:

# Transform
# ==============================================================================
window_features.transform(X=np.array([6, 12, 8]))

# Transform # ============================================================================== window_features.transform(X=np.array([6, 12, 8]))

Out[12]:

array([0.93521953])

In [13]:

# Forecaster with custom rolling features
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor       = LGBMRegressor(random_state=123, n_jobs=-1, verbose=-1),
                 lags            = 3,
                 window_features = window_features
             )

forecaster.fit(y=y)
forecaster.predict(steps=5)

# Forecaster with custom rolling features # ============================================================================== forecaster = ForecasterRecursive( regressor = LGBMRegressor(random_state=123, n_jobs=-1, verbose=-1), lags = 3, window_features = window_features ) forecaster.fit(y=y) forecaster.predict(steps=5)

Out[13]:

2008-07-01    0.769870
2008-08-01    0.824683
2008-09-01    0.819595
2008-10-01    0.799849
2008-11-01    0.802432
Freq: MS, Name: pred, dtype: float64

Create your custom window features ForecasterRecursiveMultiSeries¶

When using a ForecasterRecursiveMultiSeries, the transform method must return a numpy array that calculates all the features for all the series contained in last_window at once.

In [14]:

# Data download
# ==============================================================================
data_multiseries = fetch_dataset(name="items_sales")
data_multiseries.head()

# Data download # ============================================================================== data_multiseries = fetch_dataset(name="items_sales") data_multiseries.head()

items_sales
-----------
Simulated time series for the sales of 3 different items.
Simulated data.
Shape of the dataset: (1097, 3)

Out[14]:

	item_1	item_2	item_3
date
2012-01-01	8.253175	21.047727	19.429739
2012-01-02	22.777826	26.578125	28.009863
2012-01-03	27.549099	31.751042	32.078922
2012-01-04	25.895533	24.567708	27.252276
2012-01-05	21.379238	18.191667	20.357737

In [15]:

# Custom class to create rolling skewness features (multi-series)
# ==============================================================================
from scipy.stats import skew


class RollingSkewnessMultiSeries():
    """
    Custom class to create rolling skewness features for multiple series.
    """

    def __init__(self, window_sizes, features_names: list = 'rolling_skewness'):
        
        if not isinstance(window_sizes, list):
            window_sizes = [window_sizes]
        self.window_sizes = window_sizes
        self.features_names = features_names

    def transform_batch(self, X: pd.Series) -> pd.DataFrame:
        
        rolling_obj = X.rolling(window=self.window_sizes[0], center=False, closed='left')
        rolling_skewness = rolling_obj.skew()
        rolling_skewness = pd.DataFrame({
                               self.features_names: rolling_skewness
                           }).dropna()

        return rolling_skewness

    def transform(self, X: np.ndarray) -> np.ndarray:
        
        n_series = X.shape[1]  # Series (levels) to be predicted (present in last_window)
        n_stats = 1  # Only skewness is calculated
        rolling_skewness = np.full(
            shape=(n_series, n_stats), fill_value=np.nan, dtype=float
        )
        for i in range(n_series):
            if len(X) > 0:
                rolling_skewness[i, :] = skew(X[:, i], bias=False)
            else:
                rolling_skewness[i, :] = np.nan        
        
        return rolling_skewness

# Custom class to create rolling skewness features (multi-series) # ============================================================================== from scipy.stats import skew class RollingSkewnessMultiSeries(): """ Custom class to create rolling skewness features for multiple series. """ def __init__(self, window_sizes, features_names: list = 'rolling_skewness'): if not isinstance(window_sizes, list): window_sizes = [window_sizes] self.window_sizes = window_sizes self.features_names = features_names def transform_batch(self, X: pd.Series) -> pd.DataFrame: rolling_obj = X.rolling(window=self.window_sizes[0], center=False, closed='left') rolling_skewness = rolling_obj.skew() rolling_skewness = pd.DataFrame({ self.features_names: rolling_skewness }).dropna() return rolling_skewness def transform(self, X: np.ndarray) -> np.ndarray: n_series = X.shape[1] # Series (levels) to be predicted (present in last_window) n_stats = 1 # Only skewness is calculated rolling_skewness = np.full( shape=(n_series, n_stats), fill_value=np.nan, dtype=float ) for i in range(n_series): if len(X) > 0: rolling_skewness[i, :] = skew(X[:, i], bias=False) else: rolling_skewness[i, :] = np.nan return rolling_skewness

In [16]:

# Forecaster with custom rolling features
# ==============================================================================
window_features = RollingSkewnessMultiSeries(window_sizes=3)
forecaster = ForecasterRecursiveMultiSeries(
                 regressor       = LGBMRegressor(random_state=123, verbose=-1),
                 lags            = 3,
                 window_features = window_features,
                 encoding        = 'ordinal'
             )

forecaster.fit(series=data_multiseries)
forecaster.predict(steps=5, levels=None)  # Predict all levels

# Forecaster with custom rolling features # ============================================================================== window_features = RollingSkewnessMultiSeries(window_sizes=3) forecaster = ForecasterRecursiveMultiSeries( regressor = LGBMRegressor(random_state=123, verbose=-1), lags = 3, window_features = window_features, encoding = 'ordinal' ) forecaster.fit(series=data_multiseries) forecaster.predict(steps=5, levels=None) # Predict all levels

Out[16]:

	item_1	item_2	item_3
2015-01-02	13.530901	19.709656	19.378026
2015-01-03	14.252429	18.608582	20.457204
2015-01-04	14.230324	19.536894	21.238047
2015-01-05	16.339439	19.174911	21.039226
2015-01-06	17.231169	18.821718	19.032785

In [17]:

# Transform output multi-series shape (n_levels, n_stats)
# ==============================================================================
window_features.transform(pd.DataFrame(forecaster.last_window_).to_numpy())

# Transform output multi-series shape (n_levels, n_stats) # ============================================================================== window_features.transform(pd.DataFrame(forecaster.last_window_).to_numpy())

Out[17]:

array([[-1.7304836 ],
       [ 1.69345527],
       [ 0.27695088]])