8

Forecasting Time Series with Machine Learning Models

In the previous chapter, we started looking at machine learning as a tool to solve the problem of time series forecasting. We talked about a few techniques such as time delay embedding and temporal embedding, both of which cast a time series forecasting problem as a classical regression problem from the machine learning paradigm. In this chapter, we’ll look at these techniques in detail and go through them in a practical sense using the London Smart Meters dataset we have been working 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 library 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-/tree/main/notebooks/Chapter08.

Training and predicting with machine learning models

In Chapter 5, Time Series Forecasting as Regression, we talked about a schematic for supervised machine learning (Figure 5.2). In the schematic, we mentioned that the purpose of a supervised learning problem is to come up with a function, $\hat{y} = h(X, \phi)$, where $\hat{y}$ is the predicted value, X is the set of features as the input, $\phi$ is the model parameters, and h is the approximation of the ideal function. In this section, we are going to talk about h in more detail and see how we can use different machine learning models to estimate it.

h is any function that approximates the ideal function, but it can be thought of as an element of all possible functions from a family of functions. More formally, we can say the following:

$$\hat{y} = h(X, \phi), \quad \text{where } h \in H$$

Here, H is a family of functions that we also call a model. For instance, linear regression is a type of model or a family of functions. For each value of the coefficients, the linear regression model gives you a different function and H becomes the set of all possible functions a linear regression model can produce.

There are many families of functions, or models, available. For a more complete understanding of the space, we will need to refer to other machine learning resources. The Further reading section contains a few resources that may help you start the journey. As for the scope of this book, we narrowly define it as the application of machine learning models for forecasting, rather than machine learning in general. And although we can use any regression model, we will only review a few popular and useful ones for time series forecasting and see them in action. We leave it to you to strike out on your own and explore the other algorithms to become familiar with them as well. But before we look at the different models, we need to generate a few baselines again.

Generating single-step forecast baselines

We reviewed and generated a few baseline models back in Chapter 4, Setting a Strong Baseline Forecast. But there is a small issue – the prediction horizon. In Chapter 6, Feature Engineering for Time Series Forecasting, we talked about how the machine learning model can only predict one target at a time and that we are sticking with a single-step forecast. The baselines we generated earlier were not single-step, but multi-step. Generating a single-step forecast for baseline algorithms such as ARIMA or ETS requires us to fit on history, predict one step ahead, and then fit again using one more day. Predicting in such an iterative fashion for our test or validation period requires us to do this iteration ~1,440 times (48 data points a day for 30 days) and repeat this for all the households in our selected dataset (150, in our case). This would take quite a long time to compute.

We have chosen the naïve method and seasonal naïve (Chapter 4, Setting a Strong Baseline Forecast), which can be implemented as native pandas methods, as two baseline methods to generate single-step forecasts.

Naïve forecasts perform unreasonably well for single-step-ahead forecasts and can be considered a strong baseline. In the Chapter08 folder, there is a notebook named 00-Single_Step_Backtesting_Baselines.ipynb that generates these baselines and saves them to disk. Let’s run the notebook now. The notebook generates the baselines for both the validation and test datasets and saves the predictions, metrics, and aggregate metrics to disk. The aggregate metrics for the test period are as follows:

Figure 8.1 – Aggregate metrics for a single-step baseline

Figure 8.1: Aggregate metrics for a single-step baseline

To make training and evaluating these models easier, we have used a standard structure throughout. Let’s quickly review that structure as well so that you can follow along with the notebooks closely.

Standardized code to train and evaluate machine learning models

There are two main ingredients while training a machine learning model – data and the model itself. Therefore, to standardize the pipeline, we defined three configuration classes (FeatureConfig, MissingValueConfig, and ModelConfig) and another wrapper class (MLForecast) over scikit-learn-style estimators (.fit - .predict) to make the process smooth. Let’s look at each of them.

Notebook alert:

To follow along with the code, use the 01-Forecasting_with_ML.ipynb notebook in the Chapter08 folder and the code in the src folder.

FeatureConfig

FeatureConfig is a Python dataclass that defines a few key attributes and functions that are necessary while processing the data. For instance, continuous, categorical, and Boolean columns need separate kinds of preprocessing before being fed into the machine learning model. Let’s see what FeatureConfig holds:

In addition to a bit of validation on the inputs, there is also a helpful method called get_X_y in the class, with the following parameters:

The function returns a tuple of (features, target, original_target).

All we need to do is initialize the class, like any other class, with the feature names separated into the parameters of the class. The entire code that contains all the features is available in the accompanying notebook.

After setting the FeatureConfig data class, we can pass any DataFrame with the features defined to the get_X_y function to get the features, target, and original target:

train_features, train_target, train_original_target = feat_config.get_X_y(
    sample_train_df, categorical=False, exogenous=False
)

As you can see, we are not using categorical features or exogenous features here, as I want to focus on the core algorithms and show how they can be drop-in replacements for other classical time series models we saw earlier. We will talk about how to handle categorical features in Chapter 15, Strategies for Global Deep Learning Forecasting Models.

MissingValueConfig

Another key setting is how to deal with missing values. We saw a few ways to fill in missing values from a time series context in Chapter 3, Analyzing and Visualizing Time Series Data, and we have already filled in missing values and prepared our datasets. But a few missing values will be created in the feature engineering required to convert a time series into a regression problem. For instance, when creating lag features, the earliest date in the dataset will not have enough data to create a lag and will be left empty.

Best practice:

Although filling with zero or mean is the default or go-to method for the majority of the data scientist community, we should always make an effort to fill in missing values as intelligently as possible. In terms of lag features, filling with zero can distort the feature. Instead of filling with zero, a backward fill (using the earliest value in the column to fill backward) might be a much better fit.

Some machine learning models handle empty or NaN features naturally, while for other machine learning models, we will need to deal with such missing values before training. It’s helpful if we can define a config in which we set for a few columns where we expect NaN information on how to fill those. MissingValueConfig is a Python dataclass that does just that. Let’s see what it holds:

