17

Probabilistic Forecasting and More

Throughout the book, we have learned different techniques to generate a forecast, including some classical methods, using machine learning, and some deep learning architectures as well. But we have been focusing on one typical type of forecasting problem—generating a point forecast for continuous time series with no hierarchy and a good amount of history. We have been doing that because this is the most popular kind of problem you will face. But in this chapter, we will take some time to look at a few niche topics that, although less popular, are no less important.

In this chapter, we will focus on these topics:

• Probabilistic forecasting
  • Probability Density Functions
  • Quantile functions
  • Monte Carlo Dropout
  • Conformal Prediction
• Intermittent/sparse time series forecasting
• Interpretability
• Cold-start forecasting
  • Pre-trained models like TimeGPT
  • Similarity-based forecasting
• Hierarchical forecasting

Probabilistic forecasting

So far, we have been talking about the forecast as a single number. We have been projecting our DL models to a single dimension or training our machine learning models to output a single number. Subsequently, we were training the model using a loss, such as mean squared loss. This paradigm is what we call a point forecast. But we are not considering one important aspect. We are using the history to train our model to make the best guess. But how sure is the model about its prediction? Those of you who are aware of machine learning and classification problems would recognize that for classification problems, besides getting a prediction of which class the sample belongs to, we also get a notion of the uncertainty of the model. But our forecasting is a regression problem and we don’t get the uncertainty for free.

But why is quantifying uncertainty important in forecasting? Any forecast is created for some purpose, some downstream task for which the forecasted information is being used. In other words, there is some decision that has to be taken using the forecast we generate. And when making a decision, we usually would like to have the maximum amount of information available to us.

Let’s look at an example to really drive home the point. You have recorded your monthly grocery consumption for the last 5 years and, using the techniques in the book, created a super accurate forecast and an app that tells you how much to shop for in any month. You open up the app and it tells you that you need to buy two bread loaves this month. You head out to the supermarket, bag two loaves of bread, and return home. And a week before the end of the month, the bread loaves ran out, and you starved the rest of the time. At that peak of starvation, you started questioning your decisions and the forecast. You analyzed the data to figure out where you went wrong and realized that your consumption of bread per month varied a lot. In some months you consumed 4 loaves and, in some other months, it was just 1. So, there is a good chance that this forecast will leave you with no bread for a few months and with excess bread for a few other months.

Then, you read this chapter and convert your forecast to a probabilistic forecast and now it tells you that 50% of the time your next month’s consumption of bread is 2 loaves. But now, there is an additional feature in the app, which asks you whether you prefer to starve or have excess bread at the end of the day. So, depending on your appetite to starve or save money, you decide on an option. Let’s say you don’t want to starve, but are okay if bread runs out like 10% of the time. Once you enter this preference into the app, it revises its forecast and tells you that you should get 3 loaves of bread, and you never starve again (also because you smartened up and bought other stuff from the supermarket).

The ability to use the uncertainty in the forecast and revise it according to our appetite for risk is one of the main utilities for probabilistic forecasting. It also helps our forecast to be more transparent and trustworthy to the users.

Now, let’s quickly look at the types of uncertainty one would have in a prediction problem using learned models.

Types of Predictive Uncertainty

We saw in Chapter 5 that supervised machine learning is nothing but learning a function, $\hat{y} = h(X, \phi)$, where h, along with $\phi$, is the model that we learn and is the input data. So, if we think about the sources of uncertainty, it can be from two of these components.

The model, h, we learned is an approximation using a dataset, X, which may or may not cover all the cases completely and some uncertainty can be introduced to the system from this. We call this Epistemic Uncertainty. In the context of machine learning, epistemic uncertainty can occur when the model has not been exposed to enough data, the model itself is insufficient to learn the complexity of the problem, or the data it has been trained on does not represent all possible scenarios. This is also known as systemic or reducible uncertainty because this is the portion of the total predictive uncertainty that can be actively reduced by having better models, better data, and so on; in other words, by gaining more knowledge about the system. Let’s see a few examples to make the concept clearer:

• If a weather forecast model has less data from a particular region (maybe because of a faulty sensor), there is going to be less knowledge about that region, and this increases uncertainty.
• If a linear model was used for a non-linear problem, we would be introducing some uncertainty because of the simpler model having less knowledge about the system.
• If an economic forecasting model is not trained using a few key influencing factors like change in economic policies, climate change influencing the decisions, and so on, this also creates some uncertainty.

The good thing about this kind of uncertainty is that this is something that we can actively reduce by collecting better data, training better models, and so on.

Now there is another kind of uncertainty in the total predictive uncertainty—Aleatoric Uncertainty. This refers to the inherent randomness in the data which cannot be explained away. This is also known as statistical or irreducible uncertainty. Although our universe appears deterministic to us, there is an ever-prevalent layer of uncertainty just beneath the surface.

For example, the motion of the celestial objects can be calculated accurately (thanks to General Relativity and Einstein), but still, a random asteroid can hit any of the bodies and cause alterations to the calculated trajectory. This kind of irreducible and unavoidable uncertainty is referred to as aleatoric uncertainty. Let’s see a few examples:

• Weather predictions, no matter how accurate the measurements and the models are, are still variable. There is an inherent randomness in weather that we might never be able to completely explain.
• The performance of an athlete, no matter how much they have trained and followed the rules, is still not completely deterministic. There are a lot of factors, like weather, health, and other random events during or before the game, which will affect the performance of the athlete.

Now that we understand the different types of uncertainty, and why we need uncertainty quantification, let’s see what it means in the context of forecasting.

What are probabilistic forecasts and Prediction Intervals?

A probabilistic forecast is when the forecast, instead of having a single-point prediction, captures the uncertainty of that forecast as well. Probabilistic forecasting is a method of predicting future events or outcomes by providing a range of possible values along with associated probabilities or confidence levels. This approach captures the uncertainty in the prediction process.

In the econometrics and classical time series world, the prediction intervals were already baked into the formulation. The statistical grounding and strong assumptions of those methods made sure that the output of those models was readily interpreted in a probabilistic way as well (so long as you could satisfy the assumptions that were stipulated by those models). But in the modern machine learning/deep learning world, probabilistic forecasting is not an afterthought. A combination of factors, such as fewer rigid assumptions and the way we train the models, leads to this predicament.

There are different methods with which we can add the probabilistic dimension to our forecast, and we will go through a few of them in this chapter. But before that, let’s also understand one of the most useful manifestations of probabilistic forecasts—Prediction Intervals.

A Prediction Interval is a range within which a future observation is expected to fall with a specified probability. For instance, if we have a 95% prediction interval for a time step of [5,8], we say that 95% of the time, the actual value will lie between 5 and 8. Let’s take an example of a normal distribution with mean, $\mu_t$, and variance, $\sigma_t^2$ as the forecast at timestep, t (one of the techniques we will talk about gives us just that). So, the prediction interval at time, t, with a significance level (the probability with which the forecast can fall outside the interval), $\alpha$, can be written as:

$$[\mu_t - z\sigma, \mu_t + z\sigma]$$

where z is the z-score of the normal distribution.

There are a few terms and concepts we need to clarify before we begin the discussion in detail.

Confidence levels, error rates, and quantiles

When working with prediction intervals, it’s crucial to understand the relationship between confidence levels, error rates, and quantiles. These concepts help in determining the range within which future observations are expected to fall with a certain probability.

Error Rate ($\alpha$) is the allowable probability that the prediction interval will not contain the future observation. It is typically expressed in percentages or as a decimal between 0 and 1. If we say $\alpha = 10\%$ or $\alpha = 0.1$, it means that there is a 10% chance that the future observation will not be in the prediction interval.

Confidence Level ($CL$) is the complement of error rate and is the probability that the prediction interval contains the future observation. $CL = 1 - \alpha$. If we say that the error rate is 10%, the confidence level would be 90%.

Quantiles are points that divide the data into intervals with equal probabilities. In simpler terms, a quantile shows the value below which a certain percentage of data falls. For instance, the 5th percentile or 0.05th quantile marks the point where 5% of data lies below it. Therefore, we can use quantiles as well to define the prediction intervals where we don’t have an analytical way to get the prediction intervals based on distributional assumption.

In Figure 17.1, we show the prediction intervals for a standard normal distribution.

Figure 17.1: Prediction intervals for a standard normal distribution

There is a strong link between error rates, confidence levels, and quantiles. Error rates and confidence levels are direct complements to each other and can be used interchangeably to define the confidence that we want our prediction intervals to have or the error rate we are okay with from the prediction intervals. Another way to look at it is by the area under the curve. In Figure 17.1, the area of the green-shaded region denotes the confidence level, and the area of the red area denotes the error rate.

For a standard normal distribution, we can directly get the prediction intervals by using the analytical formula:

$$[\mu - z_{\alpha/2} \cdot \sigma, \mu + z_{\alpha/2} \cdot \sigma]$$

where $\mu$ is the mean of the distribution, $\sigma$ is the standard deviation of the distribution, and $z_{\alpha/2}$ is the critical value from the standard normal distribution corresponding to the desired confidence level, $1 - \alpha$. $\alpha / 2$ is taken because we are allowing for the error rate to be spread on both sides (red shaded area on both sides of the curve in Figure 17.1).

Now let’s look at how error rates and confidence levels are linked with quantiles because if we don’t know what the distribution is (and we don’t want to assume any distribution), we can’t go by the analytical formula to get prediction intervals. In such cases, we can use quantiles to get the same. Just like we did with the analytical formula, the error rate, $\alpha$, should be equally divided across both sides. And therefore, the prediction interval would be:

$$[q_{\alpha/2}, q_{1-(\alpha/2)}]$$

where $q_t$ is the t^th quantile. So, from the definition of quantile, we know $q_{\alpha/2}$ would have $\alpha / 2$ % of data below it and $q_{1-(\alpha/2)}$ would have $\alpha / 2$ % of data above it, thus making the area outside the intervals to be $\alpha$.

Using this relation, we can go from error rates to quantiles or confidence levels to quantiles. If the error rate is $\alpha$ we have already seen what the corresponding quantiles denoting the prediction interval are. Let’s also see one more quick formula to go from confidence levels (declared as percentages) to quantiles:

$$[q_{50-(CL/2)}, q_{50+(CL/2)}]$$

In Python code, this is simply:

level = 95  # Confidence levels
qs = [50 - level / 2, 50 + level / 2] # Quantiles

Now, let’s look at how to measure the goodness of prediction intervals.

Measuring the goodness of prediction intervals

We know what a probabilistic forecast is and what prediction intervals are. But before we look at techniques to generate prediction intervals, we need a way to measure the goodness of such an interval. Standard metrics like Mean Absolute Error or Mean Squared Error hold no more because they are point forecast measuring metrics.

What do we want from a prediction interval? If we have a prediction interval with 90% confidence, we would expect the data points to lie between the interval at least 90% of the time. This can easily be obtained by having very wide prediction intervals, but that again becomes a useless prediction interval. So, we want our prediction interval to be as narrow as possible and still have the 90% confidence criteria respected. To measure these two distinct aspects, we can use two metrics—Coverage and Average Length of Prediction Intervals.

Coverage is the proportion of true values that fall within the prediction intervals. Mathematically, if we denote the Prediction Interval by $[L_i, U_i]$ for each observation i, and true value by $y_i$, coverage can be defined as:

$$Coverage = \frac{1}{N} \sum_{i=1}^{N} \mathbb{I}(L_i \leq y_i \leq U_i)$$

where $\mathbb{I}$ is an indicator function that equals 1 if the condition inside is true and 0 otherwise, and N is the total number of observations. A coverage metric close to the desired confidence level (e.g., 95% for a 95% prediction interval) indicates that the model’s uncertainty estimates are well-calibrated.

Average Length of Prediction Intervals is calculated by averaging the lengths of the prediction intervals across all observations. Using the same notations as above, it can be mathematically written as:

$$Average\ Length = \frac{1}{N} \sum_{i=1}^{N} (U_i - L_i)$$

This metric helps in understanding the trade-off between the coverage of the intervals and their precision. Let’s also look at Python functions for both of these metrics:

Coverage:

import numpy as np
def coverage(y_true, lower_bounds, upper_bounds):
    """
    Calculate the coverage of prediction intervals.
   
    Parameters:
    y_true (array-like): True values.
    lower_bounds (array-like): Lower bounds of prediction intervals.
    upper_bounds (array-like): Upper bounds of prediction intervals.
   
    Returns:
    float: Coverage metric.
    """
    y_true = np.array(y_true)
    lower_bounds = np.array(lower_bounds)
    upper_bounds = np.array(upper_bounds)
    # Check if true values fall within the prediction intervals
    coverage = np.mean((y_true >= lower_bounds) & (y_true <= upper_bounds))
    return coverage

Average Length:

def average_length(lower_bounds, upper_bounds):
    """
    Calculate the average length of prediction intervals.
   
    Parameters:
    lower_bounds (array-like): Lower bounds of prediction intervals.
    upper_bounds (array-like): Upper bounds of prediction intervals.
   
    Returns:
    float: Average length of prediction intervals.
    """
    lower_bounds = np.array(lower_bounds)
    upper_bounds = np.array(upper_bounds)
   
    # Calculate the length of each prediction interval
    lengths = upper_bounds - lower_bounds
    # Calculate the average length
    average_length = np.mean(lengths)
    return average_length

Both of these Python functions can be found in src/utils/ts_utils.py and we will be using them to measure the quality of prediction intervals generated in the chapter.

Now, let’s look at different techniques we can use to get probabilistic forecasts and how we can use them practically.

Probability Density Function (PDF)

This is one of the most common techniques for probabilistic forecasting, especially in the deep learning space, because of the ease of implementation. The forecast at time t, $\hat{y}_i$, can be seen as the realization of a probability distribution, $p(\hat{y}_t)$. And instead of estimating $\hat{y}_i$, we can estimate $p(\hat{y}_t)$. If we assume $p(\hat{y}_t)$ is one of the parametrized distributions, $\mathcal{D}(\theta_t)$, with parameters $\theta_t$, then we can estimate the parameters $\theta_t$, instead of $\hat{y}_i$, directly.

For instance, if we assume the forecast is drawn from normal distribution, we can model

$$\hat{y}_t \sim \mathcal{N}(\mu_t, \sigma_t)$$

So, instead of getting our model to output $\hat{y}_i$, we can get it to output $\mu_t$ and $\sigma_t$. And with $\mu_t$ and $\sigma_t$, we can easily calculate the prediction intervals at the given $\alpha$:

$$Lower\ Bound = \mu - Z \cdot \sigma$$

$$Upper\ Bound = \mu + Z \cdot \sigma$$

where Z is the critical value from the standard normal distribution corresponding to the desired confidence level. For a 90% confidence level, $Z \approx 1.645$. Simple enough, right? Not so fast!

Now that we are modeling the parameters of a distribution, how do we train the model? We still have the actual point forecast as the target. In the normal distribution case, the targets are still actual $y_t$ and not the means and standard deviations. We get over this problem by using a loss function like log likelihood, instead of losses like Squared Error.

For instance, let’s say we have a set of i.i.d observations (in our case, the target), $y_1, y_2, \ldots, y_n$. With the predicted parameters of the assumed distribution ($D(\theta)$), we will be able to calculate the probability of each of the target, $p(y_t)$. The i.i.d assumption means each sample is independent of each other. And high-school mathematics tells us that when two independent events happen, we can calculate their joint probability by multiplying the two independent probabilities together. Using the same logic, we can calculate the joint probability or likelihood of all n i.i.d observations (probability that all these events occur) by just multiplying all the individual probabilities together.

$$Likelihood = \prod_{i=1}^{n} p(y_i)$$

Maximizing the likelihood helps the model learn the right parameters for each sample such that the probability under the assumed distribution maximizes. We can intuitively think about this as follows. For an assumed distribution like normal distribution, maximizing the likelihood makes sure that the target falls in the center of the distribution defined by the predicted parameters for each sample.

However, this operation is not numerically stable. Since probabilities are $0 \lt p \lt 1$, multiplying them together makes the result progressively smaller and can soon lead to numerical underflow issues. Therefore, we use the log likelihood, which is nothing but the likelihood but log transformed. We do it because:

• Being a strictly monotonic transformation, optimizing a function is equivalent to optimizing the log transform of the function. Therefore, optimizing likelihood and log likelihood is the same thing.
• Log transformation converts the multiplication into an addition, which is a much more numerically stable operation.

$$Log\ Likelihood = \sum_{i=1}^{n} \log(p(y_t))$$

The major disadvantage of this approach is that this relies on parametrized probability distributions whose log likelihood computation is tractable. Therefore, we are forced to make assumptions about the output and pick a distribution that might fit ahead of time.

This is a double-edged sword. On one hand, we can inject some domain knowledge into the problem and regularize the model training, but on the other hand, if we aren’t clear on whether choosing a normal distribution is the right choice, it can lead to an unnecessarily constrained model.

Many popular distributions, such as Normal, Poisson, Negative Binomial, Exponential, LogNormal, Tweedie, and so on, can be used for generating probabilistic forecasts with this technique.

Now, we have all the components for training and learning a model and can make it predict a full probabilistic distribution instead of a point forecast. With all that theory set, let’s switch gears and see how we can use this technique.

Forecasting with PDF—machine learning models

We have seen how to use standard machine learning models for forecasting in Part 2, Machine Learning for Time Series. But all of that was for point forecast. Can we easily convert all of that into probabilistic forecasts using the PDF approach? In theory, yes. But practically, it’s not that easy. All the popular implementations of machine learning models like sci-kit learn, xgboost, lightgbm, and so on take a point prediction paradigm. And as users of such open-source libraries, it isn’t easy for us to tweak and re-write the code to make it optimize the log likelihood as a regression loss. But fear not, it is not impossible. NGBoost is a distant cousin of the very popular gradient boosting models, like xgboost and lightgbm, and it is implemented such that it predicts the PDF, instead of the point prediction.

We aren’t going too deep into what NGBoost is and how it differs from the regular gradient boosting models here, but just know that it is a model that predicts probability distributions instead of point predictions.

Further reading has a link to a blog about NGBoost, which goes into a bit more depth on what the model is.

Let’s use a sample of eight time series from the M4 competition (which has 100,000 time series, references 5) for the probabilistic forecasts. The data is easily available online and the download script is included in the notebook. We are using this simpler dataset than the one we have been working on because I want to avoid complicating the narrative with exogenous variables and such. Here is the plot of the last 100 time steps of the eight sampled time series.

Figure 17.2: Last 100 timesteps of 8 sampled time series from the M4 Competition. Test period is drawn in a dotted purple line.

Let’s use mlforecast to create some features quickly in order to convert this into a regression problem. We had a bonus notebook back in Chapter 6, which showed how to use mlforecast as an alternative for the feature engineering that is included in the book’s repository. For the detailed code, refer to the full notebook, but for now, let’s assume that we have a dataframe called data, which has all the features necessary to run a machine learning model. We have split it into train and val and then subsequently to X_train, y_train, X_val, and y_val. Now, let’s see how we can train a model, assuming the output is a Normal distribution.

from ngboost import NGBRegressor
from ngboost.distns import Normal
# Training the model
ngb = NGBRegressor(Dist=Normal).fit(X_train, Y_train)

NGBoost doesn’t have a lot of parameters to tune and therefore isn’t as flexible as other gradient-boosting decision trees (GBDT). And it’s also not as fast as the other GBDTs. This is a model that you use only for special use cases when you need probabilistic outputs.

These are a few parameters that NGBoost has:

• Dist: This is the assumed distributional form of the output. The package currently supports Normal, LogNormal, and Exponential—all of which can be imported from ngboost.distns
• Score: This is any valid scoring function that is used to compare the predicted distributions to observations. The log likelihood score that we discussed earlier is called LogScore in the package and is the default value. All scoring functions can be imported from ngboost.scores.
• n_estimators: This is the number of estimators that are used in the boosted trees.
• learning_rate: This is the learning rate used for combining the boosted trees.
• mini_batch_frac: The percent of rows which is sub-sampled for each iteration. This is set to 1.0 as the default value.

Now that we have a trained NGBoost model, let’s see how we can use it to generate predictions and prediction intervals.

To get point predictions, the syntax is exactly the same as the sci-kit learn API.

y_pred = ngb.predict(X_val)

This is nothing but a wrapper method that calculates the location parameter of the assumed distribution. For instance, for the Normal distribution, the mean of the predicted distribution is the point prediction.

Now to get the underlying probabilistic prediction and subsequently the prediction intervals, we need to use a different method:

y_pred_dists = ngb.pred_dist(X_val)

Each point in y_pred_dists is a complete distribution. If we want to instead the first five predicted points, we can do the following:

y_pred_dists[0:5].params

Now, to get the prediction interval, we can use y_pred_dist and call a method, giving the level of confidence we expect. This, in turn, calls the scipy distribution (like scipy.stats.norm), which has a method, interval, to get the intervals given the confidence level.

y_pred_lower, y_pred_upper = y_pred_dists.dist.interval(0.95)

Now, we have a window around y_pred wide enough to envelope the uncertainty expected at each data point—the prediction interval. Let’s look at the forecast and the metrics.

Figure 17.3: Forecast with prediction intervals from NGBoost

And below are the metrics we calculated for these eight time series (Mean Absolute Error, Coverage, and Average Length):

Figure 17.4: Metrics for NGBoost

Although for some time series, the intervals seem to be good, some others (like time series H103) seem to have way too narrow a prediction interval, which is also evident in the low coverage as well.

Forecasting with PDF—deep learning models

Unlike machine learning models, it’s very easy to convert all the deep learning models we have learned about in the book to their PDF version. Remember the discussion we had before we started the practical application? The major changes we needed to do were these:

• Instead of predicting point forecast (single number), we predict the parameters of a probability distribution (one or more numbers).
• Instead of using a point loss like Mean Squared Error, use a probabilistic scoring function like log likelihood.

In the deep learning paradigm, these are very simple changes to make, aren’t they? In almost all the deep learning models we have learned to use in the book, there is a linear projection at the end, which projects the output into the required dimensions. Conceptually, it’s simple enough to make these output multiple numbers (parameters of the assumed distribution) by changing the linear projection. Similarly, it’s simple enough to change the loss function as well. Notably, DeepAR (Reference 3) is a well-known deep learning model that uses this technique for probabilistic forecasting.

Let’s see how we can do this in neuralforecast (the library we were using Chapter 16). In our example here, we will take a simple model like an LSTM, but we can do the same with any model because all we are doing is switching the loss function to DistributionLoss.

And, for this example, we are going to use the M4 competition dataset, which is freely available (the code to download the dataset is included in the notebook).

Let’s start from the point where we have the data formatted the way neuralforecast expects in Y_train_df and Y_test_df. The first thing we need to do is import the necessary classes.

from neuralforecast import NeuralForecast
from neuralforecast.models import LSTM
from neuralforecast.losses.pytorch import DistributionLoss

The only class that we haven’t looked at before in the book is DistributionLoss. This is a class that wraps torch.distribution classes and implements the negative log likelihood loss we discussed earlier. At the time of writing this book, the DistributionLoss class supports these underlying distributions:

• Poisson
• Normal
• StudentT
• NegativeBinomial
• Tweedie
• Bernoulli (Temporal Classifiers)
• ISQF (Incremental Spline Quantile Function)

The choice between these different distributions is totally up to the modeler and is a key assumption in the model. If the output we are modeling is expected to be in a normal distribution, then we can choose Normal. These are the major parameters of DistributionLoss:

• distribution: This is a string that identifies which distribution we are assuming. It can be any one from the list we saw before.
• level: This is a list of floats that defines the different confidence levels we are interested in modeling. For instance, if we want to model 80% and 90% confidence, we should give the values as [80,90].
• quantiles: This is an alternate way of defining the level. Instead of 95% confidence, you can define them in quantiles → [0.1, 0.9].

  Practitioner’s tip:

  Making a distributional assumption requires a deep study of the domain and data. But if you are not completely familiar with these distributions, Normal or StudentT is a good starting point as a lot of data resembles normal distribution. But before you jump into using Normal distributions everywhere, you should do a bit of literature study in the domain of the problem and choose the distribution accordingly. For instance, intermittent demand or sparse demand, which is very common in retail, is better modeled using a Poisson distribution. If the forecast is a count data (positive integers), Negative Binomial is a good option.

Now let’s set a horizon, the levels we need, and a few hyperparameters for LSTM (we’ve chosen some small and simple hyperparameters to make the training faster. In real-world problems, it is advisable to do a hyperparameter search to find the best parameters).

horizon = 48
levels = [80, 90]
lstm_config = dict(input_size=3*horizon, encoder_hidden_size=8, decoder_hidden_size=8)

Now we need to define the models we are going to use and the NeuralForecast class. Let’s define two models—one using Normal and another using StudentT.

models = [
    LSTM(
        h=horizon,
        loss=DistributionLoss(distribution="StudentT", level=levels),
        alias="LSTM_StudentT",
        **lstm_config
    ),
    LSTM(
        h=horizon,
        loss=DistributionLoss(distribution="Normal", level=levels),
        alias="LSTM_Normal",
        **lstm_config
    ),
]
# Setting freq=1 because the ds column is not date, but instead a sequentially increasing number
nf = NeuralForecast(models=models, freq=1)

Notice how the syntax is exactly the same as for point forecast, except for the Distribution Loss we chose. Now all that’s left is to train the models.

nf.fit(df=Y_train_df)

Once the model is trained, we can predict using the predict method. This output will have the point forecast under the alias we have defined and the high and low intervals for all the levels we have defined.

Y_hat_df = nf.predict()

Now that we have the probabilistic forecasts, let’s take a look at them and also calculate the metrics.

Figure 17.5: Forecast with prediction intervals from LSTM with StudentT distribution as output

Figure 17.6: Metrics for the LSTM with Normal and StudentT distribution outputs

If we compare the coverages with the ones with NGBoost, we can see that the deep learning approach increased the coverage, but also has wider than necessary intervals in most cases (as evidenced by larger average widths).

The biggest disadvantage of this method is that we restrict the output to one of the parametrized distributions. In many real-world cases, the data might not conform to any parametric distributions.

Now, let’s see a method that doesn’t require the assumption of any parametric distribution, but still gets the prediction intervals.

Quantile function

If our only aim is to get prediction intervals, we can also do the same using quantiles. Let’s look at the PDF method in a slightly different light. In the PDF method, we have a full probability distribution as the output at each timestep and we use that distribution to get the quantiles, which are the prediction intervals. Although for most parametric distributions, there are analytical formulae to get the quantiles, we can also find the quantiles numerically. We just draw N samples, where N is sufficiently large, and then calculate the quantiles of the drawn samples. The point is that even though we have a full probability distribution, for prediction intervals, all we need are the quantiles. And calculating quantiles given N samples is not dependent on the kind of distribution.

So, what if we can train our models to predict the specified quantiles directly, having no assumption on the underlying probability distribution? This is exactly what we do with quantile functions.

Before talking about quantile functions, let’s spend a minute on the Cumulative Distribution Function (CDF). This, again, is high school probability. In simple words, a CDF returns the probability of some random variable, X being less than or equal to some value, x:

$$F(x) = P(X \leq x)$$

where F is the CDF. This function takes in an input, x, and returns a value between 0 and 1. Let’s call this value $p$.

A quantile function is an inverse of CDF. This function tells you what value of x would make $F(x)$ return a particular value, $p$.

$$F^{-1}(p) = x$$

This function, $F^{-1}$, is the quantile function.

From the implementation perspective, for models which are capable of multi-output prediction (like the deep learning models), we can use a single model and predict all the quantiles we want to by changing the output layer. And for models that are restricted to a single output (like machine learning models), we can learn separate quantile models for each quantile.

Now, just like before, we can’t use point losses like the mean squared error. We got over this problem in PDFs by using the log likelihood function. And here, we can use quantile loss or pinball loss.

The quantile loss can be defined as below:

$$L_q(y_t, \hat{y}_t^q) = \begin{cases} (y_t - \hat{y}_t^q) \times q, & \text{if } y_t \geq \hat{y}_t^q \\ (\hat{y}_t^q - \hat{y}_t) \times (1 - q), & \text{if } y_t \lt \hat{y}_t^q \end{cases}$$

where $y_t$ is the target value at time t, $\hat{y}_t^q$ is the quantile forecast, and q is the quantile we are forecasting. The formula looks daunting, but bear with me for a second; it’s quite simple. Let’s try and get some intuitions about the loss. We know the median is the 0.5 quantile, and that is a measure of central tendency. But if we want our predictions to approximate the 75^th percentile or 0.75 quantile, then we would have to urge the model to overestimate, right? And if we want the model to overestimate, we need to penalize the model more if it’s underestimating. And vice versa, if we want to predict the 0.25 quantile, we need to underestimate. Quantile loss does exactly that:

• For $y_t >= \hat{y}_t^q$ (under estimation), the loss is $q \times (y_t - \hat{y}_t^q)$
• For $y_t \lt \hat{y}_t^q$ (over estimation), the loss is $(1 - q) \times (\hat{y}_t^q - y_t)$

The asymmetry is derived from the term q or 1 - q. The other term is just the difference between the actual and predicted values.

Let’s try to understand this with an example. Suppose we have the true value $y_t = 100$, and we want to estimate 0.75 quantile ($q = 0.75$).

• Case 1: Overestimation: $\hat{y}_t^q = 110$. Since $y_t \lt \hat{y}_t^q$, our quantile loss would be:

$$L_{q(100,110)} = (1 - 0.75) \times (110 - 100) = 2.5$$

• Case 2: Underestimation: $\hat{y}_t^q = 90$. Since $y_t \geq \hat{y}_t^q$, our quantile loss would be:

$$L_{q(100,110)} = 0.75 \times (110 - 100) = 7.5$$

Therefore, by varying the value of q, we can make the loss more or less asymmetric and toward either side. Figure 17.7 below shows the loss curves for q = 0.5 and q = 0.75 with these example predictions marked on it.

Figure 17.7: Quantile loss curves for q=0.5 and q=0.75

We can see that the quantile loss for q = 0.5 is symmetric as it represents the median, whereas the loss curve for q = 0.75 is asymmetric, penalizing underestimation a lot more than overestimation.

Although the formula has a branching structure, a maximum operation can be easily avoided when implementing it in code. The quantile loss in Python looks like this:

def quantile_loss(q, y, y_hat_q):
    """ Calculate the quantile loss for a given quantile.
    Args:
    q (float): The quantile to be evaluated, e.g., 0.5 for median.
    y (float): The target value.
    y_hat_q (float): The quantile forecast.
    """
    error = y - y_hat_q
    return np.maximum(q * error, (q - 1) * error)

Now, let’s get into the practical side of things and learn how to use quantile loss with the different techniques (both machine learning and deep learning) we have covered in this book.

Forecasting with quantile loss (machine learning)

Regular machine learning models are typically capable of modeling one output. Therefore, we will need to train a different model for each of the quantiles we are interested in using quantile loss. So, if we want to predict the 0.5, 0.05, and 0.95 quantiles for a problem, we will have to train three separate models, one for each quantile.

Let’s see how we can do that. Just like in the PDF section, we are using mlforecast to quickly whip up a synthetic problem and create some features. For the detailed code, refer to the full notebook, but for now, let’s assume that we have a dataframe data, which has all the features necessary to run a machine learning model. We have split it into train and val and then subsequently into X_train, y_train, X_val, and y_val.

The first step is to import the LGBMRegressor and set some parameters and quantiles we want to train.

params = {
    'objective': 'quantile',
    'metric': 'quantile',
    'max_depth': 4,
    'num_leaves': 15,
    'learning_rate': 0.1,
    'n_estimators': 100,
    'boosting_type': 'gbdt'
}
# converting levels to quantiles
# For 90% Confidence - 0.05 for lower, 0.5 for median, and 0.95 for upper
quantiles = [0.5] + sum([level_to_quantiles(l) for l in levels], [])

The key parameter to note here is that the objective is set as quantile and the metric is also set as quantile. The rest of the LightGBM parameters can be tuned and tweaked for each use case. Now, let’s train all the quantile models.

# Training a model each for the quantiles
quantile_models = {}
for q in quantiles:
    model = LGBMRegressor(alpha=q, **params)
    model = model.fit(X_train, Y_train)
    quantile_models[q] = model

Now that the model is trained, we can get the point prediction and prediction intervals from the quantile models.

# Point Forecast using the 0.5 quantile model
y_pred = quantile_models[0.5].predict(X_val)
# Prediction Intervals using the 0.1 and 0.9 quantile models
y_pred_lower = quantile_models[0.1].predict(X_val)
y_pred_upper = quantile_models[0.9].predict(X_val)

Now, let’s see what the forecast and metrics look like.

Figure 17.8: Forecast with prediction intervals using Quantile Regression with LightGBM

Figure 17.9: Metrics for Quantile Regression with LightGBM

The disadvantage of this is that we are training a model for each quantile. This can become unwieldy quickly. Instead of training one model, training three models would make the total training time increase. Another issue is that since all three models are trained differently, they might have different properties, and the way they have learned to solve the problem can also be very different. And because of this inconsistency, the prediction intervals can also suffer from some issues. We can see it clearly in the jagged prediction intervals in many of the time series in Figure 17.7; they seem disconnected from the median prediction. This is a problem that we don’t have in the deep learning world.

Forecasting with quantile loss (deep learning)

In the deep learning models, we use a common learning structure and just use different linear projections on the shared projection for the different quantiles. This ensures that the underlying representation and learning are common across all the quantiles and can cause a more coherent set of quantile predictions. So, for all the deep learning models we have learned about in the book, we can make them into a Quantile Forecast model by doing two things:

1. Instead of predicting point forecast (single number), we predict parameters of a probability distribution (one or more numbers).
2. Instead of using a point loss like Mean Squared Error, use a probabilistic scoring function like log likelihood.

And just like we did in the PDF section, all we need to do is to switch out the loss function in neuralforecast.

Let’s see how we can do this in neuralforecast (the library we were using Chapter 16). Just like before, we will take a simple model like an LSTM and the M4 competition dataset, but we can do the same with any model or any dataset because all we are doing is switching the loss function to MQLoss (multi-quantile loss).

Let’s start from the point where we have the data formatted the way neuralforecast expects in a Y_train_df and Y_test_df. The first thing we need to do is import the necessary classes.

from neuralforecast import NeuralForecast
from neuralforecast.models import LSTM
from neuralforecast.losses.pytorch import MQLoss

The only class that we haven’t looked at before is MQLoss. This class just calculates the quantile loss we discussed just now for multiple quantiles (which is how you would want to typically train the model). These are the major parameters of MQLoss:

• level: This is a list of floats that defines the different confidence levels we are interested in modeling. For instance, if we want to model 80% and 90% confidence, we should give the values as [80,90].
• quantiles: This is an alternate way of defining the level. Instead of 95% confidence, you can define them in quantiles → [0.1, 0.9].

Now, let’s set a horizon, the levels we need, and a few hyperparameters for LSTM.

horizon = 48
levels = [80, 90]
lstm_config = dict(input_size=3*horizon)

Now, we need to define the models we are going to use and the NeuralForecast class. Let’s define just one model.

models = [LSTM(h=horizon, loss=MQLoss(level=levels), **lstm_config)]
# Setting freq=1 because the ds column is not date, but instead a sequentially increasing number
nf = NeuralForecast(models=models, freq=1)

Notice how the syntax is exactly the same as the point forecast, except for the multi-quantile loss we chose. Now, all that’s left is to train the models.

nf.fit(df=Y_train_df)

Once the model is trained, we can predict using the predict method. This output will have the point forecast under the alias we have defined and the high and low intervals for all the levels we have defined.

Y_hat_df = nf.predict()

Now, we have a way of getting prediction intervals without making an assumption about the output distribution and that is valuable in real-world cases where we are not sure what the underlying output distribution would be. Let’s look at the generated probabilistic forecasts and its metrics.

Figure 17.10: Forecast with prediction intervals using Quantile Regression (deep learning)

Figure 17.11: Metrics for the Quantile Regression (deep learning)

We can see that the prediction intervals are quite in sync with each other and the median prediction and not disconnected like separate models with LightGBM. This is because the same learning is happening for all the quantiles, just the final projection heads are different. The coverage, in general, is also better.

There is another way to get prediction intervals that are dead simple, but not that easy to implement because of the way PyTorch is typically used. That’s what we are going to see next.

Monte Carlo Dropout

Dropout is a hugely popular regularization layer in deep learning. Without going into details, dropout regularization is when we randomly make some part of the weights of the network while training (dropout is turned off during inference). Intuitively, this forces the model to not rely on a few weights but rather to distribute the relevance of the weights across the network. From another perspective, we are applying a sort of regularization (very similar to Ridge regularization), which makes sure none of the weights are single-handedly too high to influence the output drastically.

Technically, besides making a random part of the weights zero, we also debias each layer by normalizing by the fraction of nodes/weights that were retained (not zeroed out). Let’s formalize this layer now. If the probability of dropout is $p$, and is applied on an intermediate activation, h then the activation after dropout, $h'$ will be:

$$h' = \begin{cases} 0 & \text{with } p \text{ probability} \\ \frac{h}{1 - p} & \text{otherwise} \end{cases}$$

Now, we know what dropout is. But remember that we were using dropout only during training as a regularization. In 2015, Yarin Gal et al. proposed that the good old dropout also doubles as a Bayesian Approximation of Gaussian Processes. That’s a lot of terms that we haven’t come across in a single phrase. Let’s take a short detour to understand these terms at a high level, and I’ll include other links in Further reading.

Bayesian inference is a statistical method that updates the probability of a hypothesis as more evidence or information becomes available. It is based on Bayes’ theorem, which mathematically expresses the relationship between the prior probability, the likelihood, and the posterior probability. Formally, Bayes’ theorem is given by:

$$p(\theta | D) = \frac{p(D | \theta) \cdot p(\theta)}{p(D)}$$

where $p(\theta | D)$ is the posterior probability of hypothesis $\theta$ given the data D, $p(D | \theta)$ is the likelihood (probability of observing the data D given hypothesis $\theta$), $p(\theta)$ is the prior probability of the hypothesis $\theta$ before observing the data, and $p(D)$ is the marginal likelihood or evidence (the total probability of observing the data under all possible hypotheses).

Although it has some specific terminology, this is very intuitive and provides a structured way to update our prior beliefs in the face of evidence or data. We start with a prior distribution $p(\theta)$ that represents our initial beliefs about the parameters. As we observe data D, we update our beliefs getting the posterior distribution $p(\theta | D)$. This posterior distribution combines the prior information and the likelihood of the observed data, providing a new, updated belief about the parameters. Further reading has a more detailed explanation for those who are interested. It also has a page from Seeing Theory that helps you visualize these in an intuitive way.

Now, let’s move on to the Gaussian Process (GP). As we have seen in Chapter 5, supervised learning is learning a function $\hat{y} = h(X, \phi)$ where $\hat{y}$ is the quantity we are interested in predicting, h, is the function we learn, X is the input data, and $\phi$ represents the model parameters. So, GPs assume this function as a probability distribution and use Bayesian inference to update the posterior of the function using the data we have available for training. It’s a drastically different way of learning from data and is inherently probabilistic.

There is one more term left, which is approximation. In many cases, the complex posterior distributions in Bayesian models make them intractable. So, we have a technique called Variational Inference in which we use a known parametric distribution family, q, and find a member that is closest to the true posterior. Getting into details of Variational Inference and GPs is out of the scope of the book, but I have added a few links in Further reading for those of you who are interested.

So, coming back to dropouts, Yarin Gal et al. showed that a neural network defined with dropout layers before each of the weight layers is in fact a Bayesian approximation of a GP. So, if the model with dropout is a GP and GP is a posterior over functions, this should give us a probabilistic output, right? But here, we don’t have a well-defined parametric probability distribution like the Normal distribution for us to analytically calculate the properties (like the mean, standard deviation, or quantiles) of the distribution. How do we do that?

Remember the discussion at the start of the Quantile Function section where we said that if we have N samples drawn from a distribution and if that N is sufficiently large, we can approximate the properties of the distribution? We have a name for that, and it’s called Monte Carlo sampling. Monte Carlo sampling is a computational technique used to estimate the statistical properties of a distribution by generating many random samples from that distribution. Bringing this idea to the dropout-enabled neural network, we can assess the properties of the posterior probability distribution of the function using Monte Carlo sampling, which means we need to keep dropout turned on during inference and sample from the posterior by executing the forward pass N times.

Theoretical justification allows us to apply dropout to any neural network, getting uncertainty estimates through a simple operation. Isn’t the simplicity of that beautiful?

So, all of that boils down to these simple steps to get prediction intervals for our forecasts:

1. Pick any deep learning architecture.
2. Insert Dropout Layers before every major operation and set them to a value $p$, where $p \gt 0$.
3. After training the model, with dropout enabled, do N forward passes.
4. Using the N samples, estimate the median (for the point prediction), and the quantiles corresponding to the defined confidence levels for prediction intervals.

This sounds simple enough, but it’s complicated for just one reason. Both PyTorch and Tensorflow are designed in a way that dropouts are turned off during inference. In PyTorch, we can indicate if the model is in the training phase or inference phase by doing model.train() or model.eval(), respectively. And most popular implementations that wrap PyTorch to make training easy and automated (like PyTorch Lightning) do this model.eval() step in the backend before predicting. So, when using libraries like neuralforecast (which uses PyTorch Lightning in the background), turning on dropout during prediction isn’t easy.

Before we learn how to implement MC Dropout for neuralforecast models, let’s take a slight detour and learn how to define custom models in neuralforecast. It is useful when you want to tweak any model for your use case. And we are doing it here because of two reasons:

1. I wanted to show you how to define a new model in neuralforecast.
2. I wanted a model that is ideal for the MC Dropout technique, i.e., the model needs to have dropouts before every weight/layer.

Creating a custom model in neuralforecast

Now it’s time to have some fun and define a custom PyTorch model that works with neuralforecast. To keep things simple, let’s use a model that is a small tweak on D-Linear we learned about in Chapter 16. Besides the linear trend and seasonality, we also add a component for non-linear trend. Let’s give it a wacky name as well—D-NonLinear. The architecture would be something like this:

Figure 17.12: D-NonLinear model architecture

Now, let’s understand how to write a model that would work with neuralforecast. All models in neuralforecast are inherited from one of three classes—BaseWindows, BaseRecurrent, or BaseMultivariate. The documentation clearly explains the purpose of BaseWindows, which is exactly what we require for our use case. We need to sample windows from a time series while training.

One more thing we need to keep in mind is that neuralforecast uses PyTorch Lightning under the hood for training. This link has more details on how to define a new model for neuralforecast: https://nixtlaverse.nixtla.io/neuralforecast/docs/tutorials/adding_models.html.

The full code for the model can be found in src/dl/nf_models.py, but we will look at key portions of the model definition right here.

Let’s start by defining the __init__ function (only including relevant portions here; refer to the Python file for the full class definition).

class DNonLinear(BaseWindows):
    def __init__(
        self,
        # Inherited hyperparameters with no defaults
        h,
        input_size,
        # Model specific hyperparameters
        # Window over which the moving average operates for trend extraction
        moving_avg_window=3,
        dropout=0.1,
        # Inhereted hyperparameters with defaults
        ...
        **trainer_kwargs,
    ):
        super(DropoutDNonLinear, self).__init__(
            h=h,
            ...
            **trainer_kwargs,
        )
        # Model specific hyperparameters
        self.moving_avg_window = moving_avg_window
        self.dropout = dropout
        # Model initialization to follow

We have now defined part of the __init__ function. Now, let’s initialize the different layers needed in the rest of the method. We have a series decomposition layer that uses a moving average to split the input into a trend and seasonal component, a linear trend predictor and seasonality predictor that takes the linear trend and seasonality and projects it into the future, and a non-linear predictor that takes in the original input and projects into the future.

        # Defining a decomposition Layer
        self.decomp = SeriesDecomp(self.moving_avg_window)
        # Defining a non-linear trend predictor with dropout
        self.non_linear_block = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, 100),
            nn.ReLU(),
            nn.Dropout(self.dropout),
            nn.Linear(100, 100),
            nn.ReLU(),
            nn.Dropout(self.dropout),
            nn.Linear(100, self.h),
        )
        # Defining a linear trend predictor with dropout
        self.linear_trend = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, self.h),
        )
        # Defining a seasonality predictor with dropout
        self.seasonality = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, self.h),
        )

Now, let’s define the forward method. The forward method should have just one argument that is a dictionary of different inputs:

• insample_y: The context window of the target time series we have to predict
• futr_exog: The exogenous variables for the future
• hist_exog: The exogenous variables for the context window
• stat_exog: The static variables

For this use case, we only need the insample_y since our model doesn’t use any other information. So, this is the forward method implementation:

    def forward(self, windows_batch):
        # Parse windows_batch
        insample_y = windows_batch[
            "insample_y"
        ].clone()  # --> (batch_size, input_size)
        seasonal_init, trend_init = self.decomp(
            insample_y
        )  # --> (batch_size, input_size)
        # Non-linear block
        non_linear_part = self.non_linear_block(
            insample_y
        )  # --> (batch_size, horizon)
        # Linear trend block
        trend_part = self.linear_trend(trend_init)  # --> (batch_size, horizon)
        # Seasonality block
        seasonal_part = self.seasonality(
            seasonal_init
        )  # --> (batch_size, horizon)
        # Combine the components
        forecast = (
            trend_part + seasonal_part + non_linear_part
        )  # --> (batch_size, horizon)
        # Map the forecast to the domain of the target
        forecast = self.loss.domain_map(forecast)
        return forecast

The code is pretty straightforward. The only thing we need to ensure to align to neuralforecast models is to take the data we need from the input dictionary and call self.loss.domain_map at the end so that it is mapped to the right output size depending on the loss. Now, this model will function just like any other model in the neuralforecast library.

Now, let’s get back to MC Dropout and its implementation.

Forecasting with MC Dropout (neuralforecast)

We said it wasn’t easy to implement MC Dropout in frameworks like neuralforecast and PyTorch Lightning, but just because something isn’t easy shouldn’t stop us from doing it. All we need to do is to make sure the dropouts are enabled during prediction and take multiple samples. If you are writing your own PyTorch training code, then it’s as simple as not calling model.eval() before predicting. But the best practice is to just make the dropouts into train mode and not the whole model. There may be layers like batch normalization, which also behave differently during inference, which might be affected. Let’s see a handy method that makes all dropout layers into train mode.

def enable_dropout(model):
    """Function to enable the dropout layers during test-time"""
    for m in model.modules():
        if m.__class__.__name__.startswith("Dropout"):
            m.train()

For neuralforecast, we have put together a recipe with which we can use MC Dropout for any of their models (provided they have enough dropouts). Now, we are going to use the custom DNonLinear model we just defined. Note that each component in the definition starts with a dropout layer so that we can apply MC Dropout with no qualms.

If you remember, we used PyTorch Lightning back in Chapter 13 and explained that it’s pretty much standard PyTorch code, but organized in a specified form—training_step, validation_step, predict_step, configure_optimizers, etc. For a refresher, head back to Chapter 13 and Further reading in the chapter to learn more about how to migrate from PyTorch to PyTorch Lightning. Since neuralforecast is already using PyTorch Lightning in the backend, the BaseWindows that we are inheriting is already a PyTorch Lightning model. This information is essential because we need to essentially modify the predict_step method to implement our MC Dropout.

Using the same inheritance we used to inherit BaseWindows, we can inherit the DNonLinear class we defined earlier and make a few changes so that it becomes an MC Dropout model. And for that, all we need to re-define is the predict_step method. The predict_step method is the method PyTorch Lightning calls every time it has to get a prediction for a batch. So, instead of taking the predictions as is, we need to keep the dropout enabled, take N samples from N forward passes, calculate the prediction intervals and median (point forecast), and return it.

class MCDropoutDNonLinear(DNonLinear):
    def predict_step(self, batch, batch_idx):
        enable_dropout(self)
        pred_samples = []
        # num_samples and levels will be saved to the model in MCNeuralForecast predict method
        for i in range(self.num_samples):
            y_hat = super().predict_step(batch, batch_idx)
            pred_samples.append(y_hat)
        # Stack the samples
        pred_samples = torch.stack(pred_samples, dim=0)
        # Calculate the median and the quantiles
        y_hat = [pred_samples.quantile(0.5, dim=0)]
        if self.levels is not None:
            for l in self.levels:
                lo, hi = level_to_quantiles(l)
                y_hat_lo = pred_samples.quantile(lo, dim=0)
                y_hat_hi = pred_samples.quantile(hi, dim=0)
                y_hat.extend([y_hat_lo, y_hat_hi])
        # Stack the results
        y_hat = torch.stack(y_hat, dim=-1)
        return y_hat

Are we done yet? Not quite. There is just one little thing to be done. In Chapter 16, we used a class called NeuralForecast for all the fitting and predicting of neuralforecast models. This is like a wrapper class that does the heavy lifting of preparing the inputs and outputs in the right way before calling the underlying models. This class has to be aware that we have tweaked the predict_step and therefore, we need to make a small change there. The solution is more of a hack than a principled way of editing, but a hack is just as good if it achieves the purpose. I have done the snooping around the implementation to figure out the best way to hack NeuralForecast to enable our MCDropout inference. There is no short way of explaining the hack, but just understand that I have misused the way neuralforecast flexibly produces point forecasts and prediction intervals based on different losses. So, here is the re-defined NeuralForecast class with a hack in the predict method.

class MCNeuralForecast(NeuralForecast):
    def __init__(self, num_samples, levels=None, **kwargs):
        super().__init__(**kwargs)
        self.num_samples = num_samples
        self.levels = levels
    def predict(
        self,
        df=None,
        static_df=None,
        futr_df=None,
        sort_df=True,
        verbose=False,
        engine=None,
        **data_kwargs,
    ):
        # Adding model columns to loss output names
        # Necessary hack to get the quantiles and format it correctly
        for model in self.models:
            model.loss.output_names = ["-median"]
            for l in list(self.levels):
                model.loss.output_names.append(f"-lo-{l}")
                model.loss.output_names.append(f"-hi-{l}")
            # Setting the number of samples and levels in the model
            model.num_samples = self.num_samples
            model.levels = self.levels
        return super().predict(
            df, static_df, futr_df, sort_df, verbose, engine, **data_kwargs
        )

That’s it. We have successfully “hacked” the library to do our bidding. In addition to this being an MC Dropout tutorial, it’s also a tutorial on how to hack a library to do what you want it to do. It is important to note that this doesn’t make you a “hacker,” so stop before you update your LinkedIn title.

Now, on to training the model. This is pretty much the same as you train other models in neuralforecast, but instead of the NeuralForecast class, you need to use the new MCNeuralForecast class we defined.

horizon = len(Y_test_df.ds.unique()) # 48
levels = [80, 90]
model = MCDropoutDNonLinear (
    h= horizon,
    input_size=WINDOW,
    moving_avg_window=horizon*3,
    dropout=0.1,
    max_steps=500,
    early_stop_patience_steps=5,
)
mcnf = MCNeuralForecast(models=[model], freq=1, num_samples=100, levels=levels)
mcnf.fit(Y_train_df, val_size= horizon, verbose=True)

Once the training is finished, we can generate the predictions like this:

Y_hat_df = mcnf.predict()
Y_hat_df = Y_hat_df.reset_index()

The output is exactly like any other model from neuralforecast with the prediction intervals formatted as <ModelName>-lo-<level> and <ModelName>-hi-<level>. The point forecast can be found under <ModelName>-median. In this case, <ModelName> would be MCDropoutDNonLinear.

Let’s look at the plot of the forecasts and the metrics.

Figure 17.13: Forecast with prediction intervals using MC Dropout

Figure 17.14: Metrics for MC Dropout

Our forecasting model does decent enough on the data; it’s nothing to write home about, but decent. If we do an ablation study, we might even realize that the non-linear component we added does absolutely nothing. But, as long as we had fun doing it and learned something from it, I’m happy. Now, look at the prediction intervals. They aren’t really smooth and have quite a bit of “noise” when you look at them, right? This is because of the inherent randomness in the methodology and may be because of insufficient learning. When we do MC Dropout, we are essentially relying on N sub-models or sub-networks and calculating the quantiles based on these N forecasts. Maybe a few of these sub-networks haven’t learned very well, and those outputs can skew the quantiles and thereby the prediction intervals.

There are many criticisms of the MC Dropout method. Many in the Bayesian community don’t consider MC Dropout as Bayesian and consider the Variational Approximation that was proposed such a poor approximation that we can’t refer to what it measures as Bayesian uncertainty. There is an unpublished Arxiv paper by Loic Le Folgoc et al. called “Is MC Dropout Bayesian?” (Reference 6), which claims that it isn’t. But it still doesn’t take away the fact that MC Dropout is a cheap way of getting uncertainty quantified. But when used in fields like medical studies, where uncertainty quantification is of paramount importance, we may want to take on something more principled.

We can also notice that the coverage is quite bad across all time series. And this is, again, something that is exhibited across different studies. In 2023, Nicolas Dewolf et al. published a study comparing different ways of uncertainty quantification for regression problems (Reference 7). They found that MC Dropout has one of the worst performances in both coverage and average length, underlying the claim that MC Dropout is a very dirty approximation of the uncertainty.

Now, let’s look at another technique for probabilistic forecasting that promises theoretical guarantees for perfect coverage and has become quite the rage in the last few years.

Conformal Prediction

What if I tell you that there is a technique for generating prediction intervals that statistically guarantees perfect coverage, can work on any model, and doesn’t require us to make any assumptions about the output distribution? Conformal Prediction is just that. Conformal Prediction is a method that helps machine learning models make reliable predictions by estimating how uncertain the model is. Conformal Prediction provides robust, statistically valid measures of uncertainty for any machine learning model, ensuring reliable and trustworthy predictions in critical applications.

Although it was proposed as early as 2005 by Vladmir Vovk (Reference 8), it picked up interest in the last couple of years. Let’s first understand the basic principles of Conformal Prediction using a classification example and then see how we can do it for regression and time series examples.

Conformal Prediction for classification

Let’s start with a trained model, $\hat{f}$, which outputs estimated probabilities (softmax scores) for K output classes ($\hat{f}(x) \in [0,1]^K$). It doesn’t matter what this model is; it can be a machine learning model, a deep learning model, or even a rule-based model. We have training data, $(X_{train}, Y_{train})$, and test data, $(X_{test}, Y_{test})$. Now, we need a small amount of additional data (other than training and test) called calibration data, $(X_{calib}, Y_{calib}) = \{(x_1, y_1), ..., (x_n, y_n)\}$. Now, what do we want from this? Using $\hat{f}$ and $X_{calib}$, we want to create a prediction set of possible labels, $\mathcal{C}(X_{test})$ that makes sure that the probability that a test data point is part of that set is almost exactly the user-defined error rate, $\alpha$ (we will be talking about error rates throughout this discussion. A 10% error rate is a 90% confidence level). This is exactly what Conformal Prediction guarantees. It’s called the marginal coverage guarantee and it can be written more formally as:

$$1 - \alpha \leq \mathbb{P}(Y_{test} \in \mathcal{C}(X_{test})) \leq 1 - \alpha + \frac{1}{n+1}$$

Naturally, you may have this question in your mind. What is this $\frac{1}{n+1}$? This term signifies that the coverage guarantee is derived from a finite sample of size n. Since n is in the denominator, we know as and when n increases, this term becomes smaller and smaller. Extending this to the limit, we know that if $n = \infty$, this term would be zero and the coverage would be exactly $1 - \alpha$.

We know what we want, but how do we get it? The core idea in conformal prediction is very simple and can be laid out in four steps:

1. Identify a heuristic notion of uncertainty using the trained model. In our classification example, this can be the softmax scores.
2. Define a score function, $s(\hat{y}, y) \in \mathbb{R}$, which is also called Non-Conformity Scores. This can be a score that takes in the prediction, $\hat{y}$, and actual value, $y$, and gives a score that encodes the disagreement between them. The higher the score is, the larger the disagreement is. In the classification example, this would be something as simple as $1 - \hat{f}(x)_{Y_i}$. In simple English, this means taking the softmax score of the correct class and doing $1 - score$.
3. Compute $\hat{q}$ as the $\frac{(1 - \alpha)(n + 1)}{n}$ quantile of the calibration scores. We use the calibration data and score function to calculate calibration scores and calculate the quantile on that data. The $\frac{n + 1}{n}$ is the quantile calculation is again derived from the finite sample correction. As $n$ tends to infinity, the term tends to zero.
4. Use this quantile to form prediction sets for new examples: $\mathcal{C}(x) = \{y : s(x, y) \leq \hat{q}\}$. This means selecting all the items from the output set that has a score (according to the score function) greater than the threshold, $\hat{q}$.

This simple technique will give us prediction sets that guarantee to satisfy the marginal coverage, no matter what model is used or what the distribution of the data is. Let’s see how simple this is using Python code and assuming that the model we are talking about is a scikit-learn classifier.

We have a trained model, model, calibration data, X_calib, and test data, X_test. For full code and some visualizations, check the notebook.

# 1: Get conformal scores
n = calib_y.shape[0]
cal_smx = model.predict_proba(calib_x) # shape (n, n_classes)
# scores from the softmax score for the correct class
cal_scores = 1 - cal_smx[np.arange(n), calib_y] # shape (n,)
# 2: Get adjusted quantile
alpha = 0.1  # Confidence level (1 - alpha)
q_level = np.ceil((n + 1) * (1 - alpha)) / n
qhat = np.quantile(cal_scores, q_level, method='higher')
# 3: Form prediction sets
val_smx = model.predict_proba(test_x)
prediction_sets = val_smx >= (1 - qhat)

Now, let’s think about the prediction set, $\mathcal{C}$. We have been defining it as set-valued with discrete classes for the classification scenario. This set becomes larger or smaller based on how confident the initial heuristic estimate of uncertainty is.

Conformal Prediction for regression

Let’s extend this notion of prediction sets to regression. In regression, the output space is continuous rather than discrete, and we aim to construct continuous prediction sets, which are typically a continuous interval in $\mathbb{R}$. The idea is to maintain the same principle of coverage: the prediction interval should contain the true value with high probability. So, now the prediction set, $\mathcal{C}$, that we saw earlier is also the prediction interval in the regression context. But along with the change in the interpretation of prediction sets, we will also need to change the score function, which calculates non-conformity scores. A common score function that is used is the distance to the conditional mean, $s(x, y) = |y - \mu(x)|$ (Reference 10). When we have a trained model, $\hat{f}$, we can consider the output of the model as the conditional mean, which will make this the absolute residual value for each point.

$$s(x, y) = |y - \hat{f}(x)|$$

Note that this score satisfies the condition. The larger the deviation, the larger the “heuristic” measure of uncertainty is. The rest of the procedure remains almost the same—calculating the quantile, $\hat{q}$, and forming the prediction intervals, $\mathcal{C}(x) = [\hat{f}(x) - \hat{q}, \hat{f}(x) + \hat{q}]$

Let’s check how the Python code changes (full code is in the notebook).

# 1: Get conformal scores
calib_preds = model.predict(calib_x)
cal_scores = np.abs(calib_y - calib_preds)
# 2: Get adjusted quantile
qhat = … # Exactly the same as classification
# 3: Form prediction intervals
test_preds = model.predict(test_x)
lower_bounds = test_preds - qhat
upper_bounds = test_preds + qhat

It’s as simple as that. We can check coverage and see that it will be greater than 90%, which is the error rate we defined with $\alpha = 0.1$.

But you might have noticed something. In this method, we have the same width of the interval all throughout. But we would expect the intervals to be tighter when the model is more confident and wider when it isn’t (let’s call it adaptive prediction intervals from now on). Let’s look at another technique that has this property.

Conformalized Quantile Regression

We learned about Quantile Regression as one of the methods for probabilistic forecast (or regression in a general case). Quantile Regression is powerful in the sense that it does not require us to have any prior assumption about the underlying output distribution. But it does not enjoy the coverage guarantees that Conformal Prediction offers. In 2019, Yaniv Roano et al. (Reference 11) brought the best of both worlds into Conformalized Quantile Regression (CQR). They proposed a way to take the quantile forecasts and conformalize them such that they have the coverage guarantees that conformal prediction assures.

In this case, the model to be used has a restriction. It should be a model that outputs quantile predictions. And, from our earlier discussion, we know that if the error rate is $\alpha$, then the quantiles we need for the prediction intervals are $\hat{y}_t^{\alpha/2}$ and $\hat{y}_t^{1 - (\alpha/2)}$.

So, the quantile model predicts $\hat{y}_t^{\alpha/2}$ and $\hat{y}_t^{1 - (\alpha/2)}$. By definition, if $\hat{y}_t^{\alpha/2}$ and $\hat{y}_t^{1 - (\alpha/2)}$ are true estimations of the real quantiles, a quantile regression alone will have perfect coverage. But the model fit may not be perfect and that would result in sub-par coverage. We will use conformal prediction to correct the quantiles based on calibration data such that we get the perfect coverage promised by conformal prediction.

Yaniv Roano et al. proposed to use a new non-conformity score function for Quantile Regression.

$$s(x, y) = \max \{\hat{y}_t^{\alpha/2} - y, y - \hat{y}_t^{1 - (\alpha/2)}\}$$

Let’s take a beat and explore the score function using the diagram in Figure 17.15.

Figure 17.15: Illustration of score function for Conformalized Quantile Regression

There are two terms in the max operator. If the true value, $y$, is between the two quantiles, $\hat{y}_t^{\alpha/2}$ and $\hat{y}_t^{1 - (\alpha/2)}$, both the terms inside will be negative and will be the distance to the nearest prediction interval (see points B and C in Figure 17.15).

Now, let’s look at point A (3), which is above the higher quantile. The quantiles are [1, 2]. The two terms in the max operator would be:

$$\hat{y}_t^{\alpha/2} - y = 1 - 3 = -2$$

$$y - \hat{y}_t^{1-(\alpha/2)} = 3 - 2 = 1$$

This makes the score function 3 because of the max operator. And now, we see point D (1.6), which is below the lower quantile (quantiles are [3.5, 2]). This makes the score function max{0.4, -1.9}, which would be 0.4. So, the max operator ensures that the score is positive if it falls outside the quantiles.

Therefore, what we have is a score function that assigns positive values to points where the actual value falls outside the intervals and negative to points within. And for the points that fall outside the intervals, the way the score function is constructed will choose the worse error. This satisfies our requirements from the score function. A larger score shows larger uncertainty, and it also encodes a heuristic notion of uncertainty.

Now that we have the scores, the rest of the steps are almost identical:

1. Compute $\hat{q}$ as the $\frac{(1 - \alpha)(n + 1)}{n}$ quantile of the calibration scores.
2. Use this quantile to form prediction sets for new examples: $\mathcal{C}(X_{test}) = [\hat{y}_t^{\alpha/2} - \hat{q}, \hat{y}_t^{1-(\alpha/2)} + \hat{q}]$, i.e., we widen the exiting quantiles by $\hat{q}$ and get coverage guarantees.

Quantile regression is one of the better ways of getting adaptive prediction intervals. Let’s look at one such technique now.

Conformalizing uncertainty estimates

If we think about Conformalized Quantile Regression (from the last section) a bit deeply, we can realize that what the underlying quantile regression is doing is capturing the uncertainty estimate at each point in our prediction. And we conformalize those estimates for better coverage. If we can capture this uncertainty estimate, $u(x)$, we have some hope of conformalizing this to get better adaptive prediction intervals.

Let’s say that we have a trained model, $\hat{f}$, and some uncertainty scalar, $u(x)$, that has high values when uncertainty is high, and vice versa. We can define our non-conformity score as:

$$s(x, y) = \frac{|y - \hat{f}(x)|}{u(x)}$$

The natural interpretation of this score is that we are multiplying a correction factor to the standard $|y - \hat{f}(x)|$. Once we have this new score, the rest of the process is exactly the same—taking $\hat{q}$ from the scores and forming the prediction intervals as $\mathcal{C}(x) = [\hat{f}(x) - u(x) \cdot \hat{q}, \hat{f}(x) + u(x) \cdot \hat{q}]$.

So, what are some ways to capture this uncertainty? (Note that this uncertainty measure should be capturing it at a data point level.)

• Assume a probability distribution and modeling their parameters (Probability Density Function). In the case of Gaussian Distribution, we would have an estimate of uncertainty as the standard deviation, $\hat{\sigma}(x)$, and we consider that as $u(x)$. $u(x) = \hat{\sigma}(x)$.
• Use the MC Dropout technique to generate samples and calculate the standard deviation of $\hat{f}(x)$ from the samples: $u(x) = SD_{samples}[\hat{f}(x)]$.
• Along with the main model prediction $y|X$, train another model to predict the residuals, $\hat{r} = (\hat{y} - y)|X$, and set $u(x) = \hat{r}(x)$.
• Use an ensemble of models to generate multiple predictions for each data point and take the standard deviation of different predictions at each data point: $u(x) = SD_{ensemble}[\hat{f}(x)]$.

Although the list above is not exhaustive, it does show that we can apply conformal prediction to almost any uncertainty estimates (including the ones we have already seen in the chapter). This makes the conformal prediction paradigm a very flexible toolkit to get coverage guarantees with a wide variety of problems. Even with this flexible nature, there are cases that mess with the coverage guarantees that the framework promises. And for our time series context, this is important to understand.

Exchangeability in Conformal Prediction and time series forecasting

Exchangeability is a fundamental assumption in Conformal Prediction. Exchangeability means that the data points are identically distributed and their joint probability distribution does not change when the order of the data points is changed. This concept ensures that past data points can reliably predict future data points.

Imagine a chocolate factory producing chocolates with consistent weights. If the production process is highly controlled, with the same conditions and ingredients, the weights of the chocolates are exchangeable because their order of production does not affect their weight. You can sample 100 chocolates, measure their weights, and calculate nonconformity scores based on deviations from the predicted weight. Using these scores, you can form prediction intervals for future chocolates. Because the chocolates are exchangeable, the sample distribution represents the future distribution, making prediction intervals reliable.

However, if the production process changes over time—due to machinery wear or different ingredient batches—the weights become non-exchangeable. The order of production affects the weights, making the sample distribution unrepresentative of future weights, leading to unreliable prediction intervals.

In time series data, observations are typically dependent on previous observations, violating the exchangeability assumption. For example, in a sales forecast, today’s sales may influence tomorrow’s sales due to trends, or the sales a year ago may influence tomorrow’s sales due to seasonal effects. This dependency means that the distribution of past data does not represent the distribution of future data accurately.

But what does this mean for us? The most obvious answer is that our coverage guarantees will suffer. But can we still apply these techniques for time series? Of course we can. Empirically, the community has seen that this framework works for time series data as well, but with some loss in coverage guarantees. For most practical purposes, there shouldn’t be an issue in using regular conformal prediction for time series data. In 2023, Barber et al. (Reference 12) studied this issue and derived theoretical coverage guarantees for non-exchangeable data (like time series). They defined the Coverage Gap as the difference between expected coverage ($1 - \alpha$) and actual coverage and derived an upper bound on this gap to show how much the exchangeability assumption on the scores is violated. For this bound, they considered the scores of the calibration data with our original model, $s(z)$, and an alternate model, which was trained on the same data but after swapping one randomly selected datapoint in the training data with the test datapoint, $s(z^i)$.

The bound was shown to be directly proportional to $d_{TV}(s(z), s(z^i))$, which is the distributional distance between these two scores. In most algorithms we use, swapping one data point may not change the model drastically and therefore, we can still use conformal prediction for time series data with minimal loss in coverage.

But on the other hand, if we want to be really accurate with the prediction intervals, or if we are using a model that is particularly impacted by the swapping of a datapoint, then we would need some techniques to overcome this degradation due to a shift in distributions. There are many ways to deal with this, and it is an active area of research at the time of writing the book. There are two very simple methods that are worth mentioning here.

Weighted conformal prediction

Suppose we have slowly varying changes in the data distribution in a time series, $(x_1, x_2, ..., x_N)$, and we are using a calibration set, $X_c = (x_{N-k}, ..., x_N)$, taking the last k timesteps. And we are interested in predicting the test set, $X_{test} = (x_{N+1}, ..., x_{N+H})$, where H is the horizon of forecast.

So, it stands to reason that the most recent timestep in $X_c$ would be closest to the distribution of values we would observe in the test time period. So, what if we assign weights to the non-conformity scores in the calibration data such that the most recent time step gets higher weights, and calculates a weighted quantile instead of a regular quantile? Apparently, that’s a very good idea, with some solid theoretical backing as well.

And weighing the calibration data using recency is just one of the ways we can use the weights to tackle the distribution shift. More generally, any weight schedule, $w_1, ...w_k, w_i \in [0, 1]$ can be used here. Maybe for a strongly seasonal time series, it makes sense to use seasonal periods to define the weights, or there may be some other known criteria that makes different instances of the calibration data more or less relevant to future prediction.

Before we get into the actual mechanics of this, we need to understand what weighed quantiles are. If you are already comfortable with the concept, feel free to go ahead. If you need some intuition about what it is, I strongly advise you to check out the notebook in the chapter folder named 08-Quantiles_and_Weighted_Quantiles.ipynb.

Now, let’s come back to our Weighted Conformal Prediction method.

So, as we discussed earlier, for any normalized weight schedule, $w_i \in [0,1] \text{ and } \sum_{i=1}^{k} w_i = 1$, and calibration scores, s_i the weighted quantile can be formally defined as:

$$\hat{q} = \inf\left\{q : \sum_{i=1}^{n} w_i \mathbb{I}\{s_i \lt q\} \geq 1 - \alpha\right\}$$

where inf is the infimum and $\mathbb{I}$ is an indicator function, which is 1 when the condition is true and 0 otherwise. In this context, the infimum is the smallest value of q such that the inequality holds true. This is just a more rigorous way of defining the weighted quantile that we saw in the aforementioned notebook. The rest of the process is exactly the same as before.

Practically, there are a few different ways we can use this for time series problems. For instance:

1. We can consider a moving window of length K and have a fixed weight vector of length K. Under this scheme, we will apply the weights for each point in the time series to the last K points and calculate the prediction interval for that point. These weights can be equal weights or even decayed weights, capturing the temporal element.
2. When we model multiple time series together, we can make sure the weights reflect how close another time series is to the time series we are generating the intervals for, along with the temporal context.

The bottom line is that we can be as creative as we can when coming up with the weights. The guideline is that the weight should reflect how different the calibration datapoints are from the datapoint you are generating the prediction interval for. Remember, we talked about the upper bound earlier, $d_{TV}(s(z), s(z^i))$. The weight we choose would counteract this term. When we have smaller weights for datapoints that are “farther” away from the datapoint we care about, this brings down the upper bound of the coverage gap and makes it tighter.

Now, let’s learn about another very simple modification to account for distribution shift.

Adaptive Conformal Inference (ACI)

In 2021, Gibbs et al. (Reference 13) proposed another way to deal with distribution shift (especially in time series) in an online setting. Time series data usually comes in one datapoint at a time, and this way to deal with distribution shift relies on this online aspect as it proposes to keep adjusting the prediction intervals based on the data that keeps trickling in, thus making the prediction intervals adapt to changing distributions. This method, called Adaptive Conformal Inference (ACI), can be integrated with any prediction algorithm to provide robust prediction sets under non-stationary conditions.

In traditional Conformal Prediction, we have a score function, $s(x, y)$, and a quantile function, $\hat{q}$, which gives us the prediction intervals, $\hat{C}$. Note that the underlying uncertainty model, which was conformalized as $\hat{q}$, can be any way of estimating uncertainties, like quantile regression, PDF, MC Dropouts, and so on. When the data is exchangeable, the $\hat{q}$ we calculated on calibration data will hold good on future test datapoints. But when the distribution shifts, this $\hat{q}$ will also start to become less and less relevant. To address this, the authors propose regularly re-estimating these functions to align with the most recent data observations. Specifically, at each time point, t, a new score function, $s_t(.)$, and a new quantile function, $\tilde{q}_t$(.), are fitted based on the most recent data.

For this, they defined the miscoverage rate, $M_t(\alpha)$, as the probability that the true label, $Y_t$, lies outside the prediction intervals, $\hat{C}$, where probability is calculated over the calibration data and the test datapoint. We want the Miscoverage rate, $M_t(\alpha)$, to be equal to $\alpha$ (expected error rate). But since the data distribution is shifting, $M_t(\alpha)$ is not expected to remain constant over time, and it may not equal the target level, $\alpha$. The authors hypothesize that for each time, $t$, there may exist an optimal coverage level, $\alpha_t^*$ such that the miscoverage rate, $M_t(\alpha_t^*)$, is approximately $\alpha$.

To estimate this $\alpha_t^*$, the authors propose a simple online update equation. This update takes into consideration the empirical miscoverage rate of the previous observations and then decreases or increases our estimate of $\alpha_t^*$. Concretely, if we set $\alpha_1 = \alpha$, we can define error as:

$$err_t = \begin{cases} 1, & \text{if } Y_t \notin \hat{C}(a_t), \\ 0, & \text{otherwise}, \end{cases}$$

Now, we can define the update step recursively as:

$$\alpha_{t+1} = \alpha_t + \gamma(\alpha - err_t)$$

Here, $err_t$, serves as an estimate of the historical miscoverage rate and $\gamma \gt 0$ is the step-size (a hyperparameter; more on this later). So, when $err_t = 0$ (prediction was within the interval), $\alpha - err_t$ will be positive and thus updating $\alpha_{t+1}$ to be higher than $\alpha_t$. This, in turn, makes the prediction interval narrower (according to the $\gamma$ we have defined). With the same logic, when $err_t = 1$ (prediction was outside the interval), the prediction interval becomes wider.

A natural alternative to this update, which also takes into account the history a bit more, is using a weighted average of past timesteps:

$$\alpha_{t+1} = \alpha_t + \gamma \left( \alpha - \sum_{i=1}^{t} w_i \cdot err_i \right)$$

where $\{w_i\}_{1 \leq i \leq t} \in [0,1]$ is a sequence of increasing weights with $\sum_{i=1}^{t} w_i = 1$. Instead of looking at just the last time step for the estimate of miscoverage, this update looks at the recent history. This makes it slightly more robust, at least in theory. The paper reported no significant difference between the two strategies. One of the strategies they have used to decide the weights is:

$$w_i = \frac{0.95^{t-s}}{\sum_{s=1}^{t} 0.95^{t-s}}$$

They reported that the trajectories of prediction intervals that they obtained from the simple update and weighted update were almost the same, but the weighted one was considerably smoother with less local variation in $\alpha_t$.

Now, let’s also spend some time understanding the effect of the step-size parameter, $\gamma$. The intuition is very similar to the learning rate in deep learning models. $\gamma$ decides the magnitude with which we update the $\alpha$. The larger the value, the quicker the update, and vice versa. The paper also gives us an intuition that the greater the distributional shift greater the value of $\gamma$. For all their experiments, they used $\gamma = 0.005$ with a justification that they found this value to make the trajectories relatively smooth while still being large enough to allow $\alpha_t$ to adapt to distributional shifts. We can see this as a parameter controlling the strength of “adapting” and letting us move in the spectrum between non-adaptive intervals and strongly adaptive intervals.

Now, let’s see how to apply these in practice.

Forecasting with Conformal Prediction

We didn’t find any ready-to-use implementations of all the techniques we wanted to show here, especially one that has these properties:

• Complete separation of model layer and conformal prediction layer (being model-agnostic is one of the most exciting features of Conformal Prediction)
• Out-of-the box compatibility with neuralforecast predictions
• Time series focus
• Pedagogical ease

Therefore, we have included a file (src/conformal/conformal_predictions.py) with the necessary implementations that would work with neuralforecast forecasts and have a unified API. It is also simple enough to understand. We will go through major parts of the code, but to see how it all fits together, you should just take a look at the file.

All the methods we discussed, like Conformal Prediction for Regression, Conformalized Quantile Regression, and Conformalizing uncertainty estimates, have been coded out in the same API. Let’s look at the most basic Conformal Prediction for Regression to understand the API. It can be found in the ConformalPrediction class in the file. The rest of the techniques inherit this class and make slight tweaks. And all these classes are coded in such a way that they take in the prediction dataframe from neuralforecast (or statsforecast) and use the same naming conventions to conformalize those predictions. In theory, any forecast that can be made into the expected format can be used with these classes. The expected columns in the format are:

• ds: This column should contain dates or the numerical equivalent of time.
• y: For train and calibration datasets, this column is necessary and it represents the actual value of that time series.
• unique_id: This column is the unique identifier for different time series.

In addition to these columns, we would also have a column (or multiple columns) of forecast, named accordingly.

Before we start generating conformal predictions, we also need some data and forecasts. We are using the same data that we have been using in this chapter (M4), with one additional split (calibration) created. And using the new train data, we have created these three forecasts with level = 90:

1. LSTM point forecast (LSTM)
2. LSTM with Quantile Regression (LSTM_QR)
3. LSTM with PDF (normal distribution) (LSTM_PDF)

The notebook has the entire code, but we can start from the point where we have already split the data into Y_train_df, Y_calib_df, and Y_test_df, and generated and stored forecasts in a dictionary, prediction_dict. Let’s take a look at the top five rows of the prepared data frame to see what kind of data we are working with.

Figure 17.16: Top five rows of the Y_calib_df we are working with. This is the format that the conformal prediction classes we have coded expect

Now, let’s get down to business and start creating prediction intervals.

Conformal prediction for regression

The ConformalPrediction class provides a structured way to calculate prediction intervals based on a chosen model’s predictions from a calibration dataset. It includes the following input parameters:

• model (str): The name of the column with the forecast you want to conformalize. This is a required parameter.
• level (float): The confidence level for the prediction intervals, expressed as a percentage (e.g., 95 for a 95% confidence interval). The level must be between 1 and 100. This is a required parameter.
• alias (str, optional): An optional string to provide an alias for the model. This can be useful when working with multiple models or versions and you want to call the output something else other than the model. If not provided, model is used.

The major functions for using the class are:

• fit(Y_calib_df): This method orchestrates the entire calibration process. It first calculates the calibration scores using calculate_scores and then determines the quantiles for each unique_id using get_quantile. The resulting quantiles (q_hat) are stored as an attribute of the class, making them available for subsequent prediction intervals.
• predict(Y_test_df): This method applies the prediction intervals to the test data. It uses the calc_prediction_interval method to compute the intervals and then adds them to the DataFrame as new columns. The new columns are created in this format: f"{self.alias or self.model}-{self._mthd}-lo-{self.level}". For example, the higher interval for Conformal Prediction using LSTM would have LSTM-CP-hi-90 as the column name in the dataframe.

These classes are the external API. Internally, there are some methods that actually define how it’s done. Let’s look at the major methods in the context of regular Conformal Prediction.

The calculate_scores method is defined below, where we just calculate the absolute residuals using the calibration dataset as the scores:

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
    Y_calib_df["calib_scores"] = np.abs(Y_calib_df["y"] - Y_calib_df[self.model])
    return Y_calib_df

The get_quantile method calculates the quantile using the defined $\alpha$ for each unique_id:

def get_quantile(self, Y_calib_df):
    def get_qhat(Y_calib_df):
        n_cal = len(Y_calib_df)
        q_level = np.ceil((n_cal + 1) * (1 - self.alpha)) / n_cal
        return np.quantile(
                Y_calib_df["calib_scores"].values, q_level, method="higher"
            )
    return Y_calib_df.groupby(
"unique_id").apply(get_qhat).to_dict()

The calc_prediction_interval method uses the calculated q_hat and mean prediction to generate the prediction intervals. For regular Conformal Prediction, it goes like this:

def calc_prediction_interval(self, Y_test_df, q_hat):
    return (
            Y_test_df[self.model] - Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.model] + Y_test_df["unique_id"].map(q_hat),
        )

Now, let’s use it for forecasting.

from src.conformal.conformal_predictions import ConformalPrediction
Y_calib_df, Y_test_df = prediction_dict['LSTM']
# Y_calib_df & Y_test_df have forecasts in column named "LSTM"
cp = ConformalPrediction(model="LSTM", level=level)
# Calibrating the model
cp.fit(Y_calib_df=Y_calib_df)
# Generating Prediction intervals
Y_test_df_cp = cp.predict(Y_test_df=Y_test_df)

The generated dataframe (Figure 17.15) with prediction intervals will have two columns—LSTM-CP-lo-90 and LSTM-CP-hi-90 for the lower and upper prediction interval, respectively. CP is the method tag we have assigned the conformal prediction class.

We can check the method name of any object as follows:

>> cp.method_name
'Vanilla Conformal Prediction (CP)

Let’s see what the forecast dataframe looks like (CP is the method tag we have assigned the conformal prediction class. We can check the method name of any object by doing cp.method_name):

Figure 17.17: The generated dataframe with prediction intervals

Now, we also have to calculate the coverage and average length of the intervals to assess how these prediction intervals are. We use the same methods we used earlier to do that. Instead of looking at the performance of each method, let’s save the discussion for the end and look at creating the prediction intervals for now.

So, let’s move on to the next technique.

Conformalized Quantile Regression

The first condition of applying this method of conformal prediction is that there should already be a set of prediction intervals from the underlying Quantile Regression. Therefore, we use the LSTM_QR model we trained here.

The main difference between vanilla Conformal Prediction and CQR is the way the scores are calculated and the prediction intervals. So, we can inherit ConformalPrediction and just redefine these two methods.

Let’s look at the calculate_scores method.

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
    lower_bounds = Y_calib_df[self.lower_quantile_model]
    upper_bounds = Y_calib_df[self.upper_quantile_model]
    Y_calib_df["calib_scores"] = np.maximum(
            lower_bounds - Y_calib_df["y"], Y_calib_df["y"] - upper_bounds
        )
    return Y_calib_df

We have just implemented the formula we saw earlier. self.lower_quantile_model and self.upper_quantile_model are the column names of the already-generated intervals from CQR.

Now, we also need to define the calc_prediction_interval method.

def calc_prediction_interval(self, Y_test_df, q_hat):
    return (
            Y_test_df[self.lower_quantile_model] - Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.upper_quantile_model] + Y_test_df["unique_id"].map(q_hat),
        )

q_hat is a dictionary with the quantiles calculated for each unique_id. So, all we do here is take the existing prediction interval from CQR and adjust it by mapping the quantile using the unique_id in the input dataframe.

Now, let’s use this for forecasting. The API is exactly the same as before.

from src.conformal.conformal_predictions import ConformalizedQuantileRegression
Y_calib_df, Y_test_df = prediction_dict['LSTM_QR']
# Forecast in column "LSTM_QR"
cp = ConformalizedQuantileRegression(model="LSTM_QR", level=level)
cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_cqr = cp.predict(Y_test_df=Y_test_df)

Now, let’s look at the third technique we discussed.

Conformalizing uncertainty estimates

If you remember the discussion we had earlier, to use this technique, we need an estimate of uncertainty that can be further conformalized. This is why we picked one of the other techniques we used earlier, PDF. But we can also do this with the MC Dropout just as easily. All we need is the standard deviation or something similar that captures the uncertainty at each data point.

We are using the LSTM_PDF forecast that we generated earlier for this purpose. Although the model predicts the mean and standard deviation of the normal distribution, it is used internally to generate the prediction intervals. So, the output from the PDF model we defined earlier would be the prediction intervals, but we want the standard deviation. Fear not. We know the prediction intervals were created with the normal distribution. So, it’s not hard to re-engineer the standard deviation from the prediction intervals.

$$Lower\ Bound = \mu - Z \cdot \sigma$$

$$Upper\ Bound = \mu + Z \cdot \sigma$$

Using basic math, we can derive:

$$\sigma = \frac{Upper \ Bound - Lower \ Bound}{Z}$$

And Z is very straightforward to get. We can use scipy.stats.norm for this. Below is a method that will get you the standard deviation from the prediction intervals (do keep in mind that this is only for PDFs created using a normal distribution).

from scipy.stats import norm
def calculate_standard_deviation(upper_bound, point_prediction, confidence_level):
    # Calculate the Z-value from the confidence level
    z_value = norm.ppf((1 + confidence_level) / 2)
    # Calculate the standard deviation
    sigma = (upper_bound - point_prediction) / z_value
   
    return sigma
def reverse_engineer_sd(X, model_tag, level):
    X["std"] = calculate_standard_deviation(
        X[f"{model_tag}-hi-{level}"], X[model_tag], level / 100
    )
    return X

Now, we add this to our Y_calib_df and Y_test_df.

Y_calib_df = reverse_engineer_sd(Y_calib_df, "LSTM_Normal", level)
Y_test_df = reverse_engineer_sd(Y_test_df, "LSTM_Normal", level)

Now, let’s look at how we can define the class. We need additional information in here that we didn’t need earlier—the column name of the uncertainty estimate. So, we define our new class (still inheriting ConformalPrediction) as below:

class ConformalizedUncertaintyEstimates(ConformalPrediction):
    def __init__(
        self,
        model: str,
        uncertainty_model: str,
        level: Optional[float] = None,
        alias: str = None,
    ):
        super().__init__(model, level, alias)
        self.method = "Conformalized Uncertainty Intervals"
        self._mthd = "CUE"
        self.uncertainty_model = uncertainty_model

We define an additional parameter, uncertainty_model, and pass on the other parameters to the parent class.

Now, it’s pretty straightforward. We need to define how the scores are calculated:

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
        uncertainty = Y_calib_df[self.uncertainty_model]
    Y_calib_df["calib_scores"] = (
            np.abs(Y_calib_df["y"] - Y_calib_df[self.model]) / uncertainty
        )
    return Y_calib_df

And the calc_prediction_interval method:

def calc_prediction_interval(self, Y_test_df, q_hat:
    uncertainty = Y_test_df[self.uncertainty_model]
    return (
            Y_test_df[self.model] - uncertainty * Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.model] + uncertainty * Y_test_df["unique_id"].map(q_hat),
        )

That’s it. Now, we have a new class that conformalizes uncertainty estimates. Let’s use it to get the forecast for the dataset we have been working with.

from src.conformal.conformal_predictions import ConformalizedUncertaintyEstimates
# We have saved uncertainty estimates in "std"
cp = ConformalizedUncertaintyEstimates(model="LSTM_Normal", uncertainty_model="std", level=level)
cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_pdf = cp.predict(Y_test_df=Y_test_df)

Now, let’s also look at the two techniques we saw, which were more suited for time series problems where there is a distribution shift.

Weighted conformal prediction

We saw earlier that weighted conformal prediction was just about applying the right kind of weights to the calibration data such that the points similar to the test point get more weight than the dissimilar ones. And the key difference occurs only in the way the quantiles are calculated.

What this means is that we can use any conformal prediction technique underneath, but instead of calculating a simple quantile, we need to calculate a weighted one. So, from an implementation perspective, we can look at this class as a wrapper class around other techniques we have defined and convert them into weighted conformal predictions.

Although the weighted conformal prediction can be implemented in many ways, taking in different kinds of weights (across time, across unique_id, and so on), we are going to implement a simpler look-back window-based weighted conformal prediction. We choose the last K timesteps and calculate the weighted quantile on these K scores using the weights given. The weights can either be simple uniform weights across all K steps, or have decayed weights giving the highest weightage to the most recent score. They can even have a completely custom weight.

So, let’s define the __init__ of the class as follows:

class WeightedConformalPredictor:
    def __init__(
        self,
        conformal_predictor: ConformalPrediction,
        K: int,
        weight_strategy: str,
        custom_weights: list = None,
        decay_factor: float = 0.5,
    ):
    …

Here, K is the window, and conformal_predictor is the underlying conformal prediction class we should be using (this should be one of the three classes we have defined). We can define the weight strategy to be either uniform, decay, or custom for uniform weights, decayed weights, or custom weights, respectively. decay_factor decides how fast the decay is applied to the weights for the decayed weighting strategy, and custom_weights lets you specify exactly the weight on these K timesteps.

While we won’t look at the entire code in the text here, we will go through the usual suspects so that you can understand what’s happening. But I would definitely urge you to take some time to digest the code in the file.

First up, we have our fit method. In this method, we just use the score calculation of the underlying conformal predictor and store the calibration dataframe.

def fit(self, Y_calib_df):
    self.calib_df = self.conformal_predictor.calculate_scores(
            Y_calib_df.sort_values(["unique_id", "ds"])
        )

Now, let’s look at the main parts of the predict method.

def predict(self, Y_test_df):
    # Groupby unique_id
    …   
    # Calculate quantiles for each unique_id
    self.q_hat = {}
    for unique_id, group in grouped_calib:
        # Take the last K timesteps
        group = group.iloc[-self.K :]
        scores = group["calib_scores"].values
        # Calculate weights based on the last K timesteps
        total_timesteps = len(scores)
        weights = self._calculate_weight(total_timesteps)
        normalized_weights = weights / weights.sum()
        # Calculate quantile for the current unique_id
        quantile = self.get_weighted_quantile(
                scores, normalized_weights, self.conformal_predictor.alpha
            )
        self.q_hat[unique_id] = quantile
    # Calculate prediction intervals using the underlying conformal predictor's method
    lo, hi = self.conformal_predictor.get_prediction_interval_names()
    Y_test_df[lo], Y_test_df[hi] = (
 self.conformal_predictor.calc_prediction_interval(Y_test_df, self.q_hat)
        )
    return Y_test_df

Let’s now take one of the methods we saw earlier and apply the weighted conformal predictor wrapper on it. For our example, let’s choose the simple ConformalPrediction. Let’s see how we can use this class:

from src.conformal.conformal_predictions import WeightedConformalPredictor
Y_calib_df, Y_test_df = prediction_dict['LSTM']
# Defining an underlying conformal predictor
cp = ConformalPrediction(model="LSTM", level=level)
# using the defined conformal predictor in weighted version
weighted_cp = WeightedConformalPredictor(
    conformal_predictor=cp,
    K=50,
    weight_strategy="uniform",
)
weighted_cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_wcp = weighted_cp.predict(Y_test_df=Y_test_df)

This would create the prediction intervals with the tag CP_Wtd. We can always check the tag by doing weighted_cp.method_name.

Now, there is one small shortcoming with the implementation. Although we are considering the temporal order in the scores, we still have a fixed calibration set. So, until we “re-fit” or calibrate with the latest data points, we will still be working on the same calibrated dataset. So, if you think about it, this should ideally be applied in an online way where each time we predict the new timestep, the previous timestep (with actual value) should be added to the calibration data. We have also provided an alternative implementation that is able to do this in an online fashion. We won’t go into details of the implementation because the core logic is the same, but the API is different, making it possible to update the calibration data. The full implementation can be found in the OnlineWeightedConformalPredictor class in the file.

Let’s see how it can be used. First, we define the setup, initialize the class, and fit the calibration.

from src.conformal.conformal_predictions import OnlineWeightedConformalPredictor
cp = ConformalPrediction(model="LSTM", level=level)
online_weighted_cp = OnlineWeightedConformalPredictor(
    conformal_predictor=cp,
    K=50,
    weight_strategy="uniform",
)
online_weighted_cp.fit(Y_calib_df=Y_calib_df)
joblib.dump(online_weighted_cp, "path/to/saved/file.pkl")

Now, during inference, we can do something like this for each timestep:

# Loading the saved model
online_weighted_cp = joblib.load("path/to/saved/file.pkl")
# current timestep data = current
# past timestep actuals = last_timestep_actuals
prediction = online_weighted_cp.predict_one(current_test)
# updating the calibration data using the last timestep actuals
online_weighted_cp.update(last_timestep_actuals)

For our special case where we are evaluating on test data where we know the actuals, there is another method that does similar online prediction for the data: offline_predict.

Y_test_df_wcpo = online_weighted_cp.offline_predict(Y_test_df=Y_test_df)

Now, let’s look at one last method.

Adaptive Conformal Inference

Lastly, we look at the Adaptive Conformal Inference. This can also be implemented as a wrapper over other conformal prediction methods because this technique involves updating $\alpha$ such that coverage is maintained in spite of the shift in distribution. And because of the nature of the technique, we can only apply this in an online manner, i.e., updating alpha at every time step with the available data. So, this will have the same API as the OnlineWeightedConformalPredictor we saw earlier.

The full class is available in src/conformal/conformal_predictions.py, but here, we will look at some main parts so that you understand it. Let’s take a look at the __init__ function first:

class OnlineAdaptiveConformalInference:
    def __init__(
        self,
        conformal_predictor: ConformalPrediction,
        gamma: float = 0.005,
        update_method: str = "simple",
        momentum_bw: float = 0.95,
        per_unique_id: bool = True,
    ):
        …

Similar to WeightedConformalPredictor, we take in an underlying conformal predictor (conformal_predictor). In addition, we have gamma, which is the step size ($\gamma$), update_method, which can either be simple (taking only the last timestep for update) or momentum (taking the running average of the trajectory of errors). Finally, we can also define momentum_bw, which is the momentum backweight that is used to calculate the weighted average of past errors. A higher momentum (e.g., 0.95) makes the trajectory smoother, having the effect of past miscoverage rates decay slower. Going to the other extreme (0.05) makes the weighted average more reactive and closer to the “simple” method. Lastly, we also have a parameter to do the error calculation for each unique_id separately or pool all errors together.

As usual, we have a fit method that uses a calibration dataset to calculate the scores and keep them ready. We can also see this as a warm-up period in an online implementation. The $\alpha$ update will use the calibration data as an initial history to start off.

def fit(self, Y_calib_df):
    """
    Fit the conformal predictor model with calibration data.
    """
    self.calib_df = self.conformal_predictor.calculate_scores(Y_calib_df)
    self.scores_by_id = (
        self.calib_df.groupby("unique_id")["calib_scores"].apply(list).to_dict()
    )
    # Some more code to initialize necessary data structures     
    …
    return self

Now, let’s see the predict_one method, which predicts one timestep.

def predict_one(self, current_test):
    unique_ids = current_test["unique_id"].unique()
    predictions = []
    for unique_id in unique_ids:
        group_scores = self.scores_by_id.get(unique_id, [])
        if group_scores:
            # Determine the appropriate alpha to use
            alpha = (
                self.alphat[unique_id]
                if self.per_unique_id
                else self.alphat_global
            )
            # Calculate quantile for the current unique_id
            self.q_hat = {
                unique_id: np.quantile(group_scores, 1 - alpha, method="higher")
            }
            # Calculating prediction intervals using conformal_predictor
            …
            current_test[lo] = lower.values
            current_test[hi] = upper.values
            # Storing most recent prediction
            …
        # Collecting and returning concatenated predictions
        …

The code is well-commented and you should be able to understand it. For each unique_id, we get the scores, get the appropriate $\alpha$, calculate the quantile, use this information to calculate the prediction interval using the underlying conformal predictor, and store the prediction for later use.

Now, we look at the update method, which we can use once the actual value for the timestep becomes available.

def update(self, new_data):
    new_scores = self.conformal_predictor.calculate_scores(new_data)
    for unique_id, score in zip(
        new_scores["unique_id"], new_scores["calib_scores"]
    ):
        # Updating score trajectory with new score
    …
        # Retrieve stored predictions and calculate adapt_err
        if self.per_unique_id:
            lower, upper = self.predictions[unique_id]
            actual_y = new_data.loc[new_data["unique_id"] == unique_id, "y"].values[0]
            adapt_err = int(actual_y < lower or actual_y > upper)
            # Update alpha updates the alpha using simple or momentum method
            self.update_alpha(unique_id, adapt_err)
        else:
            # Do the same update at a global error-pooled way

Let’s see how it can be used. First, we define the setup, initialize the class, and fit the calibration.

from src.conformal.conformal_predictions import OnlineAdaptiveConformalInference
cp = ConformalPrediction(model="LSTM", level=level)
aci_cp = OnlineAdaptiveConformalInference(
    conformal_predictor=cp,
    gamma=0.005,
    update_method="simple",
)
aci_cp.fit(Y_calib_df=Y_calib_df)
joblib.dump(aci_cp, "path/to/saved/file.pkl")

Now, during inference, we can do something like below for each timestep:

# Loading the saved model
aci_cp = joblib.load("path/to/saved/file.pkl")
# current timestep data = current
# past timestep actuals = last_timestep_actuals
prediction = aci_cp.predict_one(current_test)
# updating the calibration data using the last timestep actuals
aci_cp.update(last_timestep_actuals)

For our special case where we are evaluating on test data where we know the actuals, there is another method that does similar online prediction for the data: offline_predict.

Y_test_df_aci = aci_cp.offline_predict(Y_test_df=Y_test_df)

Now that we have seen all the techniques in action, let’s take a look at how they have been performing.

Evaluating the results

If you have been following the notebooks, you would know that we have been calculating coverage and average length for each of these techniques. To recap, coverage measures the percentage of time the actual value fell between the intervals we predicted and average length measures the average width of the intervals we predicted. For level=90, we expect the coverage to be around 90% or 0.9 and the average length to be as small as possible, without compromising on coverage.

Let’s take a look at the following figures summarizing the coverage and average length:

Figure 17.18: Coverage for all the conformal techniques (for each unique_id), colored according to how close they are to 0.9 (we had set level=90)

Figure 17.19: Average length for all the conformal techniques (for each unique_id), colored according to how small they are

Here’s a quick recap of the legend to understand the different columns, which are just a combination of these tags based on the application:

• LSTM: Underlying model that we trained to get point prediction
• LSTM_QR: Underlying Quantile Regression model we trained
• LSTM_Normal: Underlying PDF model we trained using Normal distribution assumption
• CP: Regular Conformal Prediction
• CQR: Conformalized Quantile Regression
• CUE: Conformalizing uncertainty Estimates
• CP_Wtd: Weighted corrections over Conformal Prediction
• CP_Wtd_O: Online Weighted corrections over Conformal Prediction
• ACI: Adaptive Conformal Inference

Right off the bat, we can see that the techniques that corrected for distribution shift are performing the best. There are more “greens” on the right side of both tables (closer to 0.9 coverage and lower average lengths). Remember that we started with regular conformal predictions and then corrected them for distributional shifts. We can note that for most unique_id, regular Conformal Prediction has wider than necessary intervals for level = 90. The coverages are greater than 0.9, in most cases 1.0. But both Weighted Correction (CP_Wtd) and Adaptive Conformal Inference (ACI) made the prediction intervals narrower and got closer to our expected level. There is no clear winner between the two and that has to be evaluated on your dataset.

And this concludes our discussion about Probabilistic Forecasts. We hope that with this, you can steer into this lesser-known territory with confidence and deliver more value to whoever you are forecasting for.

Now, let’s take a very high-level look at some niche topics in time series forecasting that get very little attention, but are quite relevant in many domains.

Road less traveled in time series forecasting

In the spirit of Robert Frost’s The Road Not Taken, this section explores the lesser-known, yet highly impactful, techniques in time series forecasting. Just as Frost chooses a path that is less traveled, we delve into niche methods that, while not mainstream, offer unique insights and potential breakthroughs in various domains.

Intermittent or sparse demand forecasting

Intermittent time series forecasting handles data with sporadic, irregular events, often seen in retail, where products may have infrequent sales. It’s crucial for managing inventory and avoiding stockouts or overstocking, especially for slow-moving items. Traditional methods struggle with these patterns because of the fact that for most of these items, the expectation of demand falls to zero, but specialized techniques are needed to improve accuracy, making them essential for retail forecasting.

Here, we will quickly name-drop a few alternative algorithms that were designed for intermittent forecasting and where you can find the implementations for it.

• Croston and its friends: The Croston model, developed in 1972, is used to forecast intermittent demand by separately estimating two components: the demand rate (when sales occur) and the time interval between sales, ignoring periods with zero sales. These estimates are combined to forecast future demand, making it useful for industries with infrequent sales. However, the model doesn’t account for external factors or changes in demand patterns, which can affect its accuracy in more complex situations. Although this was a method developed all the way back in 1972, it has stayed with us in many forms and led to many modifications and extensions. In Nixtla’s statsforecast, we can find CrostonClassic, CrostonOptimized, CrostonSBA, and TSB as four models that can be used just like the other models we saw in Chapter 4, Setting a Strong Baseline Forecast.
• ADIDA: The Aggregate-Disaggregate Intermittent Demand Approach (ADIDA) uses Simple Exponential Smoothing (SES) combined with temporal aggregation to forecast intermittent demand. The method starts by aggregating demand data into non-overlapping time buckets of a size equal to the mean inter-demand interval (MI). SES is then applied to these aggregated values to generate forecasts, which are subsequently disaggregated back to the original time scale, providing demand predictions for each time period. This model is available in statsforecast as ADIDA.
• IMAPA: The Intermittent Multiple Aggregation Prediction Algorithm (IMAPA) forecasts intermittent time series by aggregating values at regular intervals and applying optimized Simple exponential smoothing (SES) to predict future values. IMAPA is efficient, robust to missing data, and effective across various intermittent time series, making it a practical and easy-to-implement choice. This model is available in statsforecast as IMAPA.

Interpretability

We directed you toward a few interpretability techniques for machine learning models back in Chapter 10, Global Forecasting Models. While some of them, such as SHAP and LIME, can still be applied to deep learning models, none of them consider the temporal aspect by design. This is because all those techniques were developed for more general purposes, such as classification and regression. That being said, there has been some work in the interpretability of DL models and time series models. Here, I’ll list a few promising papers that tackle the temporal aspect head-on:

• TimeSHAP: This is a model-agnostic recurrent explainer that builds upon KernelSHAP and extends it to the time series domain. Research paper: https://dl.acm.org/doi/10.1145/3447548.3467166. GitHub: https://github.com/feedzai/timeshap.
• SHAPTime: This technique uses the Shapley Value to provide stable explanations in the temporal dimension, improving forecasting performance. ShapTime’s model-agnostic nature allows it to be deployed on any forecasting model at a lower cost, demonstrating significant performance improvements across various datasets. Research paper: https://link.springer.com/chapter/10.1007/978-3-031-47721-8_45. GitHub: https://github.com/Zhangyuyi-0825/ShapTime.
• Instance-wise Feature Importance in Time (FIT): This is an interpretability technique that relies on the distribution shift between the predictive distribution and a counterfactual, where all but the feature under inspection is unobserved. Research paper: https://proceedings.neurips.cc/paper/2020/file/08fa43588c2571ade19bc0fa5936e028-Paper.pdf. GitHub: https://github.com/sanatonek/time_series_explainability.
• Time Interpret: time_interpret is a library extending Captum, focused on temporal data, providing feature attribution methods to explain predictions from any PyTorch model. It includes synthetic and real-world time series datasets, various PyTorch models, and evaluation methods, with some components also applicable to language model predictions. Research paper: https://arxiv.org/abs/2306.02968. GitHub: https://github.com/josephenguehard/time_interpret.

While this is not an exhaustive list, these are a few works that we feel are important and promising. This is an area of active research, and new techniques will come up as time goes on.

Cold-start forecasting

Cold-start forecasting addresses the challenge of predicting demand for products with no historical data, a common issue in industries like retail, manufacturing, and consumer goods. It arises when launching new products, onboarding brands, or expanding into new regions. Traditional statistical forecasting models like ARIMA or exponential smoothing will not be able to tackle this as they require historical data, a good amount of it as well.

But all hope is not lost. We do have some ways to handle such cases:

• Manual Substitution Mapping: If we know that a new product is substituting another product, we can do manual alignment, which moves the history of the other product to the new product and uses regular techniques to forecast.
• Global Machine Learning Models: We saw how we can model multiple time series using global models in Chapter 10. There, we used some constructed features like lags, rolling aggregations, and so on. But if we train such a model without such features that rely on history and with features characterizing the item (item features that will enable us to cross-learn with other longer time series), we can utilize cross-learning to learn the forecast from similar items. But this works best for new launches that substitute another product.
• Launch Profile Models: If the new product is brand-new, we can use a very similar setup to Global Machine Learning Models but with a small tweak. We can convert our time series into a date/time agnostic manner and consider each time series to start from the date of launch. For instance, the first timestep after the launch of every product can be 1, the next can be 2, and so on. Once we convert all our time series in this manner, we have a dataset that considers how a new item is launched and ramps up. And this can be used to train models that consider the initial ramp-up during launch.
• Foundational Time Series Models: Foundation models for time series leverage large-scale pre-trained models to capture generalized patterns across various time series tasks. These models are highly effective in applications like cold start forecasting, where little or no historical data is available. By using foundational models, practitioners can apply pre-trained knowledge to new scenarios, improving accuracy in forecasting tasks across industries such as retail, healthcare, and finance. The models offer flexibility, allowing fine-tuning or zero-shot applications, making them particularly useful for complex, sparse, or intermittent data. Many foundational models for time series are already available for use—some are commercial and some others are open source. This is a nice survey of the foundational models at the time of writing the book: Foundation Models for Time Series Analysis: A Tutorial and Survey: https://arxiv.org/abs/2403.14735. The performance of such models is average at best, but in the absence of any data (cold-start), this presents as a good option to try out. A few popular ways of doing this practically are:
  • TimeGPT (Commercial) from Nixtla: https://nixtlaverse.nixtla.io/nixtla/docs/getting-started/introduction.html
  • Chronos from Amazon: https://auto.gluon.ai/stable/tutorials/timeseries/forecasting-chronos.html
  • Moirai from Salesforce: https://huggingface.co/Salesforce/moirai-1.0-R-large
  • TimesFM from Google: https://github.com/google-research/timesfm

Hierarchical forecasting

Hierarchical forecasting deals with time series that can be disaggregated into nested levels, such as product categories or geographic regions. These structures require forecasts that maintain coherence, meaning predictions at lower levels should sum up to higher levels, reflecting the aggregation. Grouped time series, which combine different hierarchies (e.g., product type and geography), add complexity. The goal is to produce consistent, accurate forecasts that align with the data’s natural aggregation, making hierarchical forecasting essential for businesses managing large collections of time series across multiple dimensions.

There are some techniques of disaggregating, aggregating, or reconciling all forecasts so that they add up in a logical manner. Chapter 10 from Rob Hyndman’s free Bible of time series forecasting (Forecasting: Principles and Practice) talks about this at length and is a very good resource to get up to speed (fair warning: it’s quite math-intensive). You can find the chapter here: https://otexts.com/fpp2/hierarchical.html. For a more practical approach, you can check out Nixtla’s hierarchicalforecast library here: https://nixtlaverse.nixtla.io/hierarchicalforecast/index.html.

And with that, one of the longest chapters of the book comes to an end. Congrats on making it and digesting all the information. Do feel free to use the notebooks and play around with the code, different options, and so on to get a better grasp of what’s happening.

Summary

In this chapter, we explored different techniques for generating probabilistic forecasts like Probability Density Functions, Quantile functions, Monte Carlo Dropout, and Conformal Prediction. We deep-dived into each of them and learned how to apply them to real data. For Conformal Predictions, an actively researched field, we learned different ways to conformalize prediction intervals with different underlying mechanisms, like Conformalized Quantile Regression, conformalizing uncertainty estimates, and so on. And finally, we saw a few tweaks to make conformal prediction work even better for time series problems.

To wrap it up, we also looked at a few topics that are on the road less traveled, like intermittent demand forecasting, interpretability, cold-start forecasting, and hierarchical forecasting.

In the next part of this book, we will look at a few mechanics of forecasting, such as multi-step forecasting, cross-validation, and evaluation.

References

The following are the references for this chapter:

1. Tony Duan, A. Avati, D. Ding, S. Basu, A. Ng, and Alejandro Schuler. (2019). NGBoost: Natural Gradient Boosting for Probabilistic Prediction. International Conference on Machine Learning. https://proceedings.mlr.press/v119/duan20a/duan20a.pdf.
2. Y. Gal and Zoubin Ghahramani. (2015). Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning. International Conference on Machine Learning. https://proceedings.mlr.press/v48/gal16.html.
3. Valentin Flunkert, David Salinas, and Jan Gasthaus. (2017). DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks. International Journal of Forecasting. https://www.sciencedirect.com/science/article/pii/S0169207019301888.
4. Koenker, Roger. (2005). Quantile Regression. Cambridge University Press. pp. 146–7. ISBN 978-0-521-60827-5. http://www.econ.uiuc.edu/~roger/research/rq/QRJEP.pdf.
5. Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos. (2020). The M4 Competition: 100,000 time series and 61 forecasting methods. International Journal of Forecasting. https://www.sciencedirect.com/science/article/pii/S0169207019301128.
6. Loic Le Folgoc and Vasileios Baltatzis and Sujal Desai and Anand Devaraj and Sam Ellis and Octavio E. Martinez Manzanera and Arjun Nair and Huaqi Qiu and Julia Schnabel and Ben Glocker. (2021). Is MC Dropout Bayesian?. arXiv preprint arXiv: Arxiv-2110.04286. https://arxiv.org/abs/2110.04286.
7. Nicolas Dewolf and Bernard De Baets and Willem Waegeman. (2023). Valid prediction intervals for regression problems. Artif. Intell. Rev. https://doi.org/10.1007/s10462-022-10178-5.
8. V. Vovk, A. Gammerman, and G. Shafer. (2005). Algorithmic learning in a random world. Springer. https://link.springer.com/book/10.1007/b106715.
9. Anastasios N. Angelopoulos and Stephen Bates. (2021). A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification. arXiv preprint arXiv: Arxiv-2107.07511. https://arxiv.org/abs/2107.07511.
10. Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. (2002). Inductive confidence machines for regression. Machine Learning: ECML 2002. ECML 2002. Lecture Notes in Computer Science(), vol 2430. Springer, Berlin, Heidelberg. https://doi.org/10.1007/3-540-36755-1_29.
11. Romano, Yaniv and Patterson, Evan and Candes, Emmanuel. (2019). Conformalized Quantile Regression. Advances in Neural Information Processing Systems. https://proceedings.neurips.cc/paper_files/paper/2019/file/5103c3584b063c431bd1268e9b5e76fb-Paper.pdf.
12. R. Barber, E. Candès, Aaditya Ramdas, and R. Tibshirani. (2022). Conformal prediction beyond exchangeability. Annals of Statistics. https://projecteuclid.org/journals/annals-of-statistics/volume-51/issue-2/Conformal-prediction-beyond-exchangeability/10.1214/23-AOS2276.full.
13. Gibbs, Isaac and Candes, Emmanuel. (2021). Adaptive Conformal Inference Under Distribution Shift. Advances in Neural Information Processing Systems. https://proceedings.neurips.cc/paper_files/paper/2021/file/0d441de75945e5acbc865406fc9a2559-Paper.pdf.

Further reading

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

• The Gradient Boosters VI(B): NGBoost: https://deep-and-shallow.com/2020/06/27/the-gradient-boosters-vib-ngboost/
• Dive into Deep Learning, Chapter 5.6: https://d2l.ai/chapter_multilayer-perceptrons/dropout.html
• Bayesian Inference by Marco Taboga: https://www.statlect.com/fundamentals-of-statistics/Bayesian-inference
• Seeing Theory: Bayesian Inference: https://seeing-theory.brown.edu/bayesian-inference/index.html
• A Tutorial on Sparse Gaussian Processes and Variational Inference by Felix Leibfried et al.: https://arxiv.org/pdf/2012.13962
• A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification by Anastasios N. Angelopoulos and Stephen Bates. (2021): https://arxiv.org/abs/2107.07511

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