4

Setting a Strong Baseline Forecast

In the previous chapter, we saw some techniques we can use to understand time series data, do some Exploratory Data Analysis (EDA), and so on. But now, let’s get to the crux of the matter—time series forecasting. The point of understanding the dataset and looking at patterns, seasonality, and so on was to make the job of forecasting that series easier. And with any machine learning exercise, one of the first things we need to establish before going further is a baseline.

A baseline is a simple model that provides reasonable results without requiring a lot of time to come up with them. Many people think of a baseline as something that is derived from common sense, such as an average or some rule of thumb. But as a best practice, a baseline can be as sophisticated as we want it to be, so long as it is quickly and easily implemented. Any further progress we want to make will be in terms of the performance of this baseline.

In this chapter, we will look at a few classical techniques that can be used as baselines, and strong baselines at that. Some may feel that the forecasting techniques we will be discussing in this chapter shouldn’t be baselines, but we are keeping them in here because these techniques have stood the test of time—and for good reason. They are also very mature and can be applied with very little effort, thanks to the awesome open source libraries that implement them. There can be many types of problems/datasets where it is difficult to beat the baseline techniques we will discuss in this chapter, and in those cases, there is no shame in just sticking to one of these baseline techniques.

In this chapter, we will cover the following topics:

Technical requirements

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

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

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

Setting up a test harness

Before we start forecasting and setting up baselines, we need to set up a test harness. In software testing, a test harness is a collection of code and inputs that have been configured to test a program under various situations. In terms of machine learning, a test harness is a set of code and data that can be used to evaluate algorithms. It is important to set up a test harness so that we can evaluate all future algorithms in a standard and quick way.

The first thing we need is holdout (test) and validation datasets.

Creating holdout (test) and validation datasets

As a standard practice, in machine learning, we set aside two parts of the dataset, name them validation data and test data, and don’t use them at all to train the model. The validation data is used in the modeling process to assess the quality of the model. To select between different model classes, tune the hyperparameters, perform feature selection, and so on, we need a dataset. Test data is like the final test of your chosen model. It tells you how well your model is doing in unseen data. If validation data is like the mid-term exams, the test data is your final exam.

In regular regression or classification, we usually sample a few records at random and set them aside. But while dealing with time series, we need to respect the temporal aspect of the dataset. Therefore, a best practice is to set aside the latest part of the dataset as the test data. Another rule of thumb is to set equal-sized validation and test datasets so that the key modeling decisions we make based on the validation data are as close as possible to the test data. The dataset that we introduced in Chapter 2, Acquiring and Processing Time Series Data, the London Smart Energy dataset, contains the energy consumption readings of households in London from November 2011 to February 2014. So, we are going to put aside January 2014 as the validation data and February 2014 as the test data.

Let’s open 01-Setting_up_Experiment_Harness.ipynb from the Chapter04 folder and run it. In the notebook, we must create the train-test split both before and after filling the missing values with SeasonalInterpolation and save them accordingly. Once the notebook finishes running, you will have created the following files in the pre-processed folder with the 2014 data saved separately:

Now that we have a fixed dataset that can be used to fairly evaluate multiple algorithms, we need a way to evaluate the different forecasts.

Choosing an evaluation metric

In machine learning, we have a handful of metrics that can be used to measure continuous outputs, mainly Mean Absolute Error and Mean Squared Error. But in the time series forecasting realm, there are scores of metrics with no real consensus on which ones to use. One of the reasons for this overwhelming number of metrics is that no one metric measures every characteristic of a forecast. Therefore, we have a whole chapter devoted to this topic (Chapter 19, Evaluating Forecast Errors—A Survey of Forecast Metrics). For now, we will just review a few metrics, all of which we are going to use to measure the forecasts. We are just going to consider them at face value:

$$MAE = \frac{1}{N \times L} \times \sum_{i}^{N} \sum_{j}^{L} |f_{i,j} - y_{i,j}|$$

Here, N is the number of time series, L is the length of time series (in this case, the length of the test period), and f and y are the forecast and observed values, respectively.

$$MSE = \frac{1}{N \times L} \times \sum_{i}^{N} \sum_{j}^{L} (f_{i,j} - y_{i,j})^2$$

$$MASE = \frac{\frac{1}{L} \times \sum_{i}^{L} |f_i - y_i|}{\frac{1}{L-1} \times \sum_{j=2}^{L} |y_j - y_{j-1}|}$$

We calculate forecast bias as the difference between the sum of the forecast and the sum of the observed values, expressed as a percentage over the sum of all actuals:

$$FB = \frac{\sum_{i}^{N} \sum_{j}^{L} f_{i,j} - \sum_{i}^{N} \sum_{j}^{L} y_{i,j}}{\sum_{i}^{N} \sum_{j}^{L} y_{i,j}}$$

Now, our test harness is ready. We also know how to evaluate and compare forecasts that have been generated from different models on a single, fixed holdout dataset with a set of predetermined metrics. Now, it’s time to start forecasting.

Generating strong baseline forecasts

Time series forecasting has been around since the early 1920s, and through the years, many brilliant people have come up with different models, some statistical and some heuristic-based. I refer to them collectively as classical statistical models or econometrics models, although they are not strictly statistical/econometric.