The order in which the missing values are filled is bfill_columns, then ffill_columns, and then zero_fill_columns. As the default strategy, the data class uses the column mean to fill in missing values so that even if you have not defined any strategy for a column, the missing value will be filled in using a column mean. There is a method called impute_missing_values that takes in the DataFrame and fills the empty cells with a value according to the specified strategy.

ModelConfig

ModelConfig is a Python dataclass that holds a few details regarding the modeling process, such as whether to normalize the data, whether to fill in missing values, and so on. Let’s take a detailed look at what it holds:

Let’s see how we can define the ModelConfig data class:

model_config = ModelConfig(
    model=LinearRegression(),
    name="Linear Regression",
    normalize=True,
    fill_missing=True,
)

This has just one method, clone, that clones the estimator, along with the config, into a new instance.

MLForecast

Last but not least, we have the wrapper class around a scikit-learn-style model. It uses the different configurations we have discussed to encapsulate the training and prediction functions. Let’s see what parameters are available when initializing the model:

MLForecast has a few functions that can help us manage the life cycle of a model, once initialized. Let’s take a look.

The fit function

The fit function is similar in purpose to the scikit-learn fit function but does a little extra by handling the standardization, categorical encoding, and target transformations using the information in the three configs. The parameters of the function are as follows:

The predict function

The predict function handles inferencing. It wraps around the predict function of the scikit-learn estimator, but like fit, it does a few other things, such as standardization, categorical encoding, and reversing the target transformation. There is only one parameter for this function:

The feature_importance function

The feature_importance function retrieves the feature importance from the model, if available. For linear models, it extracts the coefficients, while for tree-based models, it extracts the built-in importance and returns it in a sorted DataFrame.

Helper functions for evaluating models

While the other functions we saw earlier deal with core training and predicting, we also want to evaluate the model, plot the results, and so on. We have also defined these functions in the notebooks or in the code base. The below function is to evaluate the models in the notebook:

def evaluate_model(
    model_config,
    feature_config,
    missing_config,
    train_features,
    train_target,
    test_features,
    test_target,
):
    ml_model = MLForecast(
        model_config=model_config,
        feature_config=feat_config,
        missing_config=missing_value_config,
    )
    ml_model.fit(train_features, train_target)
    y_pred = ml_model.predict(test_features)
    feat_df = ml_model.feature_importance()
    metrics = calculate_metrics(test_target, y_pred, model_config.name, train_target)
    return y_pred, metrics, feat_df

This provides us with a standard way of evaluating all the different models, as well as automating the process at scale. We also have a function for calculating the metrics, calculate_metrics, defined in src/forecasting/ml_forecasting.py.

The standard implementation that we have provided with this book is in no way a one-size-fits-all approach, but rather something that works best with the flow and dataset of this book. Please do not consider it as a robust library, but rather a good starting point and guide to help you develop your own code.

Now that we have the baselines and a standard way to apply different models, let’s get back to what the different models are. For the discussion ahead, let’s keep time out of our minds because we have converted a time series forecasting problem into a regression problem and factored in time as a feature of the problem (the lags and rolling features).

Linear regression

Linear regression is a family of functions that takes the following form:

$$\hat{y} = \beta_0 + \sum_{i=1}^{k} X_i \beta_i$$

Here, k is the number of features in the model and $\beta$ are the parameters of the model. There is a $\beta$ for each feature, as well as a $\beta_0$, which we call the intercept, which is estimated from data. Essentially, the output is a linear combination of the feature vectors, Xi. As the name suggests, this is a linear function.

The model parameters can be estimated from data,D(Xi, yi), using an optimization method and loss, but the most popular method of estimation is using ordinary least squares (OLS). Here, we find the model parameters, $\beta$, which minimizes the residual sum of squares (mean squared error (MSE)):

$$RSS = \sum_{i=1}^{N} (y_i - \hat{y}_i)^2$$

The loss function here is very intuitive. We are essentially minimizing the distance between the training samples and our predicted points. The square term acts as a technique that does not cancel out positive and negative errors. Apart from the intuitiveness of the loss, another reason why this is widely chosen is that an analytical solution exists for least squares and because of that, we don’t need to resort to more compute-intensive optimization techniques such as gradient descent.

Linear regression has one foot firmly planted in statistics and with the right assumptions, it can be a powerful tool. Commonly, five assumptions are associated with linear regression, as follows:

But unless you are concerned about using linear regression to come up with prediction intervals (a band in which the prediction would lie with some probability), we can disregard all but the first assumption to some extent.

The linearity assumption (the first assumption) is relevant because if the variables are not linearly related, it will result in an underfit and thus poor performance. We can get around this problem to some extent by projecting the inputs into a higher dimensional space. Theoretically, we can project a non-linear problem into a higher-dimensional space, where the problem is linear. For instance, let’s consider a non-linear function, $y = 3x_1^2 + 2x_2^2 + 6x_1 x_2$. If we run linear regression in the input space of $x_1$ and $x_2$, we know the resulting model will be highly underfitting. But if we project the input space from $x_1$ and $x_2$ to $x_1^2$, $x_2^2$, and $xy$ by using a polynomial transform, the function for y becomes a perfect linear fit.

The multi-collinearity assumption (the final assumption) is partly relevant to the fit of the linear function because when we have highly correlated independent variables, the estimated coefficients are highly unstable and difficult to interpret. The fitted function would still be working well, but because we have multi-collinearity, even small changes in the inputs would make the coefficients change magnitude and sign. It is a best practice to check for multi-collinearity if you are using a pure linear regression. This is typically a problem in time series because the features we have extracted, such as the lag and rolling features, may be correlated with each other. Therefore, we will have to be careful while using and interpreting linear regression on time series data.

Now, let’s see how we can use linear regression and evaluate the fit of a sample household from our validation dataset:

