6

Feature Engineering for Time Series Forecasting

In the previous chapter, we started looking at machine learning (ML) as a tool to solve the problem of time series forecasting. We also discussed a few techniques, such as time delay embedding and temporal embedding, which cast time series forecasting problems as classical regression problems from the ML paradigm. In this chapter, we’ll look at those techniques in detail and go through them in a practical sense, using the dataset we have worked with throughout this book.

In this chapter, we will cover the following topics:

Technical requirements

You will need to set up the Anaconda environment, following the instructions in the Preface of the book, to get a working environment with all the libraries and datasets required for the code in this book. Any additional libraries needed will be installed while running the notebooks.

You will need to run the following notebooks before using the code in this chapter:

The code for this chapter can be found at https://github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-2E/tree/main/notebooks/Chapter06.

Understanding feature engineering

Feature engineering, as the name suggests, is the process of engineering features from data, mostly using domain knowledge, to make the learning process smoother and more efficient. In a typical ML setting, engineering good features is essential to get good performance from any ML model. Feature engineering is a highly subjective part of ML, where each problem at hand has a different path of solution—one that is handcrafted for that problem. Suppose you have a dataset of house prices and you have a feature, Year Built, which tells you the year the house was built. Now, to make the information better, we can create another feature, House Age, from the Year Built feature. This may give the model better information, and this is called feature engineering.

When we are casting a time series problem as a regression problem, there are a few standard techniques that we can apply. This is a key step in the process, as how well an ML model acquires an understanding of time is dependent on how well we engineer features to capture time. The baseline methods we covered in Chapter 4, Setting a Strong Baseline Forecast, are the methods that are created for the specific use case of time series forecasting, and because of that, the temporal aspect of the problem is built into those models. For instance, ARIMA doesn’t need any feature engineering to understand time because it is built into the model. However, a standard regression model has no explicit understanding of time, so we need to create good features to embed the temporal aspect of the problem.

In the previous chapter (Chapter 5, Time Series Forecasting as Regression), we talked about two main ways to encode time in the regression framework: time delay embedding and temporal embedding. Although we touched on these concepts at a high level, it is time to dig deeper and see them in action.

Notebook alert:

To follow along with the complete code, use the 01-Feature_Engineering.ipynb notebook in the Chapter06 folder.

We have already split the dataset that we were working on into train, validation, and test datasets. However, since we are generating features that are based on previous observations, operationally, it is better when we have the train, validation, and test datasets combined. It will be clearer why shortly, but for now, let’s take it on faith and move ahead. Now, let’s combine the two datasets:

# Reading the missing value imputed and train test split data
train_df = pd.read_parquet(preprocessed / "selected_blocks_train_missing_imputed.parquet")
val_df = pd.read_parquet(preprocessed / "selected_blocks_val_missing_imputed.parquet")
test_df = pd.read_parquet(preprocessed / "selected_blocks_test_missing_imputed.parquet")
#Adding train, validation and test tags to distinguish them before combining
train_df['type'] = "train"
val_df['type'] = "val"
test_df['type'] = "test"
full_df = pd.concat([train_df, val_df, test_df]).sort_values(["LCLid", "timestamp"])
del train_df, test_df, val_df

Now, we have a full_df, which combines the train, validation, and test datasets. Some of you may already have alarm bells ringing in your head at combining the train and test sets. What about data leakage? Let’s check it out.

Avoiding data leakage

Data leakage occurs when a model is trained with some information that would not be available at the time of prediction. Typically, this leads to high performance in the training set but very poor performance in unseen data. There are two types of data leakage:

When we work with time series forecasting problems, the biggest and most common mistake that we can make is target leakage. We will have to think hard about each of the features to ensure that we don’t use any data that will not be available during prediction. The following diagram will help us remember and internalize this concept:

Figure 6.1 – Usable and not-usable information to avoid data leakage

Figure 6.1: Usable and unusable information to avoid data leakage

To make this concept clearer and more relevant to the time series forecasting context, let’s look at an example. Let’s say we are forecasting sales for shampoo, and we are using sales for conditioner as a feature. We developed the model, trained it on the training data, and tested it on the validation data. The model does very well. The moment we start predicting for the future, we can see a problem. We don’t know what the sales for conditioner are in the future either. While this example is pretty straightforward, there will be times when this becomes not so obvious. And that is why we need to exercise a fair amount of caution while creating features and always evaluate the features through the lens of will this feature be available at the time of prediction?