In this section, we are going to review a few such models that can form really strong baselines when we want to try modern techniques in forecasting. As an exercise, we are going to use an excellent open source library for time series forecasting—NIXTLA (https://github.com/Nixtla). The 02-Baseline_Forecasts_using_NIXTLA.ipynb notebook contains the code for this section so that you can follow along.

Before we start looking at forecasting techniques, let’s quickly understand how to use the NIXTLA library to generate forecasts. We are going to pick one consumer from the dataset and try out all the baseline techniques on the validation dataset one by one.

The first thing we need to do is select the consumer we want using the unique ID for each customer, the LCLid column (from the expanded form of data), and set the timestamp as the index of the DataFrame:

ts_train = train_df.loc[train_df.LCLid=="MAC000193",['LCLid',"timestamp","energy_consumption"]]
ts_val = val_df.loc[val_df.LCLid=="MAC000193", ['LCLid',"timestamp","energy_consumption"]]
ts_test = test_df.loc[test_df.LCLid=="MAC000193", ['LCLid',"timestamp","energy_consumption"]]

NIXTLA has the flexibility to work directly with either pandas or Polars DataFrames. By default, NIXTLA looks for three columns:

This is very convenient as there is no need for further manipulation to go from data to modeling. NIXTLA follows the scikit-learn style with .fit() and .predict() and also adopts a .forecast() method, which is a memory-efficient method that doesn’t store the partial model outputs, whereas the scikit-learn interface stores the fitted models:

sf = StatsForecast(
    models=[model],
    freq=freq,
    n_jobs=-1,
    fallback_model=Naive()
)
sf.fit(df = _ts_train,          id_col = 'LCLid',
       time_col = 'timestamp',
       target_col = 'energy_consumption',
)
baseline_test_pred_df = sf.predict(len(ts_test) )

NIXTLA also has a .forecast() method, which is a memory-efficient method that doesn’t store the partial model outputs, whereas the scikit-learn interface stores the fitted models:

# Efficiently fit and predict without storing memory
y_pred = sf.forecast(
    h=len(ts_test),
    df=ts_train,
    id_col = 'LCLid',    time_col = 'timestamp',    target_col = 'energy_consumption',
)

When we call .predict or .forecast, we have to tell the model how long into the future we have to predict. This is called the horizon of the forecast. In our case, we need to predict our test period, which we can easily do by just taking the length of the ts_test array.

We can also calculate the metrics we discussed earlier in the test harness easily using NIXTLA’s classes. For added flexibility, we can loop through a list of metrics to get multiple measurements for each forecast:

# Calculate metrics
metrics = [mase, mae, mse, rmse, smape, forecast_bias]
for metric in metrics:
    metric_name = metric.__name__
    if metric_name == 'mase':
        evaluation[metric_name] = 	
metric(results[target_col].values,            results[model_name].values,
ts_train[target_col].values, seasonality=48)
    else:
        evaluation[metric_name] =
metric(results[target_col].values,
results[model_name].values)

Notice that, for MASE, the training set is also included.

For ease of experimentation, we have encapsulated all of this into a handy function, evaluate_performance, in the notebook. This returns the predictions and the calculated metrics in a DataFrame.

Now, let’s start looking at a few very simple methods of forecasting.

Naïve forecast

A naïve forecast is as simple as you can get. The forecast is just the last/most recent observation in a time series. If the latest observation in a time series is 10, then the forecast for all future timesteps is 10. This can be implemented as follows using the Naive class in NIXTLA:

from statsforecast.models import Naive
models = Naive()

Once we have initialized the model, we can call our helpful evaluate_performance function in the notebook to run and record the forecast and metrics.

Let’s visualize the forecast we just generated:

Figure 4.1: Naïve forecast

Here, we can see that the forecast is a straight line and completely ignores any pattern in the series. This is by far the simplest way to forecast, hence why it is naïve. Now, let’s look at another simple method.

Moving average forecast

While a naïve forecast memorizes the most recent past, it also memorizes the noise at any timestep. A moving average forecast is another simple method that tries to overcome the pure memorization of the naïve method. Instead of taking the latest observation, it takes the mean of the latest n steps as the forecast. Moving average is not one of the models present in NIXTLA, but we have implemented a NIXTLA-compatible model in this book’s GitHub repository in the Chapter04 folder:

from src.forecasting.baselines import NaiveMovingAverage
#Taking a moving average over 48 timesteps, i.e, one day
naive_model = NaiveMovingAverage(window=48)

Let’s look at the forecast we generated:

Figure 4.2 – Moving average forecast

Figure 4.2: Moving average forecast

This forecast is also almost a straight line. Now, let’s look at another simple method, but one that considers seasonality as well.

Seasonal naive forecast

A seasonal naive forecast is a twist on the simple naive method. In the naive method, we took the last observation (Yt-1), whereas in seasonal naïve, we take the Yt-k observation. So, we look back k steps for each forecast. This enables the algorithm to mimic the last seasonality cycle. For instance, if we set k=48*7, we will be able to mimic the latest seasonal weekly cycle.

This method is implemented in NIXTLA and we can use it like so:

from statsforecast.models import SeasonalNaive
seasonal_naive = SeasonalNaive(season_length=48*7)

Let’s see what this forecast looks like:

Figure 4.3: Seasonal naïve forecast

Here, we can see that the forecast is trying to mimic the seasonality pattern. However, it’s not very accurate because it is blindly following the last seasonal cycle.

Now that we’ve looked at a few simple methods, let’s look at a few statistical models.

Exponential smoothing

Exponential smoothing is one of the most popular methods for generating forecasts. It has been around since the late 1950s and has proved its mettle and stood the test of time. There are a few different variants of ETS—single exponential smoothing, double exponential smoothing, Holt-Winters’ seasonal smoothing, and so on. But all of them have one key idea that has been used in different ways. In the naïve method, we were just using the latest observation, which is like saying only the most recent data point in history matters and no data point before that matters. On the other hand, the moving average method considers the last n observations to be equally important and takes the mean of them.

ETS combines both these intuitions and says that all the history is important, but the recent history is more important. Therefore, the forecast is generated using a weighted average where the weights decrease exponentially as we move farther into the history:

$$f_t = \alpha \cdot y_{t-1} + \alpha \cdot (1 - \alpha) \cdot y_{t-2} + \alpha \cdot (1 - \alpha)^2 \cdot y_{t-3} + \cdots$$

Here, $0 \leq \alpha \leq 1$ is the smoothing parameter that lets us decide how fast or slow the weights should decay, yt is the actuals at timestep t, and ft is the forecast at timestep t.

Simple exponential smoothing (SES) is when you simply apply this smoothing procedure to the history. This is more suited for time series that have no trends or seasonality, and the forecast is going to be a flat line. The forecast is generated using the following formula:

Forecast Equation$f_t = \alpha y_{t-1} + (1 - \alpha) f_{t-1}$

Double exponential smoothing (DES) extends the smoothing idea to model trends as well. It has two smoothing equations—one for the level and the other for the trend. Once you have the estimate of the level and trend, you can combine them. This forecast is not necessarily flat because the estimated trend is used to extrapolate it into the future. The forecast is generated according to the following formula:

Forecast Equation$f_{t+h} = l_t + h \cdot b_t$
Level Equation$l_t = \alpha \cdot y_t + (1 - \alpha) \cdot (l_{t-1} + b_{t-1})$
Trend Equation$b_t = \beta \cdot (l_t - l_{t-1}) + (1 - \beta) \cdot b_{t-1}$

First, we estimate the level (lt) using the Level Equation with the available observations. Then, we estimate the trend using the Trend Equation. Finally, to get the forecast, we combine lt and bt using the Forecast Equation.

Researchers have found empirical evidence that this kind of constant extrapolation can result in over-forecasts over the long-term forecast. This is because, in the real world, time series data doesn’t increase at a constant rate forever. Motivated by this, an addition to this has also been introduced that dampens the trend by a factor of $0 \lt \phi \lt 1$, such that when $\phi = 1$, there is no damping, and it is identical to DES.

Triple exponential smoothing or Holt-Winters’ (HW) takes this one step forward by including another smoothing term to model the seasonality. This has three parameters ( $\alpha$, $\beta$, $\gamma$) for the smoothing and uses a seasonality period (m) as input parameters. You can also choose between additive or multiplicative seasonality. The forecast equations for the additive model are as follows:

Forecast Equation$f_{t+h} = l_t + h \cdot b_t + s_{t+h-m(k+1)}$
Level Equation$l_t = \alpha \cdot y_t - s_{t-m} + (1 - \alpha) \cdot (l_{t-1} + b_{t-1})$
Trend Equation$b_t = \beta \cdot (l_t - l_{t-1}) + (1 - \beta) \cdot b_{t-1}$
Seasonality Equation$s_t = \gamma (y_t - l_{t-1} - b_{t-1}) + (1 - \gamma) s_{t-m}$

These formulae are also used like in the double exponential case. Instead of estimating level and trend, we estimate level, trend, and seasonality separately.

The family of ETS methods is not limited to the three that we just discussed. A way to think about the different models is in terms of the trend and seasonal components of these models. The trend can either be no trend, additive, or additive damped. The seasonality can be no seasonality, additive, or multiplicative. Every combination of these parameters is a different technique in the family, as shown in the following table:

Trend component

Seasonal component

N (None)

A (Additive)

M (Multiplicative)

N (None)

Simple Exponential Smoothing

-

-

A (Additive)

Double Exponential Smoothing

Additive Holt-Winters

Multiplicative Holt-Winters

Ad (Additive damped)

Damped Double Exponential Smoothing

-

Damped Holt-Winters

Table 4.1: Exponential smoothing family

NIXTLA has an entire family of ETS methods.

Let’s see how we can initialize the ETS model in NIXTLA:

from statsforecast.models import (SimpleExponentialSmoothing, Holt, HoltWinters, AutoETS)
exp_smooth = HoltWinters(error_type = 'A', season_length = 48)]

Here, error_type = 'A' refers to additive error. The user has the option for either additive error or multiplicative error, which could be called using error_type = 'M'. NIXTLA models have an option to use AutoETS(). This model will automatically choose which exponential smoothing model is the best option: simple exponential smoothing, double exponential smoothing (Holt’s method), or triple exponential smoothing (Holt-Winters method). It will also choose which parameters and error types are best for each individual time series. Refer to the GitHub notebooks for examples of how to use AutoETS().

Let’s see what the forecast using ETS looks like in Figure 4.4:

Figure 4.4: Exponential smoothing forecast

The forecast has captured the seasonality but has failed to capture the peaks. But we can see the improvement in MAE already.

Now, let’s look at one of the most popular forecasting methods out there.

AutoRegressive Integrated Moving Average (ARIMA)

ARIMA models are the other class of methods that, like ETS, have stood the test of time and are one of the most popular classical methods of forecasting. The ETS family of methods is modeled around trend and seasonality, while ARIMA relies on autocorrelation (the correlation of yt with yt-1, yt-2, and so on).

The simplest in the family are the AR (p) models, which use linear regression with p previous timesteps or, in other words, p lags. Mathematically, it can be written as follows:

$$y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \epsilon_t$$

Here, c is the intercept, and $\epsilon_t$ is the noise or error at timestep t.

The next in the family are MA (q) models, in which, instead of past observed values, we use the past q errors in the forecast (which is assumed to be pure white noise) to come up with a forecast:

$$y_t = c + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \cdots + \theta_q \epsilon_{t-q}$$

Here, $\epsilon$ is white noise and c is the intercept.

This is not typically used on its own but in conjunction with AR (p) models, which makes the next one on our list ARMA (p, q) models. ARMA (AutoRegressive Moving Average) models are defined as yt = AR (p) + MA (q).

In all the ARIMA models, there is one underlying assumption—the time series is stationary (we talked about stationarity in Chapter 1, Introducing Time Series, and will elaborate on this in Chapter 6, Feature Engineering for Time Series Forecasting). There are many ways to make the series stationary but taking the difference of successive values is one such technique. This is known as differencing. Sometimes, we need to do differencing once, while other times, we have to perform successive differencing before the time series becomes stationary. The number of times we do the differencing operation is called the order of differencing. The I in ARIMA, and the final piece of the puzzle, stands for Integrated. It defines the order of differencing we need to do before the series becomes stationary and is denoted by d.

So, the complete ARIMA (p, d, q) model says that we do the dth order of differencing and then consider the last p terms in an autoregressive manner, and then include the last q moving average terms to come up with the forecast.

The ARIMA models we have discussed so far only handle non-seasonal time series. However, using the same concepts we discussed, but on a seasonal cycle, we get seasonal ARIMA. p, d, and q are slightly tweaked so that they work on the seasonal period, m. To differentiate them from the normal p, d, and q, we call the seasonal values P, D, and Q. For instance, if p means taking the last p lags, P means taking the last P seasonal lags. If p1 is yt-1, P1 would be yt-m. Similarly, D means the order of seasonal differencing.

Picking the right p, d, and q and P, D, and Q values is not very intuitive, and we will have to resort to statistical tests to find them. However, this becomes a bit impractical when you are forecasting many time series. An automatic way of iterating through the different parameters and finding the best p, d, and q, and P, D, and Q values for the data is called AutoARIMA. In Python, NIXTLA has implemented this method, AutoARIMA(). NIXTLA also has a normal ARIMA implementation as well, which is much faster but requires p, d, and q to be entered manually.

Practical considerations:

Although ARIMA and AutoARIMA can give you good-performing models in many cases, they can be quite slow when you have long seasonal periods and a long time series. In our case, where we have almost 27K observations in the history, ARIMA becomes very slow and a memory hog. Even when subsetting the data, a single AutoARIMA fit takes around 60 minutes. Letting go of the seasonal parameters brings down the runtime drastically, but for a seasonal time series such as energy consumption, it doesn’t make sense. AutoARIMA includes many such fits to identify the best parameters and, therefore, becomes impractical for long time series datasets. Almost all the implementations in the Python ecosystem suffer from this drawback. NIXLTA claims to have the fastest and most accurate version of AutoARIMA, faster than the original R method as well.

Let’s see how we can apply ARIMA and AutoARIMA using NIXTLA:

from statsforecast.models import (ARIMA, AutoARIMA)
#ARIMA model by specifying parameters
arima_model = ARIMA(order = (2,1,1), seasonal_order = (1,1,1), season_length = 48)
#AutoARIMA model by specifying max limits for parameters and letting the algorithm find the best ones
auto_arima_model = AutoARIMA( max_p = 2, max_d=1, max_q = 2, max_P=2, max_D = 1, max_Q = 2, stepwise = True, season_length=48)

For the entire list of parameters for AutoARIMA, head over to the NIXTLA documentation at https://nixtlaverse.nixtla.io/statsforecast/docs/models/autoarima.html.

Let’s see what the ETS and ARIMA forecasts look like for the households we were experimenting with:

Figure 4.5: ETS and ARIMA forecasts

With NIXTLA, both ETS and ARIMA have done a good job of capturing both the seasonality and the peaks. The resulting MAE scores are also very similar, with 0.191 and 0.203, respectively. Now, let’s look at another method—the Theta forecast.

Theta forecast

The Theta forecast was the top-performing submission in the M3 forecasting competition that was held in 2002. The method relies on a parameter, $\theta$, that amplifies or smooths the local curvature of a time series, depending on the value chosen. Using $\theta$, we smooth or amplify the original time series. These smoothed lines are called Theta lines. V. Assimakopoulos and K. Nikolopoulos proposed this method as a decomposition approach to forecasting. Although, in theory, any number of Theta lines can be used, the originally proposed method used two Theta lines, $\theta = 0$ and $\theta = 2$, and took an average of the forecast of the two Theta lines as the final forecast.

The M competitions are forecasting competitions organized by Spyros Makridakis, a leading forecasting researcher. They typically curate a dataset of time series, lay down the metrics with which the forecasts will be evaluated, and open these competitions to researchers all around the world to get the best forecast possible. These competitions are considered to be some of the biggest and most popular time series forecasting competitions in the world. At the time of writing, six such competitions have already been completed. To learn more about the latest competition, visit this website: https://mofc.unic.ac.cy/the-m6-competition/.

In 2002, Rob Hyndman et al. simplified the Theta method and showed that we can use ETS with a drift term to get equivalent results to the original Theta method, which is what is adapted into most of the implementations of the method that exist today. The main steps that are involved in the Theta forecast (which is implemented in NIXTLA) are as follows:

  1. Deseasonalization: Apply a classical multiplicative decomposition to remove the seasonal component from the time series (if it exists). This focuses the analysis on the underlying trend and cyclical components. Deseasonalization is done using statsmodels.tsa.seasonal.seasonal_decompose. This step creates a new deseasonalized time series.
  2. Theta Coefficients Application: Decompose the deseasonalized series into two “Theta” lines using coefficients $\theta_1$ and $\theta_2$. These coefficients modify the second difference of the time series to either dampen ($\theta \lt 1)$ or accentuate $(\theta \gt 1)$ local fluctuations.
  3. Extrapolation of Theta Lines: Treat each Theta line as a separate series and forecast them into the future. This is done using linear regression for the Theta line where $\theta_1 = 0$, producing a straight line, and simple exponential smoothing for the Theta line where $\theta_2 = 2$.
  4. Recomposition: Combine the forecasts from the two Theta lines. The original method uses equal weighting for both lines, which integrates the long-term trend and short-term movements effectively.
  5. Reseasonalize: If the data was deseasonalized in the beginning.

NIXTLA has very different variations of the Theta method. More information on the specifics of the NIXTLA implementation can be found here: https://nixtlaverse.nixtla.io/statsforecast/docs/models/autotheta.html.

Let’s see how we can use it practically:

theta_model = Theta(season_length =48, decomposition_type = 'additive' )

The key parameters here are as follows: season_length and decomposition_type. These parameters are used for the initial seasonal decomposition. If left empty, the implementation automatically tests for seasonality and deseasonalizes the time series automatically using multiplicative decomposition. It is recommended to set these parameters with our domain knowledge if we know them. The decomposition type can be multiplicative (default) or additive.

Let’s visualize the forecast we just generated using the Theta forecast:

Figure 4.6: The Theta forecast

Reference check:

The research paper in which V. Assimakopoulos and K. Nikolopoulos proposed the Theta method is cited as reference 1 in the References section, while subsequent simplification by Rob Hyndman is cited as reference 2.

The seasonality pattern is captured, but it’s not hitting the peaks. Let’s look at another very strong method, TBATS.

TBATS

Sometimes, a time series has more than one seasonality pattern or a non-integer seasonal period, commonly referred to as complex seasonality. An example would be an hourly forecast that could have a daily seasonality for the time of day, a weekly seasonality for the day of the week, and a yearly seasonality for the day of the year. Additionally, most time series models are designed for smaller integer seasonal periods, such as monthly (12) or quarterly (4) data, but yearly seasonality can pose a problem since a year is 364.25 days. TBATS was meant to combat these many challenges that pose problems for many forecasting models. However, with any automated approach, at times it is susceptible to poor forecasts.

TBATS stands for:

This model was first introduced by Rob J. Hyndman, Alysha M. De Livera, and Ralph D. Snyder in 2011. There is also another variant of TBATS, referred to as BATS, which is without the trigonometric seasonality component. TBATS is from the state space model family. In state space forecasting models, the observed time series is assumed to be a combination of the underlying state variables and a measurement equation that relates the state variables to the observed data. The state variables capture the underlying patterns, trends, and relationships in the data.

BATS has parameters $(\omega, \phi, p, q, m_1, m_2, \ldots m_t)$ indicating the Box-Cox parameter, damping parameter, ARMA parameters (p, q) and the seasonal periods (m1, m2, …, mt). Due to its flexibility, the BATS model can be considered a family of models encompassing many other models we have seen earlier. For instance:

BATS has the flexibility for multiple seasonality; however, it is limited to only integer-based seasonal periods, and with multiple seasonalities, it can have a large number of states resulting in increasing model complexity. This is what TBATS was meant to address.

For reference, the TBATS parameter space is:

$$TBAT(\omega, \phi, p, q, \{m_1, k_1\}, ..., \{m_T, k_T\})$$

The main advantages of TBATS are as follows:

To better understand the inner workings of TBATS, let’s break down each step.

The order in which operations are done (unlike the order in the acronym) using TBATS is:

  1. Box-Cox transformation
  2. Exponentially smoothed trend
  3. Seasonal decomposition using Fourier series (trigonometric seasonality)
  4. AutoRegressive Moving Average (ARMA)
  5. Parameter estimation through a likelihood-based approach

Box-Cox transformation

Box-Cox is a transformation in the family of power transformations.

In time series, making data stationary is an important step before forecasting (as discussed in Chapter 1). Stationarity ensures that our data does not statistically change over time, and thus more accurately resembles a probability distribution. There are several possible transformations that could be applied. More details on various target transformations, including Box-Cox, can be found in Chapter 7.

As a preview, here is a sample output from a Box-Cox transformation. After the transformation, our data more closely resembles that of a normal distribution. Box-Cox transformations can only be used with positive data, but in practice, this is often the case.

Figure 4.7 shows an example of how a time series might look before and after a Box-Cox transformation.

Figure 4.7: Box-Cox transformation

Exponentially smoothed trend

Using Locally Estimated Scatterplot Smoothing (LOESS), a smoothed trend is extracted from the time series:

$$l_t = l_{t-1} + \phi b_{t-1} + \alpha d_t \; (Local)$$

$$b_t = (1 - \phi) b + \phi b_{t-1} + \beta d_t \; (Global)$$

LOESS works by applying a locally weighted, low-degree polynomial regression over the data points to create a smooth, flowing line through them. This technique is highly effective in capturing local trend variations without assuming a global form for the data, which makes it particularly useful for data with varying trends or seasonal variations. This is the same LOESS that we used to decompose a time series into a trend back in Chapter 3.

Seasonal decomposition using Fourier series (trigonometric seasonality)

The remaining residuals are then modeled using Fourier terms (discussed in Chapter 3) to decompose the seasonality component.

$$y_t^{(\omega)} = l_{t-1} + \phi b_{t-1} + \sum_{i=1}^{M} S_{t-m_i}^{(i)} + d_t$$

$$S_t^{(i)} = \sum_{j=1}^{k_i} S_{j,t}^{(i)}$$

$$S_{j,t}^{(i)} = S_{j,t-1}^{(i)} \cos \lambda_j^{(i)} + S_{j,t-1}^{*(i)} \sin \lambda_j^{(i)} + \gamma_1^{(i)} d_t$$

$$S_{j,t}^{*(i)} = -S_{j,t-1}^{(i)} \sin \lambda_j^{(i)} + S_{j,t-1}^{*(i)} \cos \lambda_j^{(i)} + \gamma_2^{(i)} d_t$$

The main advantage of using Fourier to model seasonality is its ability to model multiple seasonalities, as well as non-integer seasonality, such as yearly seasonality with daily data since there are 364.25 days in a year. Most other decomposition methods cannot handle the non-integer period and have to resort to rounding to 365, which can fail to identify the true seasonality. An example of what a decomposed time series would like using Fourier is below. The observed time series in this example is hourly data. Therefore, our seasonal periods are:

daily = 24

weekly = 24 * 7 = 168

Here, you can clearly see the defined seasonal patterns, the trend, and the remaining residuals. Figure 4.8 shows the decomposition of the trend and seasonality, after which the residuals are modeled using an ARMA process.

Figure 4.8: Decomposed time series

ARMA

ARMA was discussed earlier as a subset of the ARIMA family:

$$d_t = \sum_{i=1}^{p} \phi_i d_{t-i} + \sum_{j=1}^{q} \phi_j \varepsilon_{t-j} + \varepsilon_t$$

The ARMA model in TBATS is used to model the remaining residuals to capture any autocorrelations of the lagged variables. The autoregressive (AR) component captures the correlation between an observation and several lagged observations. This deals with the momentum or continuation of the series. The moving average (MA) component models the error terms as a linear combination of errors at previous time periods, capturing information not explained by the AR part alone.

Parameter optimization

To select the optimal parameter space, TBATS will fit several models and automatically select the best parameters. A few of the models TBATS fits internally are:

The final model is chosen by which combination of parameters minimizes the Akaike Information Criterion (AIC), and AutoARIMA is used to determine the ARMA parameters.

As with all forecasting methods, there are benefits and trade-offs to different models. While TBATS offers some enhancements on many other models’ shortcomings, the trade-off is the need to build many models, which results in longer computation times. This can pose a problem if you have to model multiple time series. Additionally, TBATS does not allow for the inclusion of exogenous variables.

Practitioner’s note:

TBATS cannot handle exogenous regression since it is related to ETS models, as per Hyndman himself, who suggests it is unlikely to include covariates (Hyndman, 2014; Reference 7). If external regressors are to be used, other methods such as ARIMAX or SARIMAX should be used. If the time series has complex seasonality, you can add Fourier features as covariates to your ARIMAX or SARIMAX model to help capture the seasonal patterns.

This is implemented in NIXLA, and we can use the implementation shown here:

TBATS_model = TBATS(seasonal_periods  = 48, use_trend=True, use_damped_trend=True)

In NIXTLA, you can also use AutoTBATS to let the system optimize how to handle the various parameters.

Let’s see what the TBATS forecast looks like:

Figure 4.9: TBATS forecast

Again, the seasonality pattern has been replicated and is capturing most of the peaks in the forecast. Now let’s take a look at another method that is well suited for highly seasonal time series (even if it has multiple seasonalities like our case).

Multiple Seasonal-Trend decomposition using LOESS (MSTL)

Remember the time series decomposition we did back in Chapter 3? What if we can use the same techniques to forecast? That’s exactly what MSTL does. Let’s look at the components of a time series again:

Trend and cyclical components can be extracted using LOESS regression. If we fit a simple model on the trend values, we can use it to extrapolate to the future. And the seasonality component can easily be extrapolated because it is supposed to be a repeating pattern. Combining these, we get a forecasting model that works pretty well.

The MSTL method in NIXTLA applies the LOESS technique to decompose a time series into its various seasonal components. Following this decomposition, it employs a specialized non-seasonal model to forecast the trend, and a Seasonal Naive model to predict each of the seasonal components. This approach allows for the detailed analysis and forecasting of time series with complex seasonal patterns:

MSTL_model = MSTL(season_length  = 48)

Let’s see what the MSTL forecast looks like:

Figure 4.10: MSTL forecast

Let’s also take a look at how the different metrics that we chose did for each of these forecasts for the household we were experimenting with (from the notebook):

A screenshot of a computer Description automatically generated

Figure 4.11: Summary of all the baseline algorithms

Out of all the baseline algorithms we tried, AutoETS is performing the best on MAE as well as MSE. ARIMA was the second-best model followed by TBATS. However, if you look at the Time Elapsed column, TBATS stands out taking just 7.4 seconds vs. 19 seconds for ARIMA. Since they had similar performance, we will choose TBATS over ARIMA, along with AutoETS as our baseline, and run them on all 399 households in the dataset (both validation and test) we’ve chosen (the code for this is available in the 02-Baseline_Forecasts_using_NIXTLA.ipynb notebook).

Evaluating the baseline forecasts

Since we have the baseline forecasts generated from ETS as well as TBATS, we should also evaluate these forecasts. The aggregate metrics for all the selected households for both these methods are as follows:

A screenshot of a graph

Figure 4.12: The aggregate metrics of all the selected households (both validation and test)

It looks like AutoETS is performing much better in all three metrics. We also have these metrics calculated at a household level. Let’s look at the distribution of these metrics in the validation dataset for all the selected households:

Figure 4.13: The distribution of MASE and forecast bias of the baseline forecast in the validation dataset

The MASE histogram of ETS seems to have a smaller spread than TBATS. ETS also has a lower median MASE than TBATS. We can see a similar pattern for forecast bias as well, with the forecast bias of ETS centered around zero and much less spread.

Back in Chapter 1, Introducing Time Series, we saw why every time series is not equally predictable and saw three factors to help us think about the issue—understanding the Data Generating Process (DGP), the amount of data, and adequately repeating the pattern. In most cases, the first two are pretty easy to evaluate, but the third one requires some analysis. Although the performance of baseline methods gives us some idea about how predictable any time series is, they still are model-dependent. So, instead of measuring how well a time series is forecastable, we might be better measuring how well the chosen model can approximate the time series. This is where a few more fundamental techniques (relying on the statistical properties of a time series) come in.

Assessing the forecastability of a time series

Although there are many statistical measures that we can use to assess the predictability of a time series, we will just look at a few that are easier to understand and practical when dealing with large time series datasets. The associated notebook (02-Forecastability.ipynb) contains the code to follow along.

Coefficient of variation

The Coefficient of Variation (CoV) relies on the fact that the more variability that you find in a time series, the harder it is to predict it. And how do we measure variability in a random variable? Standard deviation.

In many real-world time series, the variation we see in the time series is dependent on the scale of the time series. Let’s imagine that there are two retail products, A and B. A has a mean monthly sale of 15, while B has 50. If we look at a few real-world examples like this, we will see that if A and B have the same standard deviation, B, which has a higher mean, is much more forecastable than A. To accommodate this phenomenon and to make sure we bring all the time series in a dataset to a common scale, we can use the CoV:

$$CoV_n = \frac{\sigma_n}{\mu_n}$$

Here, $\sigma_n$is the standard deviation, and $\mu_n$is the mean of the time series, n.

The CoV is the relative dispersion of data points around the mean, which is much better than looking at the pure standard deviation.

The larger the value for the CoV, the worse the predictability of the time series. There is no hard cutoff, but a value of 0.49 is considered a rule of thumb to separate time series that are relatively easier to forecast from the hard ones. Depending on the general hardness of the dataset, we can tweak this cutoff. Something I have found useful is to plot a histogram of CoV values in a dataset and derive cutoffs based on that.

Even though the CoV is widely used in the industry, it suffers from a few key issues:

To overcome these shortcomings, we propose another derived measure.

Residual variability

The thought behind residual variability (RV) is to try and measure the same kind of variability that we were trying to capture with the CoV but without the shortcomings. I was brainstorming on ways to avoid the problems of using the CoV, typically the seasonality issue, and was applying the CoV to the residuals after seasonal decomposition. It was then I realized that the residuals would have a few negative values and that the CoV wouldn’t work well. Stefan de Kok, who is a thought leader in demand forecasting and probabilistic forecasting, suggested using the mean of the original actuals, which worked.

To calculate RV, you must perform the following steps:

  1. Perform seasonal decomposition.
  2. Calculate the standard deviation of the residuals or the irregular component.
  3. Divide the standard deviation by the mean of the original observed values (before decomposition).

Mathematically, it can be represented as:

$$RV_n = \frac{\sigma_n^{res}}{\mu_n^\gamma}$$

where, $\sigma_n^{res}$ is the standard deviation of the residuals after decomposition and $\mu_n^\gamma$ is the mean of the original observed values.

The key assumption here is that seasonality and trend are components that can be predicted. Therefore, our assessment of the predictability of a time series should only look at the variability of the residuals. However, we cannot use CoV on the residuals because the residuals can have negative and positive values, so the mean of the residuals loses the interpretation of the level of the series and tends to zero. When residuals tend to zero, the CoV measure tends to infinity because of the division by mean. Therefore, we use the mean of the original series as the scaling factor.

Let’s see how we can calculate RV for all the time series in our dataset (which are in a compact form):

block_df["rv"] = block_df.progress_apply(lambda x: calc_norm_sd(x['residuals'],x['energy_consumption']), axis=1)

In this section, we looked at two measures that are based on the standard deviation of the time series. Now, let’s look at assessing the forecastability of a time series.

Entropy-based measures

Entropy is a ubiquitous term in science. We see it popping up in physics, quantum mechanics, social sciences, and information theory. And everywhere, it is used to talk about a measure of chaos or lack of predictability in a system. The entropy we are most interested in now is the one from information theory. Information theory involves quantifying, storing, and communicating digital information.

Claude E. Shannon presented the qualitative and quantitative model of communication as a statistical process in his seminal paper A Mathematical Theory of Communication. While the paper introduced a lot of ideas, some of the concepts that are relevant to us are information entropy and the concept of a bit—a fundamental unit of measurement of information.

Reference check:

A Mathematical Theory of Communication by Claude E. Shannon is cited as reference 3 in the References section.

The theory in itself is quite a lot to cover, but to summarize the key bits of information, take a look at the following short glossary:

Now, if we were to transfer one bit of information, it would reduce the uncertainty of the receiver by two. To understand this better, let’s consider a coin toss. We toss the coin in the air, and as it is spinning through the air, we don’t know whether it is going to be heads or tails. But we do know it is going to be one of these two. When the coin hits the ground and finally comes to rest, we find that it is heads. We can represent whether the coin toss is heads or tails with one bit of information (0 for heads and 1 for tails). So, the information that was passed to us when the coin fell reduced the possible outcomes from two to one (heads). This transfer was possible with one bit of information.

In information theory, the entropy of a discrete random variable is the average level of information, surprise, or uncertainty inherent in the variable’s possible outcomes. In more technical parlance, it is the expected number of bits required for the best possible encoding scheme of the information present in the random variable.

Additional reading:

If you want to intuitively understand entropy, cross-entropy, Kullback-Leibler divergence, and so on, head over to the Further reading section. There are a couple of links to blogs (one of which is my own) where we try to lay down the intuition behind these metrics.

Entropy is formally defined as follows:

$$H(X) = -\sum_{i=1}^{n} P(x_i) \cdot \log P(x_i)$$

Here, X is the discrete random variable with possible outcomes, x1, x2, …, xn. Each of those outcomes has a probability of occurring, which is denoted by P(x1), P(x2), …, P(xn).

To develop some intuition around this, we can think that the more spread out a probability distribution is, the more chaos is in the distribution, and thus more entropy. Let’s quickly check this with some code:

# Creating an array with a well balanced probability distribution
flat = np.array([0.1,0.2, 0.3,0.2, 0.2])
# Calculating Entropy
print((-np.log2(flat)* flat).sum())

>> 2.2464393446710154

# Creating an array with a peak in probability
sharp = np.array([0.1,0.6, 0.1,0.1, 0.1])
# Calculating Entropy
print((-np.log2(sharp)* sharp).sum())

>> 1.7709505944546688

Here, we can see that the probability distribution that spreads its mass has higher entropy.

In the context of a time series, n is the total number of time series observations, and P(xi) is the probability for each symbol of the time series alphabet. A sharp distribution means that the time series values are concentrated on a small area and should be easier to predict. On the other hand, a wide or flat distribution means that the time series value can be equally likely across a wider range of values and hence is difficult to predict.

If we have two time series—one containing the result of a coin toss and the other containing the result of a dice throw—the dice throw would have any output between one and six, whereas the coin toss would be either zero or one. The coin toss time series would have lower entropy and be easier to predict than the dice throw time series.

However, since time series is typically continuous, and entropy requires a discrete random variable, we can resort to a few strategies to convert the continuous time series into a discrete one. Many strategies, such as quantization or binning, can be applied, which leads to a myriad of complexity measures. Let’s review one such measure that is useful and practical.

Spectral entropy

To calculate the entropy of a time series, we need to discretize the time series. One way to do that is by using Fast Fourier Transform (FFT) and power spectral density (PSD). This discretization of the continuous time series is used to calculate spectral entropy.

We learned what Fourier Transform is earlier in this chapter and used it to generate a baseline forecast. But using FFT, we can also estimate a quantity called power spectral density. This answers the question, How much of the signal is at a particular frequency? There are many ways of estimating power spectral density from a time series, but one of the easiest ways is by using the Welch method, which is a non-parametric method based on Discrete Fourier Transform. This is also implemented as a handy function with the periodogram(x) signature in scipy.

The returned PSD will have a length equal to the number of frequencies estimated, but these are densities and not well-defined probabilities. So, we need to normalize PSD to be between zero and one:

$$nPSD_i = \frac{PSD_i}{\sum_{j=1}^{F} PSD_j}$$

Here, F is the number of frequencies that are part of the returned power spectrum density.

Now that we have the probabilities, we can just plug this into the entropy formula and arrive at the spectral entropy:

$$H_{s(X)} = -\sum_{i=1}^{n} nPSD_i \cdot \log(nPSD_i)$$

When we introduced entropy-based measures, we saw that the more spread out the probability mass of a distribution is, the higher the entropy is. In this context, the more frequencies across which the spectral density is spread, the higher the spectral entropy. So, a higher spectral entropy means the time series is more complex and, therefore, more difficult to forecast.

Since FFT has an assumption of stationarity, it is recommended that we make the series stationary before using spectral entropy as a metric. We can even apply this metric to a detrended and deseasonalized time series, which we can refer to as residual spectral entropy. This book’s GitHub repository contains an implementation of spectral entropy under src.forecastability.entropy.spectral_entropy. This implementation also has a parameter, transform_stationary, which, if set to True, will detrend the series before we apply spectral entropy. Let’s see how we can calculate spectral entropy for our dataset:

from src.forecastability.entropy import spectral_entropy
block_df["spectral_entropy"] = block_df.energy_consumption.progress_apply(lambda x: spectral_entropy(x, transform_stationary=True))
block_df["residual_spectral_entropy"] = block_df.residuals.progress_apply(spectral_entropy)

There are other entropy-based measures such as approximate entropy and sample entropy, but we will not cover them in this book. They are more computationally intensive and don’t tend to work for time series that contain fewer than 200 values. If you are interested in learning more about these measures, head over to the Further reading section.

Another metric that takes a slightly different path is the Kaboudan metric.

Kaboudan metric

In 1999, Kaboudan defined a metric for time series predictability, calling it the $\eta$-metric. The idea behind it is very simple. If we block-shuffle a time series, we are essentially destroying the information in the time series. Block shuffling is the process of dividing the time series into blocks and then shuffling those blocks. So, if we calculate the sum of squared errors (SSE) of a forecast that’s been trained on a time series and then contrast it with the SSE of a forecast trained on a shuffled time series, we can infer the predictability of the time series. The formula to calculate this is as follows:

$$\eta = 1 - \frac{SSE_Y}{SSE_S}$$

Here, SSEY is the SSE of the forecast that was generated from the original time series, while SSES is the SSE of the forecast that was generated from the block-shuffled series.

If the time series contains some predictable signals, SSEY would be lower than SSES and $\eta$ would approach one. This is because there was some information or patterns that were broken due to the block shuffling. On the other hand, if a series is just white noise (which is unpredictable by definition), there would be hardly any difference between SSEY and SSES, and $\eta$ would approach zero.

In 2002, Duan investigated this metric and suggested some modifications in his thesis. One of the problems he identified, especially in long time series, is that the $\eta$ values are found in a narrow band around 1 and suggested a slight modification to the formula. We call this the modified Kaboudan metric. The measure on the lower side is also clipped to zero. Sometimes, the metric can go below zero because SSES is lower than SSEY, which is because the series is unpredictable and, by pure chance, block shuffling made the SSE lower:

$$\eta_{modified} = 1 - \sqrt{\frac{SSE_Y}{SSE_S}}$$

Reference check:

The research paper that proposed the Kaboudan metric is cited as reference 4 in the References section. The subsequent modification that Duan suggested is cited as reference 5.

This modified version, as well as the original, has been implemented in this book’s GitHub repository.

There is no restriction on the forecasting model you use to generate the forecast, which makes it a bit more flexible. Ideally, we can choose one of the classical statistical methods that is fast enough to be applied to the whole dataset. But this also makes the Kaboudan metric dependent on the model, and the limitations of the model are inherent in the metric. The metric measures a combination of how difficult a series is to forecast and how difficult it is for the model to forecast the series.

Again, both metrics are implemented in this book’s GitHub repository. Let’s see how we can use them:

from src.forecastability.kaboudan import kaboudan_metric, modified_kaboudan_metric
model = Theta(theta=3, seasonality_period=48*7, season_mode=SeasonalityMode.ADDITIVE)
block_df["kaboudan_metric"] = [kaboudan_metric(r[0], model=model, block_size=5, backtesting_start=0.5, n_folds=1) for r in tqdm(zip(*block_df[["energy_consumption"]].to_dict("list").values()), total=len(block_df))]
block_df["modified_kaboudan_metric"] = [modified_kaboudan_metric(r[0], model=model, block_size=5, backtesting_start=0.5, n_folds=1) for r in tqdm(zip(*block_df[["energy_consumption"]].to_dict("list").values()), total=len(block_df))]

Although there are many more metrics we can use for this purpose, the metrics we just reviewed for assessing forecastability cover a lot of the popular use cases and should be more than enough to gauge any time series dataset in regards to the difficulty of forecasting it. We can use these metrics to compare one time series with another time series or to profile a whole set of related time series in a dataset with another dataset for benchmarking purposes.

Additional reading:

If you want to delve a little deeper and analyze the behavior of these metrics, how similar they are to each other, and how effective they are in measuring forecastability, go to the end of the 03-Forecastability.ipynb notebook. We compute rank correlations among these metrics to understand how similar these metrics are. We can also find rank correlations with the computed metrics from the best-performing baseline method to understand how well these metrics did in estimating the forecastability of a time series. I strongly encourage you to play around with the notebook and understand the differences between the different metrics. Pick a few time series and check how the different metrics give you slightly different interpretations.

Congratulations on generating your baseline forecasts—the first set of forecasts we have generated using this book! Feel free to head over to the notebooks, play around with the parameters of the methods, and see how forecasts change. It’ll help you develop an intuition around what the baseline methods are doing. If you are interested in learning more about how to make these baseline methods better, head over to the Further reading section, where we have provided a link to the paper The Wisdom of the Data: Getting the Most Out of Univariate Time Series Forecasting, by F. Petropoulos and E. Spiliotis.

Summary

And with this, we have come to the end of Part 1, Getting Familiar with Time Series. We have come a long way from just understanding what a time series is to generating competitive baseline forecasts. Along the way, we learned how to handle missing values and outliers and how to manipulate time series data using pandas. We used all those skills on a real-world dataset regarding energy consumption. We also looked at ways to visualize and decompose time series. In this chapter, we set up a test harness, learned how to use the NIXTLA library to generate a baseline forecast, and looked at a few metrics that can be used to understand the forecastability of a time series.

For some of you, this may be a refresher, and we hope this chapter added some value in terms of some subtleties and practical considerations. For the rest of you, we hope you are in a good place foundationally to start venturing into modern techniques using machine learning in the next part of the book.

In the next chapter, we will discuss the basics of machine learning and delve into time series forecasting.

References

The following references were provided in this chapter:

  1. Assimakopoulos, Vassilis and Nikolopoulos, K. (2000). The theta model: A decomposition approach to forecasting. International Journal of Forecasting. 16. 521-530. https://www.researchgate.net/publication/223049702_The_theta_model_A_decomposition_approach_to_forecasting.
  2. Rob J. Hyndman, Baki Billah. (2003). Unmasking the Theta method. International Journal of Forecasting. 19. 287-290. https://robjhyndman.com/papers/Theta.pdf.
  3. Shannon, C.E. (1948), A Mathematical Theory of Communication. Bell System Technical Journal, 27: 379-423. https://people.math.harvard.edu/~ctm/home/text/others/shannon/entropy/entropy.pdf.
  4. Kaboudan, M. (1999). A measure of time series’ predictability using genetic programming applied to stock returns. Journal of Forecasting, 18, 345-357: http://www.aiecon.org/conference/efmaci2004/pdf/GP_Basics_paper.pdf.
  5. Duan, M. (2002). TIME SERIES PREDICTABILITY: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.68.1898&rep=rep1&type=pdf.
  6. De Livera, A. M., & Hyndman, R. J. (2009). Forecasting time series with complex seasonal patterns using exponential smoothing (Department of Econometrics and Business Statistics Working Paper Series 15/09)
  7. Hyndman, Rob. “Rob J Hyndman - TBATS with Regressors.” Rob J Hyndman, 6 Oct. 2014, http://robjhyndman.com/hyndsight/tbats-with-regressors

Further reading

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