from sklearn.linear_model import LinearRegression
model_config = ModelConfig(
    model=LinearRegression(),
    name="Linear Regression",
    # LinearRegression is sensitive to normalized data
    normalize=True,
    # LinearRegression cannot handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

The single-step forecast looks good and is already better than the naïve forecast (MAE = 0.173):

Figure 8.2 – Linear regression forecast

Figure 8.2: Linear regression forecast

The coefficients of the model, $\beta$ (which can be accessed using the coef_ attribute of a trained scikit-learn model), show how much influence each feature has on the output. So, extracting and plotting them gives us our first level of visibility into the model. Let’s take a look at the coefficients of the model:

Figure 8.3 – Feature importance of linear regression (top 15)

Figure 8.3: Feature importance of linear regression (top 15)

If we look at the Y-axis in the feature importance chart, we can see it is in billions as the coefficient for a couple of features is in orders of magnitude in billions. We can also see that those features are Fourier series-based features, which are correlated with each other. Even though we have a lot of coefficients that are in billions, we can find them on both sides of zero, so they will essentially cancel out each other in the function. This is the problem with multi-collinearity that we talked about earlier. We can go about removing multi-collinear features and then perform some sort of feature selection (forward selection or backward elimination) to make the linear model even better.

But instead of doing that, let’s look at a few modifications we can make to the linear model that are a bit more robust to multi-collinearity and feature selection.

Regularized linear regression

We briefly talked about regularization in Chapter 5, Time Series Forecasting as Regression, and mentioned that regularization, in the general sense, is any kind of constraint we place on the learning process to reduce the complexity of the learned function. One of the ways linear models can become more complex is by having a high magnitude of coefficients. For instance, in the linear fit, we have a coefficient of 20 billion. Any small change in that feature is going to cause a huge fluctuation in the resulting prediction. Intuitively, if we have a large coefficient, the function becomes more flexible and complex. One way we can fix this is to apply regularization in the form of weight decay. Weight decay is when we add a term that penalizes the magnitude of the coefficients to the loss function. The loss function, the residual sum of squares, now becomes as follows:

$$RSS = \sum_{i=1}^{N} (y_i - \hat{y}_i)^2 + \lambda \mathcal{W}$$

Here, W is the weight decay and $\lambda \lt 0$ is the strength of regularization.

W is typically the norm of the weight matrix. In linear algebra, the norm of a matrix is a measure of how large its elements are. There are many norms for a matrix, but the two most common norms that are used for regularization are the L1 and L2 norms. When we use the L1 norm to regularize linear regression, we call it lasso regression, while when we use the L2 norm, we call it ridge regression. When we apply weight decay regularization, we are forcing the coefficients to be lower, which means that it also acts as an internal feature selection because the features that don’t add a lot of value will get very low or zero (depending on the type of regularization) coefficients, which means they contribute little to nothing in the resulting function.

The L1 norm is defined as the sum of the absolute values of the matrix. For weight decay regularization, the L1 norm would be as follows:

$$W = \sum_{i=1}^{k} |\beta_i|$$

The L2 norm is defined as the sum of squared values of a matrix. For weight decay regularization, the L2 norm would be as follows:

$$W = \sum_{i=1}^{k} |\beta_i^2|$$

By adding this term to the loss function of linear regression, we are forcing the coefficients to be small because while the optimizer is reducing the RSS, it is also incentivized to reduce W.

Another way we can think about regularization is in terms of linear algebra and geometry.

The following section discusses the geometric intuition of regularization. Although it would make your understanding of regularization more solid, it is not essential to be able to follow the rest of this book. So, feel free to skip the next section and just read the Key point callout if you are pressed for time or if you want to come back to it later when you have time.

Regularization–a geometric perspective

If we look at the L1 and L2 norms from a slightly different perspective, we will see that they are measures of distance.

Let B be the vector of all the coefficients, $\beta$, in linear regression. A vector is an array of numbers, but geometrically, it is also an arrow from the origin to a point in the n-dimensional coordinate space. Now, the L2 norm is nothing but the Euclidean distance from the origin on that point in space defined by the vector, B. The L1 norm is the Manhattan distance or taxicab distance from the origin on that point in space defined by the vector, B. Let’s see this in a diagram:

Figure 8.4 – Euclidean versus Manhattan distance

Figure 8.4: Euclidean versus Manhattan distance

Euclidean distance is the length of the direct path from the origin to the point. But if we can only move parallel to the two axes, we will have to travel the distance of $\beta_0$ along the one axis first, and then a distance of $\beta_1$ along the other. This is the Manhattan distance.

Let’s say we are in a city (for example, Manhattan) where the buildings are laid out in square blocks and the straight streets intersect at right angles, and we want to travel from point A to point B. Euclidean distance is the direct distance from point A to point B, which in the real sense is only possible if we parkour along the tops of the buildings. On the other hand, the Manhattan distance is the actual distance a taxicab would take while traveling along the right-angled roads from point A to point B.

To develop further geometrical intuition about the L1 and L2 norms, let’s do one thought experiment. If we move the point, $(\beta_0, \beta_1)$, in the 2D space while keeping the Euclidean distance or the L2 norm the same, we will end up with a circle with its center at the origin. This becomes a sphere in 3D and a hypersphere in n-D. If we trace out the same but keep the L1 norm the same, we will end up with a diamond with its center at the origin. This would become a cube in 3D and a hypercube in n-D.

Now, when we are optimizing for the weights, in addition to the main objective of reducing the loss function, we are also encouraging the coefficients to stay within a defined distance (norm) from the origin. Geometrically, this means that we are asking the optimization to find a vector, $\beta$, that minimizes the loss function and stays within the geometric shape (circle or square) defined by the norm. We can see this in the following diagram:

Figure 8.5 – Regularization with the L1 Norm (lasso regression) versus the L2 Norm (ridge regression)

Figure 8.5: Regularization with the L1 norm (lasso regression) versus the L2 norm (ridge regression)

The concentric circles in the diagram are the contours of the loss function, with the innermost being the lowest. As we move outward, the loss increases. So, instead of selecting a $\beta_{\text{opt}}$, regularized regression will select a $\beta$ that intersects with the norm geometry.

This geometric interpretation also makes understanding another key difference between ridge and lasso regression easier. Lasso regression, because of the L1 norm, produces a sparse solution. Earlier, we mentioned that weight decay regularization does implicit feature selection. But depending on whether you are applying the L1 or L2 norm, the kind of implicit feature selection differs.

Key point:

For an L2 norm, the coefficients of less relevant features are pushed close to zero, but not exactly zero. The feature will still play a role in the final function, but its influence will be minuscule. The L1 norm, on the other hand, pushes the coefficients of such features completely to zero, resulting in a sparse solution. Therefore, L1 regularization promotes sparsity and feature selection, whereas L2 regularization reduces model complexity by shrinking the coefficients toward zero without necessarily eliminating any.

This can be understood better using the geometrical interpretation of regularization. In optimization, the interesting points are usually found in the extrema or corners of a shape. There are no corners in a circle, so an L2 norm is created; the minima can lie anywhere on the edge of the circle. But for the diamond, we have four corners, and the minima would lie in those corners. So, with the L2 norm, the solution can move very close to zero, but not necessarily zero. However, with the L1 norm, the solution would be on the corners, where the coefficient can be pushed to an absolute zero.

Now, let’s see how we can use ridge regression and evaluate the fit on a sample household from our validation dataset:

from sklearn.linear_model import RidgeCV
model_config = ModelConfig(
    model=RidgeCV(),
    name="Ridge Regression",
    # RidgeCV is sensitive to normalized data
    normalize=True,
    # RidgeCV does not handle missing values
    fill_missing=True
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

Let’s look at the single-step-ahead forecast from RidgeCV. It looks very similar to linear regression. Even the MAE is the same for this household:

Figure 8.6 – Ridge regression forecast

Figure 8.6: Ridge regression forecast

But it is interesting to look at the coefficients with the L2 regularized model. Let’s take a look at the coefficients of the model:

Figure 8.7 – Feature importance of ridge regression (top 15)

Figure 8.7: Feature importance of ridge regression (top 15)

Now, the Y-axis looks reasonable and small. The coefficients for the multi-collinear features have shrunk to a more reasonable level. Features such as the lag features, which should ideally be highly influential, have gained the top spots. As you may recall, in the linear regression (Figure 8.3), these features were dwarfed by the huge coefficients on the Fourier features. We have just plotted the top 15 features here, but if you look at the entire list, you will see that there will be a lot of features for which the coefficients are close to zero.

Now, let’s try lasso regression on the sample household:

from sklearn.linear_model import LassoCV
model_config = ModelConfig(
    model=LassoCV(),
    name="Lasso Regression",
    # LassoCV is sensitive to normalized data
    normalize=True,
    # LassoCV does not handle missing values
    fill_missing=True
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

Let’s look at the single-step-ahead forecast from LassoCV. Like ridge regression, there is hardly any visual difference from linear regression:

Figure 8.8 – Lasso regression forecast

Figure 8.8: Lasso regression forecast

Let’s look at the coefficients of the model:

Figure 8.9 – Feature importance of lasso regression (top 15)

Figure 8.9: Feature importance of lasso regression (top 15)

The coefficients are very similar to ridge regression, but if you look at the full list of coefficients (in the notebook), you will see that there are a lot of features where the coefficients will be zero.

Even with the same MAE, MSE, and so on, ridge or lasso regression is preferred to linear regression because of the additional stability and robustness that comes with regularized regression, especially for forecasting, where multi-collinearity is almost always present. But we need to keep in mind that all the linear regression models are still only capturing linear relationships. If the dataset has a non-linear relationship, the resulting fit from linear regression won’t be as good and, sometimes, will be terrible.

Now, let’s switch tracks and look at another class of models – decision trees.

Decision trees

Decision trees are another family of functions that is much more expressive than a linear function. Decision trees split the feature space into different sub-spaces and fit a very simple model (such as an average) to each. Let’s understand how this partitioning works with an example. Let’s consider a regression problem for predicting Y with just one feature, X, as shown in the following diagram:

Figure 8.10 – The feature space partitioned by a decision tree

Figure 8.10: The feature space partitioned by a decision tree

Right away, we can see that fitting a linear function would result in an underfit. But what decision trees do is split the feature space (here, it is just X) into different regions where the target, Y, is similar and then fit a simple function such as an average (because it is a regression problem). In this case, the decision tree has split the feature space into partitions – A, B, and C. Now, for any X that falls into partition A, the prediction function will return the average of all the points in partition A.

These partitions are formed by creating a decision tree using data. Intuitively, a decision tree creates a set of if-else conditions and tries to arrive at the best way to partition the feature space to maximize the homogeneity of the target variable within the partition. One helpful way to understand what a decision tree does is to think of data points as beads flowing down a tree, taking a path based on its features, and ending up in a final resting place. Before we talk about how to create a decision tree from data, let’s take a look at its components and understand the terminology surrounding it:

Figure 8.11 – Anatomy of a decision tree

Figure 8.11: Anatomy of a decision tree

There are two types of nodes in a decision tree—a decision node and a leaf node. A decision node is the if-else statement we mentioned previously. This node will have a condition based on whether the data points that flow down the tree take the left or right branch. The decision node that sits right at the top has a special name—the root node. Finally, the process of dividing the data points based on a condition and directing it to the right or left branch is called splitting. Leaf nodes are nodes that don’t have any other branches below them. These are the final resting points in the beads flowing down a tree analogy. These are the partitions we discussed earlier in this section.

Formally, we can define the function that’s been generated by a decision tree that has M partitions, P1, P2, …, PM, as follows:

$$\hat{y} = \sum_{m=1}^{M} c_m I(x \in P_m)$$

Here, x is the input, cm is the constant response for the region, Pm, and I is a function that is 1 if $(x \in P_m)$; otherwise, it’s 0.

For regression trees, we usually adopt the squared loss as the loss function. In that case, cm is usually set as the average of all y, where the corresponding x falls in the Pm partition.

Now that we know how a decision tree functions, the only thing left to understand is how to decide which feature to split on and where to split the feature.

Many algorithms have been proposed over the years on how to create a decision tree from data such as ID3, C4.5, CART, and so on. Using Classification and Regression Trees (CART) is one of the most popular methods out of the lot and supports regression as well. Therefore, we will just stick to CART in this book. Classification Trees are used when the target variable is categorical (e.g., predicting a class label). Regression Trees are used when the target variable is continuous (e.g., predicting a numerical value).

The most optimal set of binary partitions that minimizes the sum of squares globally is generally intractable. So, we adopt a greedy algorithm to create the decision tree. Greedy optimization is a heuristic that builds up a solution stage by stage, selecting a local optimum at each stage. Therefore, instead of finding the best feature splits globally, we will create the decision tree, decision node by decision node, where we choose the most optimal feature split at each stage. For a regression tree, we choose a split feature, f, and split point, s, so that it creates two partitions, P1 and P2, that minimize as follows:

$$\sum_{x_i \in P_1} (y_i - c_1)^2 + \sum_{x_i \in P_2} (y_i - c_2)^2$$

Here, c1 and c2 are the averages of all y, where the corresponding x falls in between P1 and P2.

Therefore, by using this criterion, we can keep splitting the regions further and further. With each level we split, we increase the depth of the tree by one. At some point, we will start overfitting the dataset. But if we don’t do enough splits, we might be underfitting the data as well. One strategy is to stop creating further splits when we reach a predetermined depth. In the scikit-learn implementation of DecisionTreeRegressor, this corresponds to the max_depth parameter. This is a hyperparameter that needs to be estimated using a validation dataset. There are other strategies to stop the splits, such as setting a minimum number of samples required to split (min_samples_split), or a minimum decrease in cost to carry out a split (min_impurity_decrease). For a complete list of parameters in DecisionTreeRegressor, please refer to the documentation at https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html.

Now, let’s see how we can use a decision tree and evaluate the fit on a sample household from our validation dataset:

from sklearn.tree import DecisionTreeRegressor
model_config = ModelConfig(
    model=DecisionTreeRegressor(max_depth=4, random_state=42),
    name="Decision Tree",
    # Decision Tree is not affected by normalization
    normalize=False,
    # Decision Tree in scikit-learn does not handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

Let’s take a look at the single-step forecast from DecisionTreeRegressor. It’s not doing as well as the linear or regularized linear regression models we have run so far:

Figure 8.12 – Decision tree forecast

Figure 8.12: Decision tree forecast

For the linear models, some coefficients helped us understand how much each feature was important to the prediction function. In decision trees, we don’t have any coefficients, but the feature importance is still estimated using the mean decrease in the loss function, which is attributed to each feature in the tree construction process. This can be accessed in scikit-learn models by using the feature_importance_ attribute of the trained model. Let’s take a look at this feature importance:

Figure 8.13 – Feature importance of a decision tree (top 15)

Figure 8.13: Feature importance of a decision tree (top 15)

Best practice:

Although the default feature importance is a quick and easy way to check how the different features are used, due diligence should be applied before using them for any other purposes, such as feature selection or making business decisions. This way of assessing feature importance gives misleadingly high values for some continuous features and high cardinality categorical features. It is recommended to use permutation importance (sklearn.inspection.permutation_importance) for an easy but better assessment of feature importance. The Further reading section contains some resources regarding the interpretability of models, which can be a good start to understanding what influences the models.

Here, we can see that the important features such as the lag and seasonal rolling features are coming up at the top.

We talked about overfitting and underfitting in Chapter 5, Time Series Forecasting as Regression. These are also referred to as high bias (underfitting) and high variance (overfitting) in machine learning parlance (the Further reading section contains links if you wish to read up more about bias and variance and the trade-off between them).

A decision tree is an algorithm that is highly prone to overfitting or high variance because, unlike the linear function, if given enough expressiveness, it can memorize the training dataset by partitioning the feature space. Another key disadvantage is a decision tree’s inability to extrapolate. Let’s consider a feature, f, that linearly increases our target variable, y. The training data we have has fmax as the maximum value for f and ymax as the maximum value for y. Since the decision tree partitions the feature space and assigns a constant value for that partition, even if we provide f > fmax, we will still only get a prediction of $\hat{y} \leq y_{\max}$.

Now, let’s look at a model that uses decision trees, but in an ensemble, and doesn’t overfit as much.

Random forest

Random Forest is an ensemble learning method that builds multiple decision trees during training and merges their results for improved accuracy and robustness. It excels in both classification and regression tasks by reducing overfitting and enhancing predictive performance through bagging and feature randomness.

Ensemble learning is a process in which we use multiple models, or experts, and combine them in a way to solve the problem at hand. It taps into the wisdom of the crowd approach, which suggests that the decision-making of a group of people is typically better than any individual in that group. In the machine learning context, these individual models are called base learners. A single model may not perform well because it’s overfitting the dataset, but when we combine multiple such models, they can form a strong learner.

Bagging is a form of ensemble learning where we use bootstrap sampling (sampling repeatedly with replacement from a population) to draw different subsets of the dataset, train weak learners on each of these subsets, and combine them by averaging or voting (for regression and classification, respectively). Bagging works best for high-variance, low-bias weak learners and the decision tree is a prime successful candidate with bagging. Theoretically, bagging maintains the same level of bias on the weak learners but reduces the variance, resulting in a better model. But if the weak learners are correlated with each other, the benefits of bagging will be limited.

In 2001, Leo Brieman proposed Random Forest, which substantially modifies standard bagging by building a large collection of decorrelated trees. He proposed to alter the tree-building procedure slightly to make sure all the trees that are grown on bootstrapped datasets are not correlated with each other.

Reference check:

The original research paper for Random Forest is cited in the References section as reference 1.

In the Random Forest algorithm, we decide how many trees to build. Let’s call that M trees. Now, for each tree, the following steps are repeated:

  1. Draw a bootstrap sample from the training dataset.
  2. Select f features at random from all the features.
  3. Pick the best split just using f features and split the node into two child nodes.
  4. Repeat steps 2 and 3 until we hit any of the defined stopping criteria.

This set of M trees is the Random Forest. The key difference here from regular trees is the random sampling of features at each split, which increases randomness and reduces the correlation in the outputs of different trees. While predicting, we use each of these M trees to get a prediction. For regression problems, we average them, while for classification problems, we take the majority vote. The final prediction function that we learn from the Random Forest for regression is as follows:

$$\hat{y} = \frac{1}{M} \sum_{t=1}^{M} T_t(x)$$

Here, Tt(x) is the output of the tth tree in the Random Forest.

All the hyperparameters that we have to control the complexity of the decision tree are applicable here as well (RandomForestRegressor from scikit-learn). In addition to those, we have two other important parameters – the number of trees to build in the ensemble (n_estimators) and the number of features randomly chosen for each split (max_features).

Now, let’s see how we can use Random Forest and evaluate the fit on a sample household from our validation dataset:

from sklearn.ensemble import RandomForestRegressor
model_config = ModelConfig(
    model=RandomForestRegressor(random_state=42, max_depth=4),
    name="Random Forest",
    # RandomForest is not affected by normalization
    normalize=False,
    # RandomForest in scikit-learn does not handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

Let’s take a look at this single-step forecast from RandomForestRegressor. It’s better than the decision tree, but it’s not as good as the linear models. However, we should keep in mind that we have not tuned the model and may be able to get better results by setting the right hyperparameters.

Now, let’s take a look at the forecast that was generated using Random Forest:

Figure 8.14 – Random Forest forecast

Figure 8.14: Random Forest forecast

Just like the feature importance in decision trees, Random Forests also have a very similar mechanism for estimating the feature importance. Since we have a lot of trees in the Random Forest, we accumulate the decrease in split criterion across all the trees in the forest and arrive at a single feature of importance for the Random Forest. This can be accessed in scikit-learn models by using the feature_importance_ attribute of the trained model. Let’s take a look at the feature importance:

Figure 8.15 – Feature importance of a decision tree (top 15)

Figure 8.15: Feature importance of a decision tree (top 15)

Here, we can see that the feature importance is very similar to decision trees. The same caveat about this kind of feature importance applies here as well. This is just a quick and dirty way of looking at what the model is using internally.

Typically, Random Forest achieves good performance on many datasets with very little tuning, so Random Forests are a very popular option in machine learning. The fact that it is difficult to overfit with a Random Forest also increases their appeal. But since Random Forest uses decision trees as the weak learners, the inability of decision trees to extrapolate is passed down to Random Forest as well.

The scikit-learn implementation of Random Forest can get a bit slow for a large number of trees and data sizes. The XGBRFRegressor from the XGBoost library offers an alternative implementation of a Random Forest that can be faster, especially on larger datasets, due to XGBoost's optimized algorithms and parallelization capabilities. Moreover, XGBRFRegressor uses similar hyperparameters to those in the scikit-learn Random Forest, making it relatively straightforward to switch between implementations while tuning the model. In most cases, this is a drop-in replacement and gives almost the same results. The minor difference is due to small implementation details. We have used this variant in the notebooks as well. This variant is preferred going forward because of obvious runtime considerations. It also handles missing values natively and saves us from an additional preprocessing step. More details about the implementation and how to use it can be found at https://xgboost.readthedocs.io/en/latest/tutorials/rf.html.

Now, let’s look at one last family of functions that is one of the most powerful learning methods and has been proven exceedingly well in a wide variety of datasets–gradient boosting.

Gradient boosting decision trees

Boosting, like bagging, is another ensemble method that uses a few weak learners to produce a powerful committee of models. The key difference between bagging and boosting is in the way the weak learners are combined. Instead of building different models in parallel on bootstrapped datasets, as bagging does, boosting uses the weak learners in a sequential manner, with each weak learner applied to repeatedly modified versions of the data.

To understand the additive function formulation, let’s consider this function:

$$F(x) = 25 + x^2 + \cos(x)$$

We can break this function into f1(x)= 25, f2(x)= x2, f3(x)= cos(x) and rewrite F(x) as follows:

F(x) = f1(x) + f2(x) +f3(x)

This is the kind of additive ensemble function we are learning in boosting. Although, in theory, we can use any weak learner, decision trees are the most popular choice. So, let’s use decision trees to explore how gradient boosting works.

Earlier, when we were discussing decision trees, we saw that a decision tree that has M partitions, P1, P2, …, PM, is as follows:

$$T(x) = \sum_{m=1}^{M} c_m I(x \epsilon P_m)$$

Here, x is the input, Cm is the constant response for the region, Pm, and I is a function that is 1 if $x \epsilon P_m$; otherwise, it is 0. A boosted decision tree model is a sum of such trees:

$$\hat{y} = \sum_{k=1}^{M} T_k(x)$$

Since finding the optimal partitions, P, and the constant value, c, for all the trees in the ensemble is a very difficult optimization problem, we usually adopt a suboptimal, stagewise solution where we optimize each step as we build the ensemble. In gradient boosting, we use the gradient of the loss to direct our optimization, hence the name.

Let the loss function we are using in the training be $L(\hat{y}, y)$. Since we are looking at a stagewise additive functional form, we can replace $\hat{y}_k$ with $\hat{y}_{k-1} + T_k(x)$, where $\hat{y}_{k-1}$ is the prediction of the sum of all trees until k-1 and Tk(x) is the prediction of the tree at stage k. Let’s look at what the gradient boosting learning procedure for training data D with N samples is:

  1. Initialize the model with a constant value by minimizing the loss function:

$$F_{0(x)} = \underset{b_0}{\operatorname{argmin}} \sum_{i=1}^{N} L(y_i, b_0)$$

  1. Now that we have the initial solution, we can start the tree-building process. For k=1 to M, we must do the following:
    1. Compute $r_k = -\left[ \frac{\delta L(y, F_{k-1}(x))}{\delta F_{k-1}(x)} \right]$ for all the training samples:
      • rk is the derivative of the loss function with respect to F(x) from the last iteration. It’s also called pseudo-residuals.
      • For squared error loss, this is just the residual, ($\hat{y} - y$).
    2. Build a regular regression tree to the rk values with Mk partitions or leaf nodes, Pmk.
    3. Compute $\rho_t = \underset{\rho}{\operatorname{arg\,min}} \sum_{i=1}^{N} L(y_i, F_{k-1}(x_i) + \rho T_k(x_i))$
      • $\rho_k$ is the scaling factor of the leaf or partition values for the current stage.
      • $T_k(x_i)$ is the function that was learned by the decision tree from the current stage.
    4. Update $F_k(x) = F_{k-1}(x) + \eta + \rho_k \times T_k(x_i)$
      • $\eta$ is the shrinkage parameter or learning rate.

This process of “boosting” the errors of the previous weak model gives the algorithm its name—gradient boosting, where the gradient here means the residual on the previous weak model.

Boosting, typically, is a high-variance algorithm. This means that the chance of overfitting the training dataset is quite high and that enough measures need to be taken to make sure it doesn’t happen. There are many ways regularization and capacity constraining have been implemented in gradient-boosted trees. As always, all the key parameters that decision trees have to reduce capacity to fit the data are valid here because the weak learner is a decision tree. In addition to that, there are two other key parameters – the number of trees, M (n_estimators in scikit-learn), and the learning rate, $\eta$ (learning_rate in scikit-learn).

When we apply a learning rate in the additive formulation, we are essentially shrinking each weak learner, thus reducing the effect of any one weak learner on the overall function. This was originally referred to as shrinkage, but now, in all the popular implementations of gradient-boosted trees, it is referred to as the learning rate. The number of trees and the learning rate are highly interdependent. For the same problem, we will need a greater number of trees if we reduce the learning rate. It has been empirically shown that a lower learning rate improves the generalization error. Therefore, a very effective and convenient way is to set the learning rate to a very low value (<0.1), set a very high value for the number of trees (>5,000), and train the gradient-boosted tree with early stopping. Early stopping is when we use a validation dataset to monitor the out-of-sample performance while training the model. We stop adding more trees to the ensemble when the out-of-sample error stops reducing.

Another key technique a lot of the implementations adopt is subsampling. Subsampling can be done on rows and columns. Row subsampling is similar to bootstrapping, where each candidate in the ensemble is trained on a subsample of the dataset. Column subsampling is similar to random feature selection in Random Forest. Both of these techniques introduce a regularization effect to the ensemble and help reduce generalization errors. Some implementations of gradient-boosted trees, such as XGBoost and LightGBM, implement L1 and L2 regularization directly in the objective function as well.

There are many implementations of regression gradient-boosted trees. A few popular implementations are as follows:

Each of these implementations offers changes that range from subtle to very fundamental regarding the standard gradient boosting algorithm. We have included a few resources in the Further reading section so that you can read up on these differences and get acquainted with the different parameters they support.

For our exercise, we are going to use LightGBM from Microsoft Research because it is one of the fastest and best-performing implementations. LightGBM and CatBoost also support categorical features out of the box and handle missing values natively.

Reference check:

The original research papers for XGBoost, LightGBM, and CatBoost are cited in the References section as 2, 3, and 4, respectively.

Now, let’s see how we can use LightGBM and evaluate the fit on a sample household from our validation dataset:

from lightgbm import LGBMRegressor
model_config = ModelConfig(
    model=LGBMRegressor(random_state=42),
    name="LightGBM",
    # LightGBM is not affected by normalization
    normalize=False,
    # LightGBM handles missing values
    fill_missing=False,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
)

Let’s take a look at the single-step forecast from LGBMRegressor. It’s already significantly better than all the other models we have tried so far:

Figure 8.16 – LightGBM forecast

Figure 8.16: LightGBM forecast

Just like the feature importance in decision trees, gradient-boosting implementations also have a very similar mechanism for estimating the feature importance. The feature importance for the ensemble is given by the average of split criteria reduction attributed to each feature in all the trees. This can be accessed in the scikit-learn API as the feature_importance_ attribute of the trained model. Let’s take a look at the feature importance:

Figure 8.17 – Feature importance of LightGBM (top 15)

Figure 8.17: Feature importance of LightGBM (top 15)

There are multiple ways of getting feature importance from the model, and each implementation has slightly different ways of calculating it. This is controlled by parameters. The most common ways of extracting it (sticking to LightGBM terminology) are split and gain. If we choose split, the feature importance is the number of times a feature is used to split nodes in the trees. On the other hand, gain is the total reduction in the split criterion. This can be attributed to any feature. Figure 8.17 shows split, which is the default value in LightGBM. We can see that the order of the feature importance is very much similar to decision trees, or Random Forests, with almost the same features taking the top three spots.

Gradient boosted decision trees (GBDTs) typically give us very good performance on tabular data and time series as regression is no exception. This very strong model has usually been part of almost all winning entries in Kaggle competitions on time series forecasting in the recent past. While it is one of the best machine learning model families, it still has a few disadvantages:

Table 8.1 – Parallelization strategies in LightGBM

Figure 8.18: Parallelization strategies in LightGBM

To summarize, let’s look at the metrics and runtime that were taken by these machine learning models. If you have run the notebook along with this chapter, then you will find the following summary table there as well:

Figure 8.18 – Summary of the metrics and runtimes for a sample household

Figure 8.18: Summary of the metrics and runtimes for a sample household

Right off the bat, we can see that all of the machine learning models we tried have performed better than the baselines in all metrics except the forecast bias. The three linear regression models perform well with almost equal performance on MAE, MASE, and MSE, with a slight increase in runtimes for regularized models. The decision tree has underperformed, but this is usually expected. Decision trees need to be tuned a little better to reduce overfitting. Random Forest (both the scikit-learn and XGBoost implementations) has improved the decision tree’s performance, which is what we would expect. One key thing to note here is that the XGBoost implementation of Random Forest is almost six times faster than the scikit-learn one. Finally, LightGWM has the best performance across all metrics and a faster runtime.

Now, this was just one household out of all the selected ones. To see how well these models are doing, we need to evaluate them on all selected households.

Training and predicting for multiple households

We have picked a few models (LassoCV, XGBRFRegressor, and LGBMRegressor) that are doing better in terms of metrics, as well as runtime, to run on all the selected households in our validation dataset. The process is straightforward: loop over all the unique combinations, inner loop over the different models to run, and then train, predict, and evaluate. The code is available in the 01-Forecasting_with_ML.ipynb notebook in Chapter08, under the Running an ML Forecast For All Consumers heading. You can run the code and take a break because this is going to take a little less than an hour. The notebook also calculates the metrics and contains a summary table that will be ready for you when you’re back.

Let’s look at the summary now:

Figure 8.19 – Aggregate metrics on all the households in the validation dataset

Figure 8.19: Aggregate metrics on all the households in the validation dataset

Here, we can see that even at the aggregated level, the different models we used perform as expected. The notebook also saves the predictions for the validation set on disk.

Notebook alert:

We also need to run another notebook, called 01a-Forecasting_with_ML_for_Test_Dataset.ipynb, in Chapter08. This notebook follows the same process, generates the forecast, and calculates the metrics on the test dataset.

The aggregate metrics for the test dataset are as follows (from the notebook):

Figure 8.20 – Aggregate metrics on all the households in the test dataset

Figure 8.20: Aggregate metrics on all the households in the test dataset

In Chapter 6, Feature Engineering for Time Series Forecasting, we used AutoStationaryTransformer (not the Transformer model, which we will learn about in Chapter 14) on all the households and saved the transformed dataset.

Using AutoStationaryTransformer

The process is really similar to what we did earlier in this chapter, but with small changes. We read in the transformed targets and joined them to our regular dataset in such a way that the original target is named energy_consumption and the transformed target is named energy_consumption_auto_stat:

#Reading the missing value imputed and train test split data
train_df = pd.read_parquet(preprocessed/"block_0-7_train_missing_imputed_feature_engg.parquet")
auto_stat_target = pd.read_parquet(preprocessed/"block_0-7_train_auto_stat_target.parquet")
transformer_pipelines = joblib.load(preprocessed/"auto_transformer_pipelines_train.pkl")
#Reading in validation as test
test_df = pd.read_parquet(preprocessed/"block_0-7_val_missing_imputed_feature_engg.parquet")
# Joining the transformed target
train_df = train_df.set_index(['LCLid','timestamp']).join(auto_stat_target).reset_index()

And while defining FeatureConfig, we used energy_consumption_auto_stat as target and energy_consumption as original_target.

Notebook alert:

The 02-Forecasting_with_ML_and_Target_Transformation.ipynb and 02a-Forecasting_with_ML_and_Target_Transformation_for_Test_Dataset.ipynb notebooks use these transformed targets to generate the forecasts for the validation and test datasets, respectively.

Let’s look at the summary metrics that were generated by these notebooks on the transformed data:

Figure 8.21 – Aggregate metrics on all the households with transformed targets in the validation dataset

Figure 8.21: Aggregate metrics on all the households with transformed targets in the validation dataset

The target transformed models are not performing as well as the original ones. This might be because the dataset doesn’t have any strong trends.

Congratulations on making it through a very heavy and packed chapter full of theory as well as practice. We hope this has enhanced your understanding of machine learning and ability to apply these modern techniques to time series data.

Summary

This was a very practical and hands-on chapter in which we developed some standard code to train and evaluate multiple machine learning models. Then, we reviewed a few key machine learning models like ridge regression, lasso regression, decision trees, Random Forest, and gradient-boosted trees and how they work behind the hood. To complete and reinforce what we learned, we applied the machine learning models we learned about to the London Smart Meters dataset and saw how well they did. This chapter sets you up to tackle the coming chapters, where we will use the standardized code and these models to go deeper into forecasting with machine learning.

In the next chapter, we will start combining different forecasts into a single forecast and explore concepts such as combinatorial optimization and stacking to achieve state-of-the-art results.

References

The following references were provided in this chapter:

  1. Breiman, L. Random Forests, Machine Learning 45, 5–32 (2001): https://doi.org/10.1023/A:1010933404324.
  2. Chen, Tianqi and Guestrin, Carlos. (2016). XGBoost: A Scalable Tree Boosting System. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD ‘16). Association for Computing Machinery, New York, NY, USA, 785–794: https://doi.org/10.1145/2939672.2939785.
  3. Ke, Guolin et.al. (2017), LightGBM: A Highly Efficient Gradient Boosting Decision Tree. Advances in Neural Information Processing Systems, pages 3149-3157: https://dl.acm.org/doi/pdf/10.5555/3294996.3295074.
  4. Prokhorenkova, Liudmila et al. (2018), CatBoost: unbiased boosting with categorical features. Proceedings of the 32nd International Conference on Neural Information Processing Systems (NIPS’18): https://dl.acm.org/doi/abs/10.5555/3327757.3327770.

Further reading

To learn more about the topics that were covered in this chapter, take a look at the following resources: