10

Global Forecasting Models

In previous chapters, we saw how we can use modern machine learning models on time series forecasting problems, essentially replacing traditional models such as ARIMA or exponential smoothing. However, before now, we were looking at the different time series in any dataset (such as households in the London Smart Meters dataset) in isolation, just as the traditional models did.

However, we will now explore a different paradigm of modeling where we use a single machine learning model to forecast a bunch of time series together. As we will learn in the chapter, this paradigm brings many benefits with it, from the perspective of both computation and accuracy.

In this chapter, we will be covering these main 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 associated code for the chapter can be found at https://github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-/tree/main/notebooks/Chapter10.

Why Global Forecasting Models?

We talked about global models briefly in Chapter 5, Time Series Forecasting as Regression, where we mentioned related datasets. We can think of many scenarios where we would encounter related time series. We may need to forecast the sales for all the products of a retailer, the number of rides requested for a cab service across different areas of a city, or the energy consumption of all the households in a particular area (which is what the London Smart Meters dataset does). We call these related time series because all the different time series in the dataset can have a lot of factors in common with each other. For instance, the yearly seasonality that might occur in retail products might be present for a large section of products, or the way an external factor such as temperature affects energy consumption may be similar for a large number of households. Therefore, one way or the other, the different time series in a related time series dataset share attributes between them.

Traditionally, we used to consider each time series an independent time series; in other words, each time series was assumed to be generated using a different data-generating process. Classical models such as ARIMA and exponential smoothing are trained for each time series. However, we can also consider all the time series in the dataset as being generated from a single data-generating process, and the subsequent modeling approach would be to train a single model to forecast all the time series in the dataset. The latter is what we refer to as Global Forecasting Models (GFMs). GFMs are models that are designed to handle multiple related time series, allowing for shared learning across those time series. In contrast, the traditional approach is referred to as Local Forecasting Models (LFMs).

Although we briefly talked about the drawbacks of LFMs in Chapter 5, Time Series Forecasting as Regression, let’s summarize them in a more concrete fashion and see why GFMs help us smooth over a lot of those drawbacks.

Sample size

In most real-world applications (especially in business forecasting), the time series we have to forecast is not very long. Adopting a completely data-driven approach to modeling such a small time series is problematic. Training a highly flexible model with a handful of data points will lead to the model memorizing the training data, resulting in an overfit.

Traditionally, this has been overcome by placing strong priors or inductive bias into the models we use for forecasting. Inductive bias loosely refers to a set of assumptions or restrictions that are built into a model that should help the model predict feature combinations it has not encountered while training. For instance, double exponential smoothing has strong assumptions about seasonality and trend. The model does not allow any other more complicated patterns to be learned from the data. Therefore, using these strong assumptions, we are restricting the model search to a small section of the hypothesis space. While this helps in low-data regimes, the flip side is that these assumptions may limit accuracy.

Recent developments in the field of machine learning have shown us without a doubt that using a data-driven approach (with much fewer assumptions or priors) on large training sets will lead to us training better models. However, conventional statistical wisdom tells us that the number of data points needs to be at least 10 to 100 times the number of parameters that we are trying to learn from those data points.

So, if we stick to LFMs, scenarios in which we can adopt a completely data-driven approach will be very rare. This is where GFMs shine. A GFM is able to use the history of all the time series in a dataset to train the model and learn a single set of parameters that work for all the time series in the dataset. Borrowing the terminology introduced in Chapter 5, Time Series Forecasting as Regression, we increase the width of the dataset, keeping the length the same (refer back to Figure 5.2). This explosion of historical information available to a single model lets us use completely data-driven techniques on time series datasets.

Cross-learning

GFMs, by design, promote cross-learning across different time series in a dataset. Imagine we have a time series that is quite new and does not have a history rich enough for teaching the model—for instance, the sales of a newly introduced retail product or the electricity consumption of a new household in a region. If we consider these time series in isolation, it will be a while before we start to get reasonable forecasts from the models we train on them, but GFMs make that process easier by enabling cross-learning. GFMs have an implicit sense of similarity between different time series and they will be able to use patterns they have seen in similar time series with a rich history to come up with a forecast on the new time series.

Another way cross-learning helps is by acting like a regularizer while estimating common parameters such as seasonality. For instance, the seasonality exhibited by similar products in a retail scenario is best estimated at an aggregate level, because each individual time series will have some sort of noise that can creep into the seasonality extraction. By enforcing common seasonality across multiple products, we are essentially regularizing the seasonality estimation and, in the process, making the seasonality estimate more robust. The good thing about GFMs is that they take a data-driven approach to define the seasonality of which products should be estimated together and which ones have different patterns. If you have different seasonality patterns in different products, a GFM may struggle to model them together. However, when provided with enough information on how to distinguish between different products, the GFM will be able to learn that difference too.

Multi-task learning

GFMs can be considered multi-task learning paradigms where a single model is trained to learn multiple tasks (as forecasting each time series is a separate task). Multi-task learning is an active area of research, and there are many benefits to using multi-task models:

Engineering complexity

LFMs pose a challenge from the engineering side as well for large datasets. If we have thousands or millions of time series to forecast, it becomes increasingly difficult to both train and manage the life cycle of these LFMs. In Chapter 8, Forecasting Time Series with Machine Learning Models, we trained LFMs for just a subset of households in the dataset. It took almost 20 to 30 minutes to train a machine learning model for all 150 households and we ran them with the default hyperparameters. In a normal machine learning workflow, we train multiple machine learning models and do hyperparameter tuning to find the best configuration of the model. However, carrying out all these steps for thousands of time series in a dataset becomes increasingly complex and time-consuming.

Equally, then there is the issue of managing the life cycle of these models. All these individual models need to be deployed to production, their performance needs to be monitored to check for model and data drift, and they need to be retrained at a set frequency. This becomes increasingly complex as we have more and more time series to forecast.

However, by shifting to a GFM paradigm, we drastically reduce the time and effort required to train and manage a machine learning model throughout its life cycle. As we will see in this chapter, training a GFM on these 150 households takes only a fraction of the time it takes to train LFMs.

Despite all the advantages of GFMs, they are not without some drawbacks. The main drawback is that we are assuming that all the time series in a dataset are generated by a single Data Generating Process (DGP). This might not be a valid assumption and this can lead to the GFM underfitting some specific types of time series patterns that are underrepresented in the dataset.

Another open issue is whether a GFM is good for use with unrelated tasks or time series. The jury is out on this one, but Montero-Manso et al. proved that there are also gains in modeling unrelated time series with a GFM. The same finding has been put forward, although from another perspective, by Oreshkin et al., who trained a global model on the M4 dataset (a set of unrelated datasets) and obtained state-of-the-art performance. They attributed it to the meta-learning capabilities of the model.

That being said, relatedness does help the GFM, as the learning task becomes easier this way. We will see a practical application of this in upcoming sections of this chapter as well.

In the larger scheme of things, the benefits we derive from a GFM paradigm far outweigh the drawbacks. On most tasks, the GFMs either perform on par with or better than local models. It has been proven theoretically as well, by Montero-Manso et al., that a GFM, in a worst-case scenario, learns the same function as a local model. We will see this clearly in the models we are going to train in the upcoming sections. Finally, the training time and engineering complexity drop drastically as you move to a GFM paradigm.