Best practices:

There are many ways of identifying target leakage, apart from thinking hard about the features:

Although we generated forecasts earlier in this book, we never explicitly discussed forecast horizons. It is an important concept and essential for what we will discuss. Let’s take a bit of time to understand forecast horizons.

Setting a forecast horizon

A forecast horizon is the number of time steps into the future we want to forecast at any point in time. For instance, if we want to forecast the next 24 hours for the electricity consumption dataset that we have worked with, the forecast horizon becomes 48 (because the data is half-hourly). In Chapter 5, Time Series Forecasting as Regression, where we generated baselines, we just predicted the entire test data at once. In such cases, the forecast horizon becomes equal to the length of the test data.

We never had to worry about this until now because, in the classical statistical methods of forecasting, this decision is decoupled from modeling. If we train a model, we can use it to predict any future point without retraining. But with time series forecasting as regression, we have a constraint on the forecast horizon, and it has its roots in data leakage. This might be unclear to you now, so we’ll revisit this point after we have learned about the feature engineering techniques. For now, let’s only look at single-step-ahead forecasting. In the context of the dataset we are working with, this means that we will be answering the question, What is the energy consumption in the next half an hour? We will talk about multi-step forecasting and other mechanics of forecasting in Part 4, The Mechanics of Forecasting.

Now that we have set some ground rules, let’s start looking at the different feature engineering techniques. To follow along with the Jupyter notebook, head over to the Chapter06 folder and use the 01-Feature_Engineering.ipynb file.

Time delay embedding

The basic idea behind time delay embedding is to embed time in terms of recent observations. In Chapter 5, Time Series Forecasting as Regression, we discussed including previous observations of a time series as lags (Figure 5.6 under the subsection Time delay embedding).

However, there are a few more ways to capture recent and seasonal information using this concept.

Let’s take a look.

Lags or backshift

Let’s assume we have a time series with time steps, YL. Consider that we are at time T and that we have a time series where the length of history is L. So our time series will have yT as the latest observation in the time series, and then yT-1, yT-2, and so on as we move back in time. So lags, as explained in Chapter 5, Time Series Forecasting as Regression, are features that include the previous observations in the time series, as shown in the following diagram:

Figure 6.2 – Lag features

Figure 6.2: Lag features

We can create multiple lags by including observations that are a timesteps before (yT-a); we will call this Lag a. In the preceding diagram, we have shown Lag 1, Lag 2, and Lag 3. However, we can add any number of lags we like. Let’s learn how to do that now in code:

df["lag_1"]=df["column"].shift(1)

Remember when we combined the train and test datasets and I asked you to take it in good faith? It’s time to repay that faith. If we consider the lag operation (or any autoregressive feature), it relies on a continuous representation along the time axis. If we consider the test dataset, for the first few rows (or earliest dates), the lags would be missing because they are part of the training dataset. So by combining the two, we create a continuous representation along the time axis where standard functions in pandas, such as shift, can be utilized to create these features easily and efficiently.

It is as simple as that, but we need to perform the lag operation for each LCLid separately. We have included a helpful method in src.feature_engineering.autoregressive_features called add_lags that adds all the lags you want for each LCLid quickly and efficiently. Let’s see how we can use that.

We are going to import the method and use a few of its parameters to configure the lag operation the way we want:

from src.feature_engineering.autoregressive_features import add_lags
# Creating first 5 lags and then same 5 lags but from previous day and previous week to capture seasonality
lags = (
    (np.arange(5) + 1).tolist()
    + (np.arange(5) + 46).tolist()
    + (np.arange(5) + (48 * 7) - 2).tolist()
)
full_df, added_features = add_lags(
    full_df, lags=lags, column="energy_consumption", ts_id="LCLid", use_32_bit=True
)

Now, let’s look at the parameters that we used in the previous code snippet:

This method returns the DataFrame with the lags added, as well as a list with the column names of the newly added features.

Rolling window aggregations

With lags, we connect the present points to single points in the past, but with rolling window features, we connect the present with an aggregate statistic of a window from the past. Instead of looking at the observation from previous time steps, we would look at an average of the observations from the last three timesteps. Take a look at the following diagram to understand this better:

Figure 6.3 – Rolling window aggregation features

Figure 6.3: Rolling window aggregation features

We can calculate rolling statistics with different windows, and each of them will capture slightly different aspects of the history. In the preceding diagram, we can see an example of a window of three and a window of four. When we are at timestep T, a rolling window of three would have yT – 3, yT – 2, yT – 1 as the vector of past observations. Once we have these, we can apply any aggregation functions, such as the mean, standard deviation, min, max, and so on. Once we have a scalar value after the aggregation function, we can include that as a feature for timestep t.

We do not include yT in the vector of past observations because that leads to data leakage.

Let’s see how we can do this with pandas:

# We shift by one to make sure there is no data leakage
df["rolling_3_mean"] = df["column"].shift(1).rolling(3).mean()

Similar to the lags, we need to do this operation for each LCLid column separately. We have included a helpful method in src.feature_engineering.autoregressive_features called add_rolling_features that adds all the rolling features you want for each LCLid quickly and efficiently. Let’s see how we can use that.

We are going to import this method and use a few of its parameters to configure the rolling operation the way we want:

from src.feature_engineering.autoregressive_features import add_rolling_features
full_df, added_features = add_rolling_features(
    full_df,
    rolls=[3, 6, 12, 48],
    column="energy_consumption",
    agg_funcs=["mean", "std"],
    ts_id="LCLid",
    use_32_bit=True,
)

Now, let’s look at the parameters that we used in the previous code snippet:

This method returns the DataFrame with the rolling features added, as well as a list with the column names of the newly added features.

Seasonal rolling window aggregations

Seasonal rolling window aggregations are very similar to rolling window aggregations, but instead of taking past n consecutive observations in the window, they take a seasonal window, skipping a constant number of timesteps between each item in a window. The following diagram will make this clearer:

Figure 6.4 – Seasonal rolling window aggregations

Figure 6.4: Seasonal rolling window aggregations

The key parameter here is the seasonality period, which is commonly referred to as M. This is the number of timesteps after which we expect the seasonality pattern to repeat. When we are at timestep T, a rolling window of three would have yT – 3, yT – 2, yT – 1, as the vector of past observations. But the seasonal rolling window would skip m timesteps between each item in the window. This means that the observations that are there in the seasonal rolling window would be yT – M, yT – 2M, yT – 3M. Also, as usual, once we have the window vector, we just need to apply the aggregation function to get a scalar value and include that as a feature.

We do not include yT as an element in the seasonal rolling window vector to avoid data leakage.

This is not an operation that you can do easily and efficiently with pandas. Some fancy NumPy indexing and Python loops should do the trick. We will use an implementation from github.com/jmoralez/window_ops/ that uses NumPy and Numba to make the operation fast and efficient.

Just like the features we saw earlier, we need to do this operation for each LCLid separately. We have included a helpful method in src.feature_engineering.autoregressive_features called add_seasonal_rolling_features that adds all the seasonal rolling features you want for each LCLid quickly and efficiently. Let’s see how we can use that.

We are going to import the method and use a few parameters of the method to configure the seasonal rolling operation the way we want:

from src.feature_engineering.autoregressive_features import add_seasonal_rolling_features
full_df, added_features = add_seasonal_rolling_features(
    full_df,
    rolls=[3],
    seasonal_periods=[48, 48 * 7],
    column="energy_consumption",
    agg_funcs=["mean", "std"],
    ts_id="LCLid",
    use_32_bit=True,
)

Now, let’s look at the parameters that we used in the previous code snippet:

As always, the method returns the DataFrame with seasonal rolling features and a list containing the column names of the newly added features.

Exponentially weighted moving average (EWMA)

With the rolling window mean operation, we calculated the average of the window, and it works synonymously with the moving average. EWMA is the slightly smarter cousin of the moving average. While the moving average considers a rolling window and considers each item in the window equally on the computed average, EWMA tries to do a weighted average on the window, and the weights decay at an exponential rate. There is a parameter, $\alpha$, that determines how fast the weights decay. Because of this, we can consider all the history available as a window and let the $\alpha$ parameter decide how much recency is included in EWMA. This can be written simply and recursively, as follows:

$$EWMA_T = \alpha \times y_T + (1 - \alpha) \times EWMA_{T-1}$$

Here, we can see that the larger the value of $\alpha$, the more the average is skewed toward recent values (see Figure 6.6 to get a visual impression of how the weights would be). If we expand the recursion, the weights of each term work out to be:

$$W_{T-k} = \alpha \times (1 - \alpha)^k$$

where k is the number of timesteps behind T. If we plot the weights, we can see them in an exponential decay; $\alpha$ determines how fast the decay happens. Another way to think about $\alpha$ is in terms of span. Span is the number of periods at which the decayed weights approach zero (not in a strictly mathematical way but intuitively). $\alpha$ and span are related through this equation:

$$\alpha = \frac{2}{1 + span}$$

This will become clearer in the following diagram, where we have plotted how the weights decay for different values of $\alpha$:

Figure 6.5 – Exponential weight decay for different values of

Figure 6.5: Exponential weight decay for different values of $\alpha$

Here, we can see that the weight becomes small by the time we reach the span.

Intuitively, we can think of EWMA as an average of the entire history of the time series, but with parameters such as $\alpha$ and span, we can make different periods of history more representative of the average. If we define a 60-period span, we can think that the last 60 time periods are what majorly drive the average. So making EWMAs with different spans or $\alpha$ s gives us representative features that capture different periods of history.

The overall process is depicted in the following diagram:

Figure 6.6 – EWMA features

Figure 6.6: EWMA features

Now, let’s see how we can do this in pandas:

df["ewma"]=df['column'].shift(1).ewm(alpha=0.5).mean()

Like the other features we discussed earlier, EWMA also needs to be done for each LCLid separately. We have included a helpful method in src.feature_engineering.autoregressive_features called add_ewma that adds all the EWMA features you want for each LCLid quickly and efficiently. Let’s see how we can use that.

We are going to import the method and use a few parameters of the method to configure EWMA the way we want to:

from src.feature_engineering.autoregressive_features import add_ewma
full_df, added_features = add_ewma(
    full_df,
    spans=[48 * 60, 48 * 7, 48],
    column="energy_consumption",
    ts_id="LCLid",
    use_32_bit=True,
)

Now, let’s look at the parameters that we used in the previous code snippet:

As always, the method returns the DataFrame containing the EWMA features, as well as a list with the column names of the newly added features.

These are a few standard ways of including time delay embedding in your ML model, but you are not restricted to just these. As always, feature engineering is a space that is not bound by rules, and we can get as creative as we want and inject domain knowledge into the model. Apart from the features we have seen, we can include the difference in lag as custom lags that inject domain knowledge, and so on. In most practical cases, we end up using more than one way of time delay embedding into our models. The lag feature is the most basic and essential in most cases, but we do end up encoding more information with seasonal lags, rolling features, and so on. As with everything in ML, there is no silver bullet. Each dataset has its own intricacies, which makes feature engineering very important and different for each case.

Now, let’s look at the other class of features we can add via temporal embedding.

Temporal embedding

In Chapter 5, Time Series Forecasting as Regression, we briefly discussed temporal embedding as a process where we try to embed time into features that an ML model can leverage. If we think about time for a second, we can see that two aspects of time are important to us in the context of time series forecasting—the passage of time and the periodicity of time.

There are a few features that we can add to help us capture these aspects in an ML model:

Let’s look at each of them.

Calendar features

The first set of features that we can extract are features based on calendars. Although the strict definition of time series is a set of observations taken sequentially in time, more often than not, we will have the timestamps of these collected observations alongside the time series. We can utilize these timestamps and extract calendar features, such as the month, quarter, day of the year, hour, minutes, and so on. These features capture the periodicity of time and help an ML model capture seasonality well. Only the calendar features that are temporally higher than the frequency of the time series make sense. For instance, an hour feature in a time series with a weekly frequency doesn’t make sense, but a month feature and week feature make sense. We can utilize in-built datetime functionalities in pandas to create these features and treat them as categorical features in the model.

Time elapsed

This is another feature that captures the passage of time in an ML model. This feature increases monotonically as time increases, giving the ML model a sense of the passage of time. There are many ways to create this feature, but one of the easiest and most efficient ways is to use the integer representation of dates in NumPy:

df['time_elapsed'] = df['timestamp'].values.astype(np.int64)/(10**9)

We have included a helpful method in src.feature_engineering.temporal_features called add_temporal_features that adds all relevant temporal features automatically. Let’s see how we can use it.

We are going to import the method and use a few parameters of this method to configure and create the temporal features:

full_df, added_features = add_temporal_features(
    full_df,
    field_name="timestamp",
    frequency="30min",
    add_elapsed=True,
    drop=False,
    use_32_bit=True,
)

Now, let’s look at the parameters that we used in the previous code snippet:

Just like the previous methods we discussed, this also returns the new DataFrame with the temporal features added and a list containing the column names of the newly added features.

Fourier terms

Previously, we extracted a few calendar features such as the month, year, and so on, and we discussed using them as categorical variables in the ML model. Another way we can represent the same information, but on a continuous scale, is by using Fourier terms. We discussed the Fourier series in Chapter 3, Analyzing and Visualizing Time Series Data. Just to reiterate, the sine-cosine form of the Fourier series is as follows:

$$S_{N(x)} = \frac{a_0}{2} + \sum_{n=1}^{N} \left( a_n \cdot \cos\left(\frac{2\pi \cdot n \cdot x}{P}\right) + b_n \cdot \sin\left(\frac{2\pi \cdot n \cdot x}{P}\right) \right)$$

Here, SN is the N-term approximation of the signal, S. Theoretically, when N is infinite, the resulting approximation is equal to the original signal. P is the maximum length of the cycle, an and bn are the coefficients of the cosine and sine term, respectively, of the nth term in the expansion, and a0 is the intercept.

We can create these cosine and sine functions as features to represent the seasonal cycle. If we encode the month, we know that it goes from 1 to 12 and then repeats itself. So P, in this case, will be 12, and x will be 1, 2, …12. Therefore, for each x, we can calculate the cosine and sine terms and add them as features to the ML model. Intuitively, we can think that the model will infer the coefficients based on the data and, thus, help the model predict the time series easier.

The following diagram shows the difference in representations between the month on an ordinal scale and as a Fourier series:

Figure 6.7 – Month as an ordinal step function (top) versus Fourier terms (bottom)

Figure 6.7: Month as an ordinal step function (top) versus Fourier terms (bottom)

The preceding diagram shows just a single Fourier term; we can add multiple Fourier terms to help capture complex seasonality.

We cannot say that continuous representation of seasonality is better than categorical because it depends on the type of model you use and the dataset. This is something we will have to find out empirically.

To make the process of adding Fourier features easy, we have made some easy-to-use methods available in src.feature_engineering.temporal_features, in a file called bulk_add_fourier_features that adds Fourier features for all the calendar features we want automatically. Let’s see how we can use that.

We are going to import the method and use a few of its parameters to configure and create the Fourier series-based features:

full_df, added_features = bulk_add_fourier_features(
    full_df,
    ["timestamp_Month", "timestamp_Hour", "timestamp_Minute"],
    max_values=[12, 24, 60],
    n_fourier_terms=5,
    use_32_bit=True,
)

Now, let’s look at the parameters that we used in the previous code snippet:

Just like the previous methods we’ve discussed, this also returns a new DataFrame with the Fourier features added, as well as a list with column names of the newly added features.

After executing the 01-Feature_Engineering.ipynb notebook in Chapter06, we will have the following feature-engineered files written to disk:

In this section, we looked at a few popular and effective ways to generate features for time series. But there are many more, and depending on your problem and the domain, many of them will be relevant.

Additional information:

The world of feature engineering is vast, and there are a few open-source libraries that make exploring that space easier. A few of them are https://github.com/Nixtla/tsfeatures, https://tsfresh.readthedocs.io/en/latest/, and https://github.com/DynamicsAndNeuralSystems/catch22. A preprint by Ben D. Fulcher titled Feature-based time-series analysis at https://arxiv.org/abs/1709.08055 also gives a nice summary of the space.

A newer library called functime (https://github.com/functime-org/functime) also provides fast feature engineering routines, written in Polars, and is worth checking out. A lot of the feature engineering discussed in the book can be made even faster using functime and Polars.

Summary

After a brief overview of the ML for time series forecasting paradigm in the previous chapter, in this chapter, we looked at this practically and saw how we can prepare the dataset with the required features to start using these models. We reviewed a few time series-specific feature engineering techniques, such as lags, rolling, and seasonal features. All the techniques we learned in this chapter are tools with which we can quickly iterate through experiments to find out what works for our dataset. However, we only talked about feature engineering, which affects one side of the standard regression equation (y = mX + c). The other side, which is the target (y) we predict, is also equally important. In the next chapter, we’ll look at a few concepts such as stationarity and some transformations that affect the target.