Now that we have explained why a GFM is a worthwhile paradigm to adopt, let’s see how we can train one.

Creating GFMs

Training a GFM is very straightforward. While we were training LFMs in Chapter 8, Forecasting Time Series with Machine Learning Models, we were looping over different households in the London Smart Meters dataset and training a model for each household. However, if we just take all the households into a single dataframe (our dataset is already that way) and train a single model on it, we get a GFM. One thing we want to keep in mind is to make sure that all the time series in the dataset have the same frequency. In other words, if we mix daily time series with weekly ones while training these models, the performance drop will be noticeable—especially if we are using time-varying features and other time-based information. For a purely autoregressive model, mixing time series in this way is much less of a problem.

Notebook alert:

To follow along with the complete code, use the notebook named 01-Global_Forecasting_Models-ML.ipynb in the Chapter10 folder.

The standard framework we developed in Chapter 8, Forecasting Time Series with Machine Learning Models, is general enough to work for GFMs as well. So, as we did in that chapter, we define FeatureConfig and MissingValueConfig in the 01-Global_Forecasting_Models-ML.ipynb notebook. We also slightly tweaked the Python function to train and evaluate the machine learning to make it work for all households. The details and exact functions can be found in the notebook.

Now, instead of looping over different households, we input the entire training dataset into the get_X_y function:

# Define the ModelConfig
from lightgbm import LGBMRegressor
model_config = ModelConfig(
    model=LGBMRegressor(random_state=42),
    name="Global LightGBM Baseline",
    # LGBM is not sensitive to normalized data
    normalize=False,
    # LGBM can handle missing values
    fill_missing=False,
)
# Get train and test data
train_features, train_target, train_original_target = feat_config.get_X_y(
    train_df, categorical=True, exogenous=False
)
test_features, test_target, test_original_target = feat_config.get_X_y(
    test_df, categorical=True, exogenous=False
)

Now that we have the data, we need to train the model. Training the model is also exactly the same as we saw in Chapter 8, Forecasting Time Series with Machine Learning Models. We will just choose LightGBM, which was the best-performing LFM model, and use functions we defined earlier to train the model and evaluate the results:

y_pred, feat_df = train_model(
        model_config,
        _feat_config,
        missing_value_config,
        train_features,
        train_target,
        test_features,
    )
agg_metrics, eval_metrics_df = evaluate_forecast(
    y_pred, test_target, train_target, model_config
)

Now, in y_pred, we will have the forecast for all the households and feat_df will have the feature importance. agg_metrics will have the aggregated metric for all the selected households.

Let’s look at how well our GFM model did:

Figure 10.1 – Aggregate metrics with the baseline GFM

Figure 10.1: Aggregate metrics with the baseline GFM

We are not doing better than the best LFM (in the first row) in terms of the metrics. However, one thing we should note is the time taken to train the model—~30 seconds. The LFM for all the selected households was taking ~30 minutes. This huge reduction in time taken gives us a lot of flexibility to iterate faster with different features and techniques.

With that said, let’s now look at a few techniques with which we can improve the accuracy of the GFMs.

Strategies to improve GFMs

GFMs have been in use in many forecasting competitions in Kaggle and outside of it. They have been battle-tested empirically, although very little work has gone into examining why they work so well from a theoretical point of view. Montero-Manso and Hyndman (2020) have a working paper titled Principles and Algorithms for Forecasting Groups of Time Series: Locality and Globality, which is an in-depth investigation, both theoretical and empirical, of GFMs and the many techniques that have been developed by the data science community collectively. In this section, we will try to include strategies to improve GFMs and, wherever possible, try to give theoretical justifications for why they would work.

Reference check:

The Montero-Manso and Hyndman (2020) research paper is cited in References as reference 1.

In the paper, Montero-Manso and Hyndman use a basic result in machine learning about generalization error to carry out the theoretical analysis, and it is worth spending a bit of time understanding that, at least at a high level. Generalization error, as we know, is the difference between out-of-sample error and in-sample error. Yaser S Abu-Mostafa has a free, online Massive Open Online Course (MOOC) and an associated book (both of which are linked in the Further reading section). It is a short course on machine learning and is a course that I would recommend to anyone in the machine learning field for developing a stronger theoretical and conceptual basis for what we do. One of the important concepts the course and book put forward is the use of Hoeffding’s inequality from probability theory to derive bounds on a learning problem. Let’s quickly look at the result to develop our understanding:

$$E_{out} \lt E_{in} + \sqrt{\frac{\log(|H|) + \log\left(\frac{2}{\delta}\right)}{2N}}$$

It has a probability of at least 1-$\delta$.

Ein is the in-sample average error and Eout is the expected out-of-sample error. N is the total number of samples in the dataset from which we are learning and H is the hypothesis class of models. It is a finite set of functions that can potentially fit the data. The size of H, denoted by |H|, is the complexity of H. Although the formula of the bound looks intimidating, let’s simplify the way we look at it to develop the necessary understanding.

We want Eout to be as close to Ein as possible, and for that, we need the terms in the square root to be as small as possible. There are two terms under the square root that are in our control, so to speak—N and |H|. Therefore, to make the generalization error (Ein- Eout) as small as possible, we either need to increase N (have more data) or decrease |H| (have a less complex model). This is a result that is applicable to all machine learning but Montero-Manso and Hyndman, with a few assumptions, made this applicable to time series models as well. It is this result that they used to give theoretical backing to the arguments put forward in their working paper.

Montero-Manso and Hyndman have taken Hoeffding’s inequality and applied it to LFMs and GFMs to compare them. We can see the result here (for a full mathematical and statistical understanding, refer to the original paper under References):

$$E_{out}^{Local} \lt E_{in}^{Local} + \sqrt{\frac{log(\prod_{i=1}^{K} |H_i|) + log\left(\frac{2}{\delta}\right)}{2NK}}$$

$$E_{out}^{Global} \lt E_{in}^{Global} + \sqrt{\frac{log(|J|) + log\left(\frac{2}{\delta}\right)}{2NK}}$$

$E_{out}^{Local}$ and $E_{out}^{Global}$ are the average in-sample errors across all the time series using the local and global approaches, respectively. $E_{out}^{Global}$ and $E_{in}^{Global}$ are the out-of-sample expectations under the local and global approaches, respectively. Hi is the hypothesis class for the i-th time series and J is the hypothesis class for the global approach (the global approach only fits a single function and hence, has just a single hypothesis class).

One of the most interesting results that comes out of this is that the complexity term for LFMs ($\log(\prod_{i=1}^{K} |H_i|)$) grows the size of the dataset. The greater the number of time series we have in the dataset, the more complexity and the worse the generalization error, whereas with GFMs, the complexity term (log(|J|)) stays constant. Therefore, for a dataset of moderate size, the overall complexity of LFMs (such as exponential smoothing) can be much higher than a single GFM, no matter how complex the GFM is. As a corollary, we can also think that with the available dataset (NK), we can afford to train a model with much higher complexity than a model for LFMs. There are many ways to increase the complexity of the model, which we will see in the following section.

Now, let’s return to the GFMs we were training. We saw that the performance of the GFM we trained was not up to the mark when we compared it with the best LFM (LightGBM), but it is better than the baseline and other models we tried, so right off the bat, we know the GFM we trained is not terrible. Now, let’s look at a few ways to improve the performance of the model.

Increasing memory

As we discussed in Chapter 5, Time Series Forecasting as Regression, the machine learning models that we discuss in this book are finite memory models or Markov models. A model such as exponential smoothing takes into account the entire history of a time series while forecasting, but models such as any of the machine learning models we discussed only take in finite memory to make their predictions. In a finite memory model, the amount of memory we allow the model to access is called the size of the memory (M) or order of autoregression (from econometrics).

Providing a greater amount of memory to the model increases the complexity of the model. Therefore, one of the ways to increase the performance of the GFM is to increase the amount of memory the model has access to. There are many ways to increase the amount of memory.

Adding more lag features

If you have prior exposure to ARIMA models, you will know that the number of Autoregressive (AR) terms are sparingly used. We usually see AR models with single-digit lags. There is nothing stopping us from running ARIMA models with larger lags, but since we do run ARIMA in the LFM paradigm, the model has to learn the parameters of all the lags using limited data and therefore, in practice, practitioners commonly choose smaller lags. However, when we are moving to GFMs, we can afford to have much larger lags. Montero-Manso and Hyndman empirically showed the benefits of adding more lags to GFMs. For highly seasonal time series, a peculiar phenomenon was observed. The accuracy improves with an increase in lags, but it then saturates and suddenly worsens when the lag becomes equal to the seasonal cycle. On further increasing the lags beyond the seasonal cycle, the accuracy shows huge gains. This may be because of the overfitting that happens because of seasonality. It becomes very easy for the model to favor the seasonal lag because it works very well in a sample, so it’s better to add a few more lags on the plus side of the seasonal cycle.

Adding rolling features

Another way to increase the memory of the model is to include rolling averages as features. Rolling averages take information from extended windows on memory and encode that information by way of descriptive statistics (such as the mean or max). This is an efficient way of including the memory because we can take very large windows for memory and include the information as a single feature in the model.

Adding EWMA features

An Exponentially Weighted Moving Average (EWMA) is a way to include infinite memory in a finite memory model. The EWMA essentially takes the average of the entire history but is weighted according to the $\alpha$ that we set. Therefore, with different values of $\alpha$, we get different kinds of memory, again encoded as a single feature. Including different EWMA features has also empirically proved beneficial.

We have already included these kinds of features in our feature engineering (Chapter 6, Feature Engineering for Time Series Forecasting), and they are part of the baseline GFM we trained, so let’s move on to the next strategy for improving the accuracy of GFMs.

Using time series meta-features

The baseline GFM we trained earlier in the Creating Global Forecasting Models (GFMs) section had lag features, rolling features, and EWMA features, but we have given no feature that helps the model distinguish between different time series in the dataset. The baseline GFM model learned a generalized function that generates a forecast provided the features. This might work well enough for homogenous datasets where all the time series are very similar in nature, but for heterogenous datasets, the information with which the model can distinguish each time series comes in handy.

So, information about the time series itself is what we call meta-features. In a retail context, it can be the product ID, the category of products, the store number, and so on. In our dataset, we have features such as stdorToU, Acorn, Acorn_grouped, and LCLid, which give some information about the time series itself. Including these meta-features in the GFM will improve the performance of the model.

However, there is just one problem—more often than not, these meta-features are categorical in nature. A feature is categorical when the values in the feature can only take discrete values. For instance, Acorn_grouped can only have one of three values—Affluent, Comfortable, or Adversity. Most machine learning models do not work well with categorical features. All the models in scikit-learn, the most popular machine learning library in the Python ecosystem, do not allow categorical features at all. To include categorical features in machine learning models, we need to encode them into numerical form, and there are many ways to encode categorical columns. Let’s review a few popular options.

Ordinal encoding and one-hot encoding

The most popular ways of encoding categorical features are ordinal encoding and one-hot encoding, but they are not always the best choices. Let’s quickly review what these techniques are and when they are suitable.

Ordinal encoding is the simplest of them all. We simply assign a numerical code to the unique values of a category and then replace the categorical value with the numerical code. To encode the Acorn_grouped feature from our dataset, all we need to do is assign codes, say 1 for Affluent, 2 for Comfortable, and 3 for Adversity, and replace all instances of the categorical values with the code we assigned. While this is really easy, this kind of encoding introduces meanings to the categorical values that we may or may not intend. When we assign numerical codes, we are implicitly saying that the categorical value that gets assigned 2 as a code is better than the categorical value with 1 as the code. This kind of encoding only works for ordinal features (features whose categorical values have an intrinsic sense of rank in their meaning) and should be sparingly used. Another way we can think about the problem is in terms of distance. When we do ordinal encoding, the distance between Comfortable and Affluent can be higher than the distance between Comfortable and Adversity, depending on the way we encode.

One-hot encoding is a better way of representing categorical features with no ordinal meaning. It essentially encodes the categorical features in a higher dimension, placing the categorical values equally distant in that space. The size of the dimension it requires to encode the categorical values is equal to the cardinality of the categorical variable. Cardinality is the number of unique values in the categorical feature. Let’s see how sample data would be encoded in a one-hot encoding scheme:

Figure 10.2 – One-hot encoding of categorical features

Figure 10.2: One-hot encoding of categorical features

We can see that the resulting encoding will have a column for each unique value in the categorical feature and the value is indicated by 1 in the column. For instance, the first row is Comfortable, and therefore, every other column except the Comfortable column will have 0 and the Comfortable column will have 1. If we calculate the Euclidean distance between any two categorical values, we can see that they are the same.

However, there are three main issues with this encoding, all of which become a problem with high cardinality categorical variables:

There is a slightly different way of one-hot encoding where we drop one of these dimensions, called dummy variable encoding. This has the added benefit of making the encoding linearly independent, which, in turn, has some advantages, especially for vanilla linear regression. The Further reading section has a link if you want to know more.

Since the categorical columns that we must encode have high cardinality (at least a few of them), we will not be doing this encoding. Instead, let’s look at a few encoding techniques that can handle high cardinality categorical variables better.

Frequency encoding

Frequency encoding is an encoding schema that does not increase the dimensions of the problem. It takes a single categorical array and returns a single numeric array. The logic is very simple—it replaces the categorical values with the number of times the value occurs in the training dataset. Although it’s not perfect, this works pretty well, as it lets the model distinguish between different categories based on how frequently they occur.

There is a popular library, category_encoders, that implements a lot of different encoding schemes in a standard scikit-learn style estimator, and we will be using that in our experiments as well. The standard framework we developed in Chapter 8, Forecasting Time Series with Machine Learning Models, also had a couple of functionalities that we didn’t use—encode_categorical and categorical_encoder.

So, let’s use them and train our model now:

from category_encoders import CountEncoder
from lightgbm import LGBMRegressor
#Define which columns names are categorical features
cat_encoder = CountEncoder(cols=cat_features)
model_config = ModelConfig(
    model=LGBMRegressor(random_state=42),
    name="Global LightGBM with Meta Features (CountEncoder)",
    # LGBM is not sensitive to normalized data
    normalize=False,
    # LGBM can handle missing values
    fill_missing=False,
    # Turn on categorical encoding
    encode_categorical=True,
    # Pass the categorical encoder to be used
    categorical_encoder=cat_encoder
)

The rest of the process is the same as what we saw in the Creating Global Forecasting Models (GFMs) section and we get the forecast using the encoded meta-features:

Figure 10.3 – Aggregate metrics with the GFM with meta-features (frequency encoding)

Figure 10.3: Aggregate metrics with the GFM with meta-features (frequency encoding)

Right away, we can see that there is a reduction in error, although it is minimal. We can also see that the training time has almost doubled. This may be because now we have an additional step of encoding the categorical features in addition to training the machine learning model.

The main issue with frequency encoding is that it doesn’t work with features that are uniformly distributed in the dataset. For instance, the LCLid feature, which is just a unique code for each household, is uniformly distributed in the dataset and when we use frequency encoding, all the LCLid features will come to almost the same frequency, and hence the machine learning model considers them almost the same.

Let’s now look at a slightly different approach.

Target mean encoding

Target mean encoding, in its most vanilla form, is a very simple concept. It is a supervised approach that uses the target in the training dataset to encode the categorical columns. Let’s look at an example:

Figure 10.4 – Target mean encoding

Figure 10.4: Target mean encoding

The vanilla target mean encoding has a few limitations. It increases the chance of overfitting the training data because we are using the mean targets directly and thereby leaking the target into the model in a way. Another problem with the approach is that when the categorical values are unevenly distributed, there may be a few categorical values with very small sample sizes, and therefore, the mean estimate becomes noisy. Extending this problem to the extreme, we get another case where an unseen categorical value comes up in test data. This is also not supported in the vanilla version. Therefore, in practice, this simple version is almost never used, but slightly more sophisticated versions of this concept are widely used and are an effective strategy for encoding categorical features.

In category_encoders, there are many variations of this concept, but let’s look at two popular and effective ones here.

In 2001, Daniele Micci-Barreca proposed a variant of mean encoding. If we consider the target as a binary variable, say 1 and 0, the mean (which is the number of 1s or number of samples) is also the probability of having 1. Using this interpretation of the means, Daniele proposed blending two probabilities—prior and posterior probabilities—as the final encoding for the categorical features.

Reference check:

The research paper by Daniele Micci-Barreca is cited in References as reference 2.

The prior probability is defined as follows:

$$P_{prior} = \frac{n_y}{n_{TR}}$$

Here, ny is the number of cases such that target = 1, and nTR is the number of samples in the training data.

The posterior probability is defined for category i as follows:

$$P_{posterior}^{i} = \frac{n_{iY}}{n_i}$$

Here, niY is the number of samples in the dataset where category = i and Y = 1, and ni is the number of samples in the dataset where category = i.

Now, the final encoding for category i is as follows:

$$S_i = \lambda(n_i) \times P_{posterior}^{i} + (1 - \lambda(n_i)) \times P_{prior}$$

Here, $\lambda$ is the weighting factor, which is a monotonically increasing function on ni that is bounded between 0 and 1. So, this function gives a larger weight to the posterior probability as the number of samples increases.

Adapting this to the regression setting, the probabilities change to expected values so that the formula becomes the following:

$$S_i = \lambda(n_i) \times \frac{\sum_{k \in TR_i} Y_k}{n_i} + (1 - \lambda(n_i)) \left\{ \frac{\sum_{k \in TR} Y_k}{n_{TR}} \right\}$$

Here, TRi is all the rows where category = 1 and $\sum_{k \in TR} Y_k$ is the sum of Y for TRi.$\sum_{k \in TR} Y_k$ is the sum of Y for all the rows in the training dataset. As with the binary variable, we are mixing the expected value of Y, given category = i (E[Y|category = i]) and the expected value of Y (E[Y]) for the final categorical encoding.

There are many functions that we can use for $\lambda$. Daniele mentions a very common functional form (sigmoid):

$$\lambda(n_i) = \frac{1}{1 + e^{-\frac{n_i - k}{f}}}$$

Here, ni is the number of samples in the dataset where, category = i and k and f are tunable hyperparameters. k determines half of the minimal sample size for which we completely trust the estimate. If k = 1, what we are saying is that we trust the posterior estimate from a category that has only two samples. f determines how fast the sigmoid transitions between the two extremes. As f tends to infinity, the transition becomes a hard threshold between prior and posterior probabilities. TargetEncoder from category_encoders has implemented this $\lambda$. The k parameter is called min_samples_leaf with a default value of 1, and the f parameter is called smoothing with a default value of 1. Let’s see how this encoding works on our problem. Using a different encoder in the framework we are working on is as simple as passing a different cat_encoder (the initialized categorical encoder) to ModelConfig:

from category_encoders import TargetEncoder
cat_encoder = TargetEncoder(cols=cat_features)

The rest of the code is exactly the same. We can find the full code in the corresponding notebook. Let’s see how well the new encoding has done:

Figure 10.5 – Aggregate metrics with the GFM with meta-features (target encoding)

Figure 10.5: Aggregate metrics with the GFM with meta-features (target encoding)

It’s not doing that well, is it? As with machine learning models, the No Free Lunch Theorem (NFLT) applies to categorical encoding as well. There is no one encoding scheme that works well all the time. Although not directly related to the topic, if you want to know more about the NFLT, head to the Further reading section.

With all these supervised categorical encoding techniques, such as target mean encoding, we have to be really careful not to induce data leakage. The encoder should be fit using training data and not using the validation or test data. Another very popular technique is to generate categorical encoding using cross-validation and use the out-of-sample encodings to absolutely avoid data leakage or overfitting.

There are many more encoding schemes, such as MEstimateEncoder (which uses additive smoothing as the $\lambda$), HashingEncoder, and so on, in category_encoders. Another very effective way of encoding categorical features is using embedding from deep learning. The Further reading section has a link to a tutorial for doing this.

Before now, all this categorical encoding was a separate step before the modeling. Now, let’s look at a technique that considers categorical features natively for model training.

LightGBM’s native handling of categorical features

A few machine learning model implementations handle categorical features natively, especially gradient-boosting models. CatBoost and LightGBM, two of the most popular GBM implementations, handle categorical features out of the box. CatBoost has a unique way of encoding categorical features into numerical ones internally using something similar to additive smoothing. The Further reading section has links to further information on how this encoding is done. category_encoders has implemented this logic as CatBoostEncoder so that we can use this type of encoding for any machine learning model as well.

While CatBoost handles this internal conversion into numerical features, LightGBM takes a more native approach to dealing with categorical features. LightGBM considers the categorical features as is while growing and splitting the trees. For a categorical feature with k unique values (cardinality of k), there are 2k-1-1 possible partitions. This soon becomes intractable, but for regression trees, Walter D. Fisher proposed a technique back in 1958 that makes the complexity of finding an optimal split much less. The essence of the method is to use average target statistics for each categorical value to order them and then find the optimal split in the ordered categorical values.

Reference check:

The research paper by Fisher is cited in References as reference 3.

LightGBM’s scikit-learn API supports this feature by taking an argument, categorical_feature, which has a list of categorical feature names, during fit. We can use the fit_kwargs argument in the fit of our MLModel that we defined in Chapter 8, Forecasting Time Series with Machine Learning Models, to pass in this parameter. Let’s see how we can do this:

from lightgbm import LGBMRegressor
model_config = ModelConfig(
    model=LGBMRegressor(random_state=42),
    name="Global LightGBM with Meta Features (NativeLGBM)",
    # LGBM is not sensitive to normalized data
    normalize=False,
    # LGBM can handle missing values
    fill_missing=False,
    # We are using inbuilt categorical feature handling
    encode_categorical=False,
)
# Training the model and passing in fit_kwargs
y_pred, feat_df = train_model(
    model_config,
    _feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    fit_kwargs=dict(categorical_feature=cat_features),
)

y_pred has the forecasts, which we evaluate as usual. Let’s also see the results:

Figure 10.6 – Aggregate metrics with the GFM with meta-features (native LightGBM)

Figure 10.6: Aggregate metrics with the GFM with meta-features (native LightGBM)

We can observe a good reduction in MAE as well as meanMASE with the native handling of categorical features. We can also see a reduction in the total training time because we don’t have a separate step for encoding the categorical feature. Empirically, the native handling of categorical features works better most of the time.

Now that we have encoded the categorical features, let’s look at another way to improve accuracy.

Tuning hyperparameters

A hyperparameter is a setting that controls how a machine learning model is trained but is not learned from the data. In contrast, model parameters are learned from the data during training. For example, in Gradient Boosting Decision Trees (GBDT), model parameters are the decision thresholds in each tree, learned from the data. Hyperparameters, like the number of trees, learning rate, and tree depth, are set before training and control the model’s structure and how it learns. While parameters adjust based on the data, hyperparameters must be tuned externally.

Although hyperparameter tuning is common practice in machine learning, we haven’t been able to do so because of the sheer number of models we had under the LFM paradigm. Now that we have a GFM that finishes training in 30 seconds, hyperparameter tuning becomes feasible. From a theoretical perspective, we also saw that GFMs can afford a larger complexity and can therefore evaluate a greater number of functions to pick the best without overfitting.

Mathematical optimization is defined as the selection of a best element, with regard to some criterion, from some set of available alternatives. In most cases, this involves finding the maximum or minimum value of some function (an objective function) from a set of alternatives (the search space) subject to some conditions (constraints). The search space can be discrete variables, continuous variables, or a mixture of both, and the objective function can be differentiable or non-differentiable. There is a large body of research that tackles these variations.

You may be wondering why we are talking about mathematical optimization now, right? Hyperparameter tuning is a mathematical optimization problem. The objective function here is non-differentiable and returns the metric for which we are optimizing—for instance, the Mean Absolute Error (MAE).

The search space comprises the different hyperparameters we are tuning—say, the number of trees or depth of the trees. It could be a mixture of continuous and discrete variables and the constraints would be any restriction on the search space we impose—for instance, a particular hyperparameter cannot be negative, or a particular combination of hyperparameters cannot occur. Therefore, being aware of the terms used in mathematical optimization will help us in our discussion.

Even though hyperparameter tuning is a standard machine learning concept, we will quickly review three main techniques (besides manual trial and error) for doing hyperparameter tuning.

Grid search

Grid search can be thought of as a brute-force method where we define a discrete grid over the search space, check the objective function at each point in the grid, and pick the best point in that grid. The grid is defined as a set of discrete points for each of the hyperparameters we choose to tune. Once the grid is defined, all the intersections of the grid are evaluated to search for the best objective value. If we are tuning 5 hyperparameters and the grid has 20 discrete values for each parameter, the total number of trials for a grid search would be 3,200,000 (205). This means training a model 3.2 million times and evaluating it. This becomes quite limiting because most modern machine learning models have many hyperparameters. For instance, LightGBM has more than 100, and out of those, at least 20 are highly impactful parameters when tuned. So, using a brute force approach such as grid search forces us to make the search space quite small so that it becomes feasible to carry out the tuning in a reasonable amount of time.

For our case, we have defined a very small grid of just 27 trials by limiting ourselves to a really small search space. Let’s see how we do that:

from sklearn.model_selection import ParameterGrid
grid_params = {
    "num_leaves": [16, 31, 63],
    "objective": ["regression", "regression_l1", "huber"],
    "random_state": [42],
    "colsample_bytree": [0.5, 0.8, 1.0],
}
parameter_space = list(ParameterGrid(grid_params))

We just tune three hyperparameters (num_leaves, objective, and colsample_bytree), and with just three options for each parameter. Performing the grid search after this is just about looping over the parameter space and evaluating the model at each combination of hyperparameters:

scores = []
for p in tqdm(parameter_space, desc="Performing Grid Search"):
    _model_config = ModelConfig(
        model=LGBMRegressor(**p, verbose=-1),
        name="Global Meta LightGBM Tuning",
        # LGBM is not sensitive to normalized data
        normalize=False,
        # LGBM can handle missing values
        fill_missing=False,
    )
    y_pred, feat_df = train_model(
        _model_config,
        _feat_config,
        missing_value_config,
        train_features,
        train_target,
        test_features,
        fit_kwargs=dict(categorical_feature=cat_features),
    )
    scores.append(ts_utils.mae(
                test_target['energy_consumption'], y_pred
            ))

This takes about 15 minutes to complete and gives us the best MAE of 0.73454, which is already a great improvement from our untuned GFM.

However, this makes us wonder whether there is an even better solution that we haven’t covered in the grid we defined. One option is to expand the grid and run the grid search again. This increases the number of trials exponentially and soon becomes infeasible.

Let’s look at a different method where we can explore a larger search space with the same number of trials.

Random search

Random search takes a slightly different route. In random search, we also define the search space, but instead of discretely defining specific points in the space, we define probability distributions over the range we want to explore. These probability distributions can be anything from a uniform distribution (which says any point in the range is equally likely) to a Gaussian distribution (which has the familiar peak in the middle), or any other esoteric distributions, such as gamma or beta distributions. As long as we can sample from the distribution, we can use it for random search. Once we define the search space, we can sample points from the distribution and evaluate each of the points to find the best hyperparameter.

While the number of trials is a function of the defined search space for grid search, it is a user input for random search, so we get to decide how much time or computational budget we need to use for hyperparameter tuning and, because of that, we can also search over a larger search space.

With this new flexibility, let’s define a larger search space for our problem and use random search:

import scipy
from sklearn.model_selection import ParameterSampler
random_search_params = {
    # A uniform distribution between 10 and 100, but only integers
    "num_leaves": scipy.stats.randint(10,100),
    # A list of categorical string values
    "objective": ["regression", "regression_l1", "huber"],
    "random_state": [42],
    # List of floating point numbers between 0.3 and 1.0 with a resolution of 0.05
    "colsample_bytree": np.arange(0.3,1.0,0.05),
    # List of floating point numbers between 0 and 10 with a resolution of 0.1
    "lambda_l1":np.arange(0,10,0.1),
    # List of floating point numbers between 0 and 10 with a resolution of 0.1
    "lambda_l2":np.arange(0,10,0.1)
}
# Sampling from the search space number of iterations times
parameter_space = list(ParameterSampler(random_search_params, n_iter=27, random_state=42))

This also runs for about 15 minutes, but we have explored a larger search space. However, the best MAE reported was just 0.73752, which is lower than with grid search. Maybe if we run the search for a greater number of iterations, we will get a better score, but that is just a shot in the dark. Ironically, that is pretty much what random search also does. It closes its eyes and throws a dart at random places on the dartboard and hopes it hits the bull’s eye.

There are two terms in mathematical optimization called exploration and exploitation. Exploration ensures the optimization algorithm reaches different regions of the search space, whereas exploitation makes sure we search more in regions that are giving us better results. Random search is purely explorative and is unaware of what is happening as it evaluates different trials.

Let’s look at one last technique that tries to balance between exploration and exploitation.

Bayesian optimization

Bayesian optimization has a lot of similarities with random search. Both define their search space as probability distributions, and in both techniques, the user decides how many trials it needs to evaluate, but where they differ is the key advantage of Bayesian optimization. While random search is randomly sampling from the search space, Bayesian optimization is doing it intelligently. Bayesian optimization is aware of its past trials and the objective values that came out of those trials so that it can adapt future trials to exploit the regions where better objective values were seen. At a high level, it does this by building a probability model of the objective function and using it to focus trials on promising areas. The details of the algorithm are worth knowing and we have linked to a couple of resources in Further reading to help you along the way.

Now, let’s use a popular library, optuna, to implement Bayesian optimization for hyperparameter tuning on the GFM we have been training.

The process is quite simple. We need to define a function that takes in a parameter called trial. Inside the function, we sample the different parameters we want to tune from the trial object, train the model, evaluate the forecast, and return the metric we want to optimize (the MAE). Let’s quickly do that:

def objective(trial):
    params = {
        # Sample an integer between 10 and 100
        "num_leaves": trial.suggest_int("num_leaves", 10, 100),
        # Sample a categorical value from the list provided
        "objective": trial.suggest_categorical(
            "objective", ["regression", "regression_l1", "huber"]
        ),
        "random_state": [42],
        # Sample from a uniform distribution between 0.3 and 1.0
        "colsample_bytree": trial.suggest_uniform("colsample_bytree", 0.3, 1.0),
        # Sample from a uniform distribution between 0 and 10
        "lambda_l1": trial.suggest_uniform("lambda_l1", 0, 10),
        # Sample from a uniform distribution between 0 and 10
        "lambda_l2": trial.suggest_uniform("lambda_l2", 0, 10),
    }
    _model_config = ModelConfig(
        # Use the sampled params to initialize the model
        model=LGBMRegressor(**params, verbose=-1),
        name="Global Meta LightGBM Tuning",
        # LGBM is not sensitive to normalized data
        normalize=False,
        # LGBM can handle missing values
        fill_missing=False,
    )
    y_pred, feat_df = train_model(
        _model_config,
        _feat_config,
        missing_value_config,
        train_features,
        train_target,
        test_features,
        fit_kwargs=dict(categorical_feature=cat_features),
    )
    # Return the MAE metric as the value
    return ts_utils.mae(test_target["energy_consumption"], y_pred)

Once we have defined the objective function, we need to initialize a sampler. optuna has many samplers, such as GridSampler, RandomSampler, and TPESampler. For all standard use cases, TPESampler is the one to use. GridSampler does grid search and RandomSampler does random search. When defining a Tree Parzen Estimator (TPE) sampler, there are two parameters that we should pay attention to:

The rest of the parameters are best left untouched for the most common use cases.

Now, we create a study, which is the object that runs the trials and stores all the details about the trials:

# Create a study
study = optuna.create_study(direction="minimize", sampler=sampler)
# Start the optimization run
study.optimize(objective, n_trials=27, show_progress_bar=True)

Here, we define the direction of optimization, and we pass in the sampler we initialized earlier. Once the study is defined, we need to call the optimize method and pass the objective function we defined, the number of trials we need to run, and some other parameters. A full list of parameters for the optimize method is available here—https://optuna.readthedocs.io/en/stable/reference/generated/optuna.study.Study.html#optuna.study.Study.optimize.

This runs slightly longer, maybe because of the additional computation required to generate new trials, but still only takes about 20 minutes for the 27 trials. As expected, this has come up with another combination of hyperparameters for which the objective value is 0.72838 (the lowest before now).

To fully illustrate the difference between the three, let’s compare how the three techniques spent their computational budget:

Figure 10.7 – Distribution of computational effort (grid versus random versus Bayesian optimization)

Figure 10.7: Distribution of computational effort (grid versus random versus Bayesian optimization)

We can see that the Bayesian optimization has a fat tail on the lower side, indicating that it spent most of its computational budget evaluating and exploiting the optimal regions in the search space.

Let’s look at how the different trials with these techniques fared as the optimization procedure progressed.

The notebook has a more detailed comparison and commentary on the three techniques.

The bottom line is that if we have unlimited computation, grid search with a well-defined and fine-grained grid is the best option, but if we value the efficiency of our computational effort, we should go for Bayesian optimization.

Let’s see how the new parameters worked out for us:

Figure 10.8 – Aggregate metrics with the tuned GFM with meta-features

Figure 10.8: Aggregate metrics with the tuned GFM with meta-features

We have had huge improvements in MAE and meanMASE, mostly because we were optimizing for the MAE when hyperparameter tuning. The MAE and MSE have slightly different priorities and we will spend more time on that in Part 4, Mechanics of Forecasting. The runtime also increased because the new parameters build more leaves for a tree than before and are more complex than the default parameters.

Now, let’s look at another strategy for improving the performance of a GFM.

Partitioning

Out of all the strategies we have discussed so far, this is the most counter-intuitive, especially if you are coming from a standard machine learning or statistics background. Normally, we would expect the model to do well with more data, but partitioning or splitting the dataset into multiple, almost equal parts has been shown (empirically) to improve the accuracy of the model. While this has been seen empirically, why this happens is something that is still not quite clear. One explanation is that the GFMs have a slightly simpler job of learning when trained on a subset of similar entities and hence, can learn specific functions to subsets of similar entities. Another explanation for the phenomenon has been put forward by Montero-Manso and Hyndman (Reference 1). They put forward that partitioning the data is another form of increasing the complexity because instead of having log(|J|) as the complexity term, we have

$$log\left(\prod_{i=1}^{P} |H_i|\right)$$

where P is the number of partitions. With this rationale, the LFMs are special cases where P is equal to the number of time series in the dataset.

There are many ways we can partition the data, each with varying degrees of complexity.

Random partition

The simplest method is to randomly split the dataset into P-equal partitions and train separate models for each partition. This method faithfully follows the explanation that Montero-Manso and Hyndman provide because we are splitting the dataset randomly, with no concern for the similarity of the different households. Let’s see how we can do that:

# Define a function which splits a list into n partitions
def partition (list_in, n):
    random.shuffle(list_in)
    return [list_in[i::n] for i in range(n)]
# split the unique LCLids into partitions
partitions = partition(train_df.LCLid.cat.categories.tolist(), 3)

Then, we just loop over these partitions and train separate models for each partition. The exact code can be found in the notebook. Let’s see how well the random partition does:

Figure 10.9 – Aggregate metrics with the tuned GFM with meta-features and random partitioning

Figure 10.9: Aggregate metrics with the tuned GFM with meta-features and random partitioning

We can see a decrease in MAE and meanMASE even with a random partition. There is even a decrease in runtime because the individual models are working on less data and hence, train faster.

Now, let’s see another way of partitioning, keeping the similarity of different time series in mind.

Judgmental partitioning

Judgmental partitioning is when we use some attribute of the time series to split the dataset, and this is called judgmental because, usually, this depends on the judgment of the person who is working on the model. There are many ways of doing this. We can use some meta-feature, or we can use some characteristics of the time series (such as volume, variability, intermittency, or a combination of them) to partition the dataset.

Let’s use a meta-feature called Acorn_grouped to partition the dataset. Again, we will just loop over the unique values in Acorn_grouped and train a model for each value. We will also not use Acorn_grouped as a feature. The exact code is in the notebook. Let’s see how well this partitioning does:

Figure 10.10 – Aggregate metrics with the tuned GFM with meta-features and Acorn_grouped partitioning

Figure 10.10: Aggregate metrics with the tuned GFM with meta-features and Acorn_grouped partitioning

This does even better than random partitioning. We can assume each of the partitions (Affluent, Comfortable, and Adversity) has some kind of similarity, which makes the learning easier, and hence, we get better accuracy.

Now, let’s look at another way to partition the dataset, again, using similarity.

Algorithmic partitioning

In judgmental partitioning, we pick some meta-features or time series characteristics for partitioning the dataset. We pick a handful of dimensions to partition the dataset because we are doing it in our minds and our mental faculties cannot handle more than two or three dimensions well, but we can see this partitioning as an unsupervised clustering approach and this approach is called algorithmic partitioning.

There are two ways we can cluster time series:

tslearn is an open source Python library that has implemented a few time series clustering approaches based on the distances between time series. There is a link in Further reading for more information on the library and how it can be used for time series clustering.

In our example, we are going to use the first method, where we derive a few time series characteristics and use them for clustering. There are many features from statistical and temporal literature, such as autocorrelation, mean, variance, entropy, and peak-to-peak distance, that we can extract from the time series.

We can use another open source Python library called the Time Series Feature Extraction Library (tsfel) to make the process easier.

The library has many classes of features—statistical, temporal, and spectral domains—that we can choose from, and the rest is handled by the library. Let’s see how we can generate these features and create a dataframe to perform clustering:

import tsfel
cfg = tsfel.get_features_by_domain("statistical")
cfg = {**cfg, **tsfel.get_features_by_domain("temporal")}
uniq_ids = train_df.LCLid.cat.categories
stat_df = []
for id_ in tqdm(uniq_ids, desc="Calculating features for all households"):
    ts = train_df.loc[train_df.LCLid==id_, "energy_consumption"]
    res = tsfel.time_series_features_extractor(cfg, ts, verbose=False)
    res['LCLid'] = id_
    stat_df.append(res)
stat_df = pd.concat(stat_df).set_index("LCLid")

The dataframe looks something like this:

Figure 10.11 – Features extracted from different time series

Figure 10.11: Features extracted from different time series

Now that we have the dataframe with each row representing a time series with different features, we can ideally apply any clustering method, such as k-means, k-medoids, or HDBSCAN, and find clusters. However, in high dimensions, a lot of the distance metrics (including Euclidean) do not work as well as they are supposed to. There is a seminal paper on the topic by Charu C. Agarwal et al. from 2001 that explores the topic. When we increase the dimensionality of the space, our common sense (which conceptualizes three dimensions) does not work as well and, as a consequence, common distance metrics such as Euclidean distance do not work very well with high dimensions. We have linked to a blog summarizing the paper (in Further reading) and the paper itself (Reference 5), which make the concept clearer. So, a common way of handling high-dimensional clustering is by performing dimensionality reduction first and then using normal clustering.

Principal Component Analysis (PCA) was the go-to tool in the field, but since PCA only captures and details linear relationships while reducing the dimensions, nowadays, another class of techniques is starting to become more popular—manifold learning.

t-distributed Stochastic Neighbor Embeddings (t-SNE) is a popular technique from this category, which is really popular for high-dimensional visualization. It is a really clever technique where we project the points from a high-dimensional space to a lower dimension, keeping the distribution of distance in the original space as close as possible to the one in lower dimensions. There is a lot to learn here that is beyond the scope of this book. There are links in the Further reading section that can help you get started.

To cut a long story short, we will be using t-SNE to reduce the dimensions of the dataset we have and then cluster the dataset with the reduced dimensions. If you really want to cluster time series and use those clusters in some other way, I would not suggest using t-SNE because it doesn’t preserve the distance between points and the density of points. The distil.pub article in Further reading throws more light on the issue. But in our case, we are using the clusters just as a grouping for training another model, so this approximation can do well. Let’s see how we do that:

from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from src.utils.data_utils import replace_array_in_dataframe
from sklearn.manifold import TSNE #T-Distributed Stochastic Neighbor Embedding
# Standardizing to make distance calculation fair
X_std = replace_array_in_dataframe(stat_df, StandardScaler().fit_transform(stat_df))
#Non-Linear Dimensionality Reduction
tsne = TSNE(n_components=2, perplexity=50, learning_rate="auto", init="pca", random_state=42, metric="cosine", square_distances=True)
X_tsne = tsne.fit_transform(X_std.values)
# Clustering reduced dimensions into 3 clusters
kmeans = KMeans(n_clusters=3, random_state=42).fit(X_tsne)
cluster_df = pd.Series(kmeans.labels_, index=X_std.index)

Since we reduced the dimensions to two, we can also visualize the clusters formed:

Figure 10.12 – Clustered time series after t-SNE dimensionality reduction

Figure 10.12: Clustered time series after t-SNE dimensionality reduction

We have three well-defined clusters formed and now we are just going to use these clusters to train a model for each cluster. As usual, we loop over the three clusters and train the models. Let’s see how we did so:

Figure 10.13 – Aggregate metrics with the tuned GFM with meta-features and clustered partitioning

Figure 10.13: Aggregate metrics with the tuned GFM with meta-features and clustered partitioning

It looks as though this is the best MAE we have seen in all our experiments, but the three partition techniques have very similar MAEs. We can’t see whether any one is better than the other just by looking at a single hold-out set. For good measure, we can run these forecasts with a test dataset using the 01a-Global_Forecasting_Models-ML-test.ipynb notebook in the Chapter08 folder. Let’s see how the aggregate metrics are on the test dataset:

Figure 10.14 – Aggregate metrics on test data

Figure 10.14: Aggregate metrics on test data

As expected, the clustered partition is still the methodology that performs the best in this case.

In Chapter 8, Forecasting Time Series with Machine Learning Models, it took us 8 minutes and 20 seconds to train an LFM for all the households in our dataset. Now, with the GFM paradigm, we finished training a model in 57 seconds (in the worst-case scenario). That’s 777% less training time and this comes with an 8.78% decrease in the MAE.

We chose to do these experiments with LightGBM. This does not mean that LightGBM or any other gradient-boosting model is the only choice for GFMs, but they are a pretty good default. A well-tuned gradient-boosted trees model is a very difficult baseline to beat, but as always in machine learning, we should check what works best using well-defined experiments.

Although there are no hard and fast rules or cutoffs for when a GFM makes more sense than an LFM, as the number of time series in a dataset increases, the GFM becomes more favorable, both from the perspective of accuracy and computation.

Although we have achieved good results using GFMs, typically the complex models that do well in this paradigm are black boxes. Let’s look at some ways to open the black box and understand and explain the model better.

Interpretability

Interpretability can be defined as the degree to which a human can understand the cause of a decision. In machine learning and artificial intelligence, that translates to the degree to which someone can understand the how and why of an algorithm and its predictions. There are two ways to look at interpretability—transparency and post hoc interpretation.

Transparency is when the model is inherently simple and can be simulated or thought about using human cognition. A human should be able to fully understand the inputs and the process a model takes to convert these inputs to outputs. This is a very stringent condition that almost none of the model machine learning or deep learning models satisfy.

This is where post hoc interpretation techniques shine. There is a wide variety of techniques that use the inputs and outputs of a model to understand why a model has made the predictions it has.

There are many popular techniques such as permutation feature importance, Shapley values, and LIME. All of these are general-purpose interpretation techniques that can be used on any machine learning model and that includes the GFMs we were discussing. Let’s talk about a few of them at a high level.

Mean decrease in impurity:

This is the regular “feature importance” that we get out of the box from tree-based models. This technique measures how much a feature reduces impurity (such as Gini impurity in classification or variance in regression) when used to split nodes in a decision tree. The higher the reduction in impurity, the more important the feature is considered. However, it’s biased towards continuous features or those with high cardinality. It is fast and readily available in libraries like scikit-learn but may give misleading results if features have varying scales or many categories.

Drop column importance (Leave One Covariate Out (LOCO)):

This method assesses feature importance by iteratively removing one feature at a time and retraining the model. The drop in performance from the baseline model indicates the importance of that feature. It is model-agnostic and captures interactions between features, but is computationally expensive since it requires retraining the model for each feature removed. It can also give misleading results if collinear features exist, as the model may compensate for the removed feature.

Permutation importance:

Permutation importance measures the drop in model performance when the values of a single feature are randomly shuffled, disrupting its relationship with the target. This technique is intuitive and model-agnostic, and it doesn’t require retraining the model, making it computationally efficient. However, it can inflate the importance of correlated features, as models can rely on related features to compensate for the permuted one.

Partial Dependence Plots (PDPs) and Individual Conditional Expectation (ICE) plots:

PDPs visualize the average effect of a feature on the model’s predictions, showing how the target variable changes as the feature’s values change, while ICE plots show the effect of a feature for individual instances. These plots help understand the feature-target relationship but assume independence between features, which can lead to misleading interpretations in the presence of correlated variables.

Local Interpretable Model-agnostic Explanations (LIME):

LIME is a model-agnostic technique that explains individual predictions by approximating a complex model locally using simpler, interpretable models, like linear regression. It works by generating perturbations of the data point in question and fitting a local model to these samples. This method is intuitive and widely applicable to both structured and unstructured data (text and images), but defining the right locality for perturbations can be challenging, especially for tabular data.

SHapley Additive exPlanations (SHAP):

SHAP unifies several interpretation methods, including Shapley values and LIME, into a single framework that attributes feature importance in a model-agnostic way. SHAP provides both local and global interpretations and benefits from fast implementations for tree-based models (TreeSHAP). It combines the theoretical strength of Shapley values with practical efficiency, although it can still be computationally intensive for large datasets.

Each technique has its strengths and trade-offs, but SHAP stands out due to its strong theoretical foundation and ability to connect local and global interpretations effectively. For more extensive coverage of such techniques, I have included a few links in Further reading. The blog series by yours truly and the free book by Christoper Molnar are excellent resources for you to get up to speed (more on interpretability in Chapter 17).

Congratulations on finishing the second part of the book! It has been quite an intensive part where we went over quite a bit of theory and practical lessons, and we hope you are now comfortable with using machine learning for time series forecasting.

Summary

To round up the second part of the book nicely, we explored GFMs in detail and saw why they are important and why they are an exciting new direction in time series forecasting. We saw how we can use a GFM using machine learning models and also reviewed many techniques to make GFMs perform better, most of which are quite frequently used in competitions and industry use cases alike. We also took a high-level look at the interpretability techniques. Now that we have wrapped up the machine learning section of the book, we will move on to a specific type of machine learning that has become well-known over the past few years—deep learning—in the next chapter.

References

The following are sources that we have referenced throughout the chapter:

  1. Montero-Manso, P., Hyndman, R.J. (2020), Principles and algorithms for forecasting groups of time series: Locality and globality. arXiv:2008.00444[cs.LG]: https://arxiv.org/abs/2008.00444.
  2. Micci-Barreca, D. (2001), A preprocessing scheme for high-cardinality categorical attributes in classification and prediction problems. SIGKDD Explor. Newsl. 3, 1 (July 2001), 27–32: https://doi.org/10.1145/507533.507538.
  3. Fisher, W. D. (1958). On Grouping for Maximum Homogeneity. Journal of the American Statistical Association, 53(284), 789–798: https://doi.org/10.2307/2281952.
  4. Fisher, W.D. (1958), A preprocessing scheme for high-cardinality categorical attributes in classification and prediction problems. SIGKDD Explor. Newsl. 3, 1 (July 2001), 27–32.
  5. Aggarwal, C. C., Hinneburg, A., and Keim, D. A. (2001). On the Surprising Behavior of Distance Metrics in High Dimensional Spaces. In Proceedings of the 8th International Conference on Database Theory (ICDT ‘01). Springer-Verlag, Berlin, Heidelberg, 420–434: https://dl.acm.org/doi/10.5555/645504.656414.
  6. Oreshkin, B. N., Carpov D., Chapados N., and Bengio Y. (2020). N-BEATS: Neural basis expansion analysis for interpretable time series forecasting. 8th International Conference on Learning Representations, ICLR 2020: https://openreview.net/forum?id=r1ecqn4YwB.

Further reading

The following are a few resources that you can explore for a detailed study:

Leave a Review!

Thank you for purchasing this book from Packt Publishing—we hope you enjoy it! Your feedback is invaluable and helps us improve and grow. Once you’ve completed reading it, please take a moment to leave an Amazon review; it will only take a minute, but it makes a big difference for readers like you.

Scan the QR or visit the link to receive a free ebook of your choice.

https://packt.link/NzOWQ