7

Target Transformations for Time Series Forecasting

In the previous chapter, we delved into how we can do temporal embedding and time delay embedding by making use of feature engineering techniques. But that was just one side of the regression equation—the features. Often, we see that the other side of the equation—the target—does not behave the way we want. In other words, the target doesn’t have some desirable properties that make forecasting easier. One of the major culprits in this area is stationarity—or more specifically, the lack of it. And it creates problems with the assumptions we make while developing a machine learning (ML)/statistical model. In this chapter, we will look at some techniques for handling such problems with the target.

In this chapter, we will cover the following topics:

Technical requirements

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

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

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

Detecting non-stationarity in time series

Stationarity is a prevalent assumption in most econometrics models and is a rigorous and mathematical concept. But without getting into a lot of math, we can intuitively think about stationarity as the state where the statistical properties of the distribution from which the time series is sampled remain constant over time. This is relevant in time series as regression as well because we are estimating a single forecasting function across time. And if the behavior of the time series changes with time, the single function that we estimate may not be relevant all the time. For instance, if we think about the number of visitors to the nearby park in a day as a time series, we know that those patterns are going to be very different for pre- and post-pandemic periods. In the ML world, this phenomenon is called concept drift.

Intuitively, we can understand that it is easier to forecast a stationary series than a non-stationary series. But here comes the punchline: in the real world, almost all time series do not satisfy the stationarity assumption—more specifically, the strict stationarity assumption. Strict stationarity is when all the statistical properties such as the mean, variance, skewness, and so on do not change with time. Many times, this strict stationarity assumption is relaxed in favor of weak stationarity, where we only stipulate that the mean and the variance of the time series do not change with time.

There are four main questions we can ask ourselves to check whether our time series is stationary or not:

Out of these questions, the first three can be ascertained using a simple visual inspection. Unit roots are more difficult to understand. We will take a deeper look at unit roots shortly. Let’s take a look at a few time series and check whether we can tell whether they are stationary or not via visual inspection (you can note your answers):

Figure 7.1 – Testing your understanding of stationarity

Figure 7.1: Testing your understanding of stationarity

Now, check your responses and see how many of them you guessed correctly. If you got at least four out of six, you are doing great with your intuition of stationarity:

When we have hundreds, or even millions, of time series, we can’t practically do a visual inspection to ascertain whether they are stationary or not. So, now, let’s look at a few ways of detecting these key properties using statistical tests and also how to try and correct them.

Although we are talking about correcting or making a time series stationary, it is not always essential to do that in the ML paradigm because some of these can be handled by using the right kind of features in the model. Whether to make a series stationary or not is a decision we will have to make after experimenting with the techniques. This is because, as you will see, while there are advantages to making a series stationary, there are also disadvantages to using some of these techniques, as we will see when we discuss each transformation in detail.

Notebook alert:

To follow along with the complete code, use the 02-Dealing_with_Non-Stationarity.ipynb notebook in the Chapter06 folder.

Detecting and correcting for unit roots

Let’s talk about unit roots first since this is what is most commonly tested for stationarity. Time series analysis has its roots in econometrics and statistics and unit root is a concept derived directly from those fields.

Unit roots

Unit roots are quite complicated to understand fully but to develop some intuition, we can look at a simplification. Let’s consider an autoregressive model of order 1(AR(1) model):

$y_t = \phi y_{t-1} + \epsilon_t$, where $\epsilon_t$ is the white noise and $\phi$ is the AR coefficient.

If we think about the different values of $\phi$ in the equation, we can come up with three scenarios (Figure 7.2):

  1. $|\phi| \gt 1$: When $|\phi|$ is greater than 1, every successive value in the time series is multiplied by a number greater than 1, which means it will have a strong and rapidly increasing/decreasing trend and thereby be non-stationary.
  2. $|\phi| \lt 1$: When $|\phi|$ is less than 1, every successive value in the time series is multiplied by a number less than 1, which means over the long term, the mean of the series trends to zero and will oscillate around it. Therefore, it is stationary.
  3. $|\phi| = 1$: When $|\phi|$ is equal to 1, things become trickier. When $|\phi| = 1$ for an AR(1) model, this is known as it having a unit root and the equation becomes $y_t = y_{t-1} + \epsilon_t$. This is called random walk in econometrics and is a very popular kind of time series in financial and economic domains. Mathematically, we can prove that such a series will have a constant mean but a non-constant variance:
Figure 7.2 – Autoregressive time series with different parameters

Figure 7.2: Autoregressive time series with different $\phi$ parameters. Top: Scenario 1, phi <1; Middle: Scenario 2, phi = 1; Bottom: Scenario 3, phi>1

While we discussed unit roots in an AR(1) process, we can extend the same intuition to multiple lags or an AR(p) model. Calculating and testing unit roots is more complicated there, but still possible.

So, now that we know what a unit root is, how can we statistically test this? This is where the Dickey-Fuller test comes in.

The Augmented Dickey-Fuller (ADF) test

The null hypothesis in this test is that the $\phi$ in an AR(1) model of the time series is equal to 1, and by extension non-stationary. The alternate hypothesis is that the $\phi$ in the AR(1) model is less than 1. The ADF test takes the Dickey-Fuller test and extends it to an AR(p) model because most time series are not defined by just one lag of the time series. This is the standard and most popular statistical test to check for unit roots. The core of the test involves running a regression on the lags of the time series and calculating the statistic on the variance of the residuals.

Let’s see how we can do this in Python using statsmodels:

from statsmodels.tsa.stattools import adfuller
result = adfuller(y)

result from adfuller is a tuple that contains the test statistic, p-value, and critical values at different confidence levels. Here, we are most interested in the p-value, which is an easy and practical way to check whether the null hypothesis is rejected or not. If p < 0.05, there is a 95% probability that the series does not have a unit root; the series is stationary from a unit root perspective.

To make this process even easier, we have included a method called check_unit_root in src.transforms.stationary_utils that does the inference for you (comparing the returned probability with the confidence and rejecting or accepting the null hypothesis) and returns a namedtuple with a Boolean attribute called stationary, and the entire results from statsmodels in results:

from src.transforms.stationary_utils import check_unit_root
# We pass the time series along with the confidence with which we need the results
check_unit_root(y, confidence=0.05)

Now that we’ve learned how to check whether a series has a unit root or not, how do we make it stationary? Let’s look at a few transforms that help us do that.

Differencing transform

The differencing transform is a very popular transform to make a time series stationary, or at least get rid of unit roots. The concept is simple: we transform the time series from the domain of observation to the domain of change in observations. The differencing transform subtracts subsequent observations from one another:

zt = yt – yt-1

Differencing helps us stabilize the mean of the time series and, with that, reduce or eliminate trend and seasonality. Let’s see how differencing can make a series stationary.

Let the time series in question be $y_t = \beta_0 + \beta_1 t + \epsilon_t$ , where $\beta_0$ and $\beta_1$ are the coefficients and $\epsilon$ is white noise. From this equation, we can see that time, t, is part of the equation, making yt a time series with a trend. So, the differenced time series z would be as follows:

$$z_t = y_t = y_{t-1}$$

$$= (\beta_0 + \beta_1 t + \epsilon_t) - (\beta_0 + \beta_1 (t - 1) + \epsilon_{t-1})$$

$$= \beta_1 + (\epsilon_t - \epsilon_{t-1})$$

What we need to look for in this new equation is that there is no mention of t. This means that the dependence on t, which created the trend, has been removed, and now the time series has constant mean and variance at any point in time.

Differencing does not remove all kinds of non-stationarity but works for the majority of time series. But there are a few drawbacks to this approach as well. One of them is that we lose the scale of the time series while modeling. Many times, the scale of the time series holds some information that is useful for forecasting. For instance, in a supply chain, SKUs with higher sales exhibit a different kind of pattern from SKUs with lower sales and when we do differencing, this information about the distinction is lost.

Another drawback is more from an operational point of view. When we use differencing for forecasting, we also need to inverse the transform after we get the differenced output from the model. This is an additional layer of complexity that we have to manage. One way is to keep the most recent observation in memory and keep adding the differences to it to inverse the transform. Another way is to have yt-1 ready for every t that we need to inverse transform and keep adding the difference to yt-1.

We have implemented the latter using the datetime index as a key to align and fetch the yt-1 observation in src.transforms.target_transformations.py in this book’s GitHub repository. Let’s see how we can use it:

from src.transforms.target_transformations import AdditiveDifferencingTransformer
diff_transformer = AdditiveDifferencingTransformer()
# [1:] because differencing reduces the length of the time series by one
y_diff = diff_transformer.fit_transform(y, freq="1D")[1:]

y_diff will have the transformed series. To get back to the original time series, we can call inverse_transform using diff_transformer. The associated notebook has examples and plots to see how the differencing changes the time series.

Here, we saw differencing as the process of subtracting subsequent values in the time series. But we can also do differencing with other operators such as division (yt/ yt-1), which is implemented in the src.transforms.target_transformations.py file as MultiplicativeDifferencingTransformer. We can also experiment with these transforms to check whether these work best for your dataset.

Although differencing solves the majority of stationarity issues, it’s not guaranteed to take care of all kinds of trends (non-linear or piecewise trends), seasonality, and so on. Sometimes, we may not want to difference the series but still handle trends and seasonality. So, let’s see how we can detect and remove trends in a time series.

Detecting and correcting for trends

In Chapter 5, Time Series Forecasting as Regression, we talked about forecasting being a difficult problem because it is intrinsically an extrapolation problem. Trends are one of the major contributors to forecasting being an extrapolation problem. If we have a time series that is trending upward, any model that attempts to forecast it needs to extrapolate beyond the range of values it has seen during training. ARIMA handles this using autoregression, whereas exponential smoothing handles it by modeling the trend explicitly. But standard regression may not be naturally suited to extrapolation. However, with suitable features, such as lags, it can start to do that.

But if we can confidently estimate and extract a trend in the time series, we can simplify the problem we have to apply regression to by detrending the time series.

But before we move ahead, it is worth learning about two major types of trends.

Deterministic and stochastic trends

Let’s take the simple AR(1) model we saw earlier to develop intuitions about this one too. Earlier, we saw that having $\phi \gt 1$ in an AR(1) model leads to a trend in the time series. But another way we can think about a trending time series is if we include time as an ordinal variable in the equation defining the time series. For instance, let’s consider two time series:

Time series 1:

$$y_t = \phi y_{t-1} + \epsilon_t; \epsilon_t \sim \mathcal{N}(0, \sigma^2)$$

Time series 2:

$$y_t = \beta_0 + \beta_1 t + \epsilon_t; \epsilon_t \sim \mathcal{N}(0, \sigma^2)$$

These can be seen in the following graphs:

Figure 7.3 – Time series 1 (stochastic trend) and Time series 2 (deterministic trend)

Figure 7.3: Top: Stochastic trend; Bottom: Deterministic trend

We saw both of these equations earlier; Time series 1 is the AR(1) model, while Time series 2 is the time series equation we chose to illustrate differencing. We already know that for $\phi \gt 1$, both Time series 1 and Time series 2 have trends. But there is a difference between the two trends.

In Time series 2, the trend is constant and can be perfectly modeled. In this case, just a linear fit would explain the trend perfectly. But in Time series 1, the trend is not something that can be explained by a simple linear fit. It is inherently dependent on the previous value of the time series that has $\epsilon_{t-1}$ and hence is stochastic. Therefore, Time series 2 has a deterministic trend and Time series 1 has a stochastic trend.

We can use the same ADF test we saw earlier in this chapter to check whether a time series has deterministic or stochastic trends. Without going into the math of the statistical test, we know that it tests for a unit root by fitting an AR(p) model to the time series. There are a few variants of this test that we can specify using the regression parameter in the statsmodels implementation. This parameter takes in the following values:

So, if we run an ADF test with regression="c", it will be non-stationary. Now, if we run the ADF test with regression="ct", it will come out as stationary. This means that when we removed a deterministic trend from the time series, it became stationary. This test is what we can use to determine whether a trend that we observe in a time series is deterministic or stochastic. In the Further reading section, we have provided a link to a blog post by Fabian Kostadinov, where he experiments with a few time series to make the distinction between the different variants of ADF tests clear.

We have implemented this test in src.transforms.stationary_utils as check_deterministic_trend, which does the inference for you and returns a namedtuple with a Boolean attribute of deterministic_trend. The namedtuple also includes the raw results from the two adfuller tests we did under adf_res and adf_ct_res if you want to investigate further. Let’s see how we can use this test:

check_deterministic_trend(y, confidence=0.05)

This will tell us whether the trend is stationary or deterministic. Now, let’s look at a couple of ways to identify and statistically test trends (irrespective of whether it is deterministic or not) in a time series.

Kendall’s Tau

Kendall’s Tau is a measure of correlation but is carried out on the ranks of the data. Kendall’s Tau is a non-parametric test and therefore does not make assumptions about the data. The correlation coefficient, Tau, returns a value between -1 and 1, where 0 shows no relationship and 1 or -1 is a perfect relationship. We will not dive into the details of how Kendall’s Tau is calculated or how the significance test is done as this is outside the scope of this book. The Further reading section contains a link that explains this well.

In this section, we will see how we can use Kendall’s Tau to measure the trend in our time series. As mentioned earlier, Kendall’s Tau calculates a rank correlation between two variables. If we chose one of those variables as the time series and set the other as the ordinal representation of time, the resulting Kendall’s Tau would represent the trend in the time series. An additional benefit is that the higher the value of Kendall’s Tau is, the stronger we expect the trend to be.

scipy has an implementation of Kendall’s Tau that we can use as follows:

import scipy.stats as stats
tau, p_value = stats.kendalltau(y, np.arange(len(y)))

We can compare the returned p-value to our required confidence (typically, this is 0.05) and say that if p_value < confidence, we conclude that the trend is statistically significant. The sign of tau tells us whether this is an increasing trend or a decreasing one.

We have made an implementation of Kendall’s Tau in src.transforms.stationary_utils as check_trend, which checks for the presence of a trend for you. The only parameters we need to provide are as follows:

A few more parameters are there, but those are for the Mann-Kendall (M-K) test, which will be explained next.

Let’s see how we can use this test:

check_trend(y, confidence=0.05)

This method also checks whether the trend that has been identified is deterministic or stochastic and calculates the direction of the trend. The result is returned as a namedtuple with the following parameters:

Now, let’s look at the Mann-Kendall test.

Mann-Kendall test (M-K test)

The Mann-Kendall test is used to check for the presence of a monotonic upward or downward trend. And since the M-K test is a non-parametric test, like Kendall’s Tau, there is no assumption of normality or linearity. The test is done by analyzing the signs between consecutive points in the time series. The crux of the test is the idea that in the presence of a trend, the sign values, if summed up, increase or decrease constantly.

Although non-parametric, there were a few assumptions in the original test:

Numerous alterations have been made to the original tests to tackle these problems over the years and a lot of such alterations, along with the original test, have been implemented at https://github.com/mmhs013/pyMannKendall. They are available in pypi as pymannkendall.

Pre-whitening is a common technique used to remove the autocorrelation in a time series. In a nutshell, the idea is as follows:

  1. Identify $\phi$ with an AR(1) model
  2. $y_t^{prewhiten} = y_t - \phi * y_{t-1}$

M. Bayazit and B. Önöz (2007) suggested not using pre-whitening before doing the M-K test if the sample size is larger than 50 and if the trend is strong enough (slope>0.01). For seasonal data, a seasonal variant of the M-K test has also been implemented in pymannkendall.

Reference check:

The research paper by M. Bayazit and B. Önöz is cited in the References section under reference 1.

The same method we discussed earlier, check_trend, also implements M-K tests that can be enabled by setting mann_kendall=True. However, one thing we need to keep in mind is that the M-K test is considerably slower than Kendall’s Tau, especially for long time series. There are a couple more parameters specific to the M-K test:

Let’s see how we can use this test:

check_trend(y, confidence=0.05, mann_kendall=True)

The result is returned as a namedtuple with the following parameters:

Let’s see an example where we applied both these tests on a time series (see 02-Dealing_with_Non-Stationarity.ipynb for the full code):

# y_unit_root is the a synthetic unit root timeseries
y_unit_root.plot()
plt.show()

Figure 7.4: M-K test

kendall_tau_res = check_trend(y_unit_root, confidence=0.05)
mann_kendall_res = check_trend(y_unit_root, confidence=0.05, mann_kendall=True)
print(f"Kendalls Tau: Trend: {kendall_tau_res.trend} | Direction: {kendall_tau_res.direction} | Deterministic: {kendall_tau_res.deterministic}")
print(f"Mann-Kendalls: Trend: {mann_kendall_res.trend} | Direction: {mann_kendall_res.direction} | Deterministic: {mann_kendall_res.deterministic}")
## Output

>> Kendalls Tau: Trend: True | Direction: decreasing | Deterministic: False
>> Mann-Kendalls: Trend: True | Direction: decreasing | Deterministic: False

It would benefit you if you can generate some time series, or even pick a few time series you have come across and use these functions to see how it works and how the results help you. The associated notebook has some examples to get you started. You can observe how the direction and slope are different for different types of trends.

Now that we know how to detect a trend, let’s look at detrending.

Detrending transform

If the trend is deterministic, removing the trend would add some value to the modeling procedure. In Chapter 3, Analyzing and Visualizing Time Series Data, we discussed detrending as it was an integral part of the decomposition we were doing. But techniques such as moving average or LOESS regression have one drawback—they can’t extrapolate. But if we are considering a deterministic linear (or even polynomial) trend, it can be easily estimated by using linear regression. The added advantage here is that the trend that is identified can easily be extrapolated.

The procedure is simple: we regress the time series on the ordinal representation of time and extract the parameters. Once we have these parameters, using the dates, we can extrapolate the trend to any point in the future. The core logic in Python is shown below:

# y is the time series we are detrending
x = np.arange(len(y))
# degree is the degree of trend we are estimating. Linear, or polynomial
# Fitting a regression on y using a linearly increasing x
linear_params = np.polyfit(x=x, y=y, deg=degree)
# Extract trend using fitted parameters
trend = get_trend(y)
# Now this extracted trend can be removed from y
detrended = y - trend

We have made and implemented this detrender as a transformer in src.transforms.target_transformations.py as DetrendingTransformer. You can see how we have implemented it in the GitHub repository. Now, let’s see how we can use it:

from src.transforms.target_transformations import DetrendingTransformer
detrending_transformer = DetrendingTransformer(degree=1)
y_detrended = detrending_transformer.fit_transform(y, freq="1D")

y_detrended will contain the detrended series. To get the original time series back, we can call inverse_transform using detrending_transformer. The associated notebook has examples and plots to see how the detrending changes the time series.

Best practice:

We have to be careful with the trend assumptions, especially if we are forecasting for the long term. Even a linear trend assumption can lead to an unrealistic forecast because trends don’t continue the same way forever in the real world. It is always advisable to dampen the trend by some factor, $\phi$, to be conservative in our extrapolation of the trend. This dampening can be as simple as $f_{t+h}^{damped} = f_{t+h} \times \phi^h, \text{ where } \phi \lt 1$.

Another key aspect that makes a time series non-stationary is seasonality. Let’s look at how to identify seasonality and remove it.

Detecting and correcting for seasonality

A vast majority of real-world time series have seasonality such as retail sales, energy consumption, and so on. And generally, the presence or absence of seasonality comes as part of the domain knowledge. But when we are working with a time series dataset, the domain knowledge becomes slightly diluted. The majority of time series may exhibit seasonality, but that doesn’t mean every time series in the dataset is seasonal. For instance, within a retail dataset, there might be items that are seasonal and some items that are not. Therefore, when working with a time series dataset, being able to determine whether a particular time series is seasonal or not has some value.

Detecting seasonality

There are two popular ways to check for seasonality, apart from just eyeballing it: autocorrelation and fast Fourier transform. Either is equally capable of identifying the seasonality period automatically. For our discussion, we’ll cover the autocorrelation method and examine how we can use that to determine seasonality.

Autocorrelation, as explained in Chapter 3, Analyzing and Visualizing Time Series Data, is the correlation of a time series to its lagged values. Typically, we expect the correlation to be higher in the immediate lags (lag 1, lag 2, and so on) and gradually die down as we move farther into the past. But for time series with seasonality, we will also see a spike in the seasonal periods.

Let’s understand this by looking at an example. Consider a synthetic time series that is just white noise combined with a sinusoidal signal with a seasonality cycle of 25 (identical to the seasonal time series we saw earlier, in Figure 7.1):

#WhiteNoise + Seasonal
y_random = pd.Series(np.random.randn(length), index=index)
t = np.arange(len(y_random))
y_seasonal = (y_random+1.9*np.cos((2*np.pi*t)/(length/4)))

If we plot the autocorrelation function (ACF) for this time series, it will look as follows (the code to calculate and plot this can be found in the 02-Dealing_with_Non-Stationarity.ipynb notebook):

Figure 7.4 – Autocorrelation plot of the synthetic time series with a seasonality cycle of 25

Figure 7.5: Autocorrelation plot of the synthetic time series with a seasonality cycle of 25

We can see that apart from the first few lags, the autocorrelation increases as we approach the seasonal cycle and peaks at the exact seasonality. We can use this property of the ACF to detect seasonality. darts, a time series forecasting library, has an implementation of this technique that identifies seasonality. But since it was designed to work for the time series data structure of darts, we have adapted the same logic to work on regular pandas series in src.transforms.stationary_utils.py under the name check_seasonality. The implementation can do two kinds of seasonality checks. It can take a seasonality_period as input and verify whether a seasonality corresponding to that seasonality_period exists in the data or not. If we do not give a seasonality_period ahead of time, it will return to you the shortest seasonality_period that is statistically significant.

The procedure, at a high level, does the following:

  1. It calculates the ACF.
  2. It finds all the relative maxima in the ACF. A relative maximum is a point where the function changes direction from increasing to decreasing.
  3. It checks whether the provided seasonal_period is a relative maximum. If not, we conclude there is no seasonality associated with seasonality_period.
  4. Now, we take the assumption that the ACF is normally distributed and compute the upper limit at the specified confidence level. The upper bound is given by:

$$UB = z_{1-\frac{\alpha}{2}} \times SE(r_h)$$

where rh is the estimated autocorrelation at lag h, SE is the standard error, and $z_{1-\frac{\alpha}{2}}$ is the quantile of the normal distribution based on the required confidence, $\alpha$. The SE is approximated using Bartlett’s formula (for the math behind this, head over to the Further reading section).

  1. Each of our candidates for seasonality_period is checked against this upper limit and the ones that are above this limit are deemed statistically significant.

There are only three parameters for this function, apart from the time series itself:

Let’s see how we can use this function on the same data we saw in Figure 7.4 (with a seasonal period of 25). This will give you a namedtuple with seasonal, a Boolean flag to indicate seasonality, and seasonal_periods, the seasonal periods with significant seasonality, as parameters:

# Running the function without specifying seasonal period to identify the seasonality
seasonality_res = check_seasonality(y_seasonal, max_lag=60, confidence=0.05)
print(f"Seasonality identified for: {seasonality_res.seasonal_periods}")
## Output

>> Seasonality identified for: 25

This function can also be used to verify if your assumption about the seasonality is right.
# Running the function specifying seasonal period to verify
seasonality_res = check_seasonality(y_seasonal, max_lag=30, seasonal_period=25, confidence=0.05)print(f"Seasonality Test for 25th lag: {seasonality_res.seasonal}")
## Output

>> Seasonality Test for 25th lag: True

Now that we know how to identify and test for seasonality, let’s talk about deseasonalizing.

Deseasonalizing transform

In Chapter 3, Analyzing and Visualizing Time Series Data, we reviewed techniques for seasonal decomposition. We can use the same techniques here as well, but with just one tweak. Earlier, we were not concerned with projecting the seasonality into the future. But when we are using deseasonalizing in forecasting, it is essential to be able to project it into the future as well. We are in luck since projecting the seasonal cycle forward is trivial. This is because we are looking at a fixed seasonality profile that will always keep repeating in the seasonal cycle. For instance, if we identified a seasonality profile for the 12 months of a year (yearly seasonality at monthly frequency data), the seasonality that’s extracted for these 12 months will just repeat itself in chunks of 12 months.

Using this property, we have implemented a transformer in src.transforms.target_transformations.py as DeseasonalizingTransformer. There are a few parameters and properties that we need to be aware of:

Let’s see how we can use it:

from src.transforms.target_transformations import DeseasonalizingTransformer
deseasonalizing_transformer = DeseasonalizingTransformer(seasonality_extraction="period_averages",seasonal_period=25)
y_deseasonalized = deseasonalizing_transformer.fit_transform(y, freq="1D")

y_deseasonalized will have the deseasonalized time series. To get back to the original time series, we can use the inverse_transform function. Typically, this can be used to add the seasonality back after making predictions.

Best practice:

Modeling seasonality can be done either separately, as discussed here, or by using the seasonal features that we discussed earlier in this chapter. Although the final evaluation on which one works better has to be found out empirically for each dataset, we can have a few rules of thumb/guidelines to decide on priority.

When we have enough data, letting the model learn seasonality as part of the main forecasting problem seems to work better. But in cases where data is not that rich, extracting seasonality separately before feeding it to an ML model works well.

When the dataset has varied seasonality (different seasonal cycles for different time series), then it should be treated accordingly. Either deseasonalize each time series separately or split the global ML model into different local models, each with its own seasonality pattern.

The last aspect that we talked about earlier is heteroscedasticity. Let’s quickly take a look at that as well.

Detecting and correcting for heteroscedasticity

Despite having a scary name, heteroscedasticity is a simple enough concept. It is derived from ancient Greek, where hetero means different and skedasis means dispersion. True to its name, heteroscedasticity is defined when the variability of a variable is different across another variable. In the context of a time series, we say a time series is heteroscedastic when the variability or dispersion of the time series varies with time. For instance, let’s think about the spending of a household over a number of years. In these years, this particular household went from being poor to middle class and finally upper middle class. When the household was poor, the spending was less and only on essentials, and because of that, the variability in spending was less. But as they approached the upper middle class, the household could afford luxuries, which created spikes in the time series and therefore higher variability. If we refer back to Figure 7.1, we can see what a heteroscedastic time series looks like.

But in addition to visual inspection, it would be neat if we could carry out an automated statistical test to ascertain heteroscedasticity.

Detecting heteroscedasticity

There are many ways to detect heteroscedasticity, but we will be using one of the most popular techniques, known as the White test, proposed by Halbert White in 1980. The White test uses an auxiliary regression task to check for constant variance. We run an initial regression using some covariates and calculate the residuals of this regression. Then, we fit another regression model with these residuals as the target and the covariates used in the first regression, and their squares and cross-products. The final statistic is estimated by using the R2 value of this auxiliary regression. For a detailed account of the test, head over to the Further reading section; for the rigorous mathematical procedure, the research paper is cited in the References section.

Reference check:

To learn more about the rigorous mathematical procedure of the White test, take a look at the research paper cited in the References section under reference 2.

In the context of a time series, we adapt this formulation by using a deterministic trend model. The initial regression is done by using time as an ordinal variable and the residuals are used to carry out the White test. The White test has an implementation in statsmodels of het_white, which we will be using to carry out this test. The het_white test returns two statistics and p-values—Lagrangian Multiplier and F-Statistic. Lagrangian Multiplier tests if there is any relationship between the variance of the residuals and the independent variables in the regression model. F-Statistic compares the fit of your original model to a model allowing for varying error variance. A p-value less than confidence in either of these tests indicates heteroscedasticity. But to be conservative, we can also use both tests and mark something as heteroscedastic only when both of the p-values are less than confidence.

We have wrapped all of this in a helpful function in src.transforms.stationary_utils as check_heteroscedasticity, which has only one additional parameter—confidence. Let’s see the core implementation of the method in Python:

import statsmodels.api as sm
# Fitting a linear trend regression
x = np.arange(len(y))
x = sm.add_constant(x)
model = sm.OLS(y,x)
results = model.fit()
# Using the het_white test on residuals
lm_stat, lm_p_value, f_stat, f_p_value = het_white(results.resid, x)
# Checking if both p values are less than confidence
if lm_p_value<confidence and f_p_value < confidence:
        hetero = True
    else:
        hetero = False

Now, let’s see how we can use this function:

from src.transforms.stationary_utils import check_heteroscedastisticity
check_heteroscedastisticity(y, confidence=0.05)

This returns a namedtuple with the following parameters:

Best practice:

The heteroscedasticity test we are doing only considers a trend in the regression and therefore, in the presence of seasonality, may not work very well. It is advised to deseasonalize the data before applying the function.

Detecting heteroscedasticity was the easier part. There are a few transforms that attempt to remove heteroscedasticity but with advantages and disadvantages. Let’s take a look at a few such transforms.

Log transform

Log transform, as the name suggests, is about applying a logarithm to the time series. There are two main properties of a log transform—variance stabilization and reducing skewness—thereby making the data distribution more normal. And out of these, we are more interested in the first property because that is what combats heteroscedasticity.

Log transforms are typically known to reduce the variance of the data and thereby remove heteroscedasticity in the data. Intuitively, we can think of a log transform as something that pulls in the extreme values on the right of the histogram, at the same time stretching back the very low values on the left of the histogram.

But it has been shown that the log transform does not always stabilize the variance. In addition to that, the log transform poses another challenge in ML. The optimization of loss now happens on the log scale. Since the log transformation compresses the lower end of the value range more than the higher one, the learned model can be less sensitive to errors in the lower range as compared to the higher one. Another key disadvantage is that the log transform can only be applied to strictly positive data. And if any of your data is zero or less than zero, then you will need to offset the whole distribution by adding some constant, M, and then applying the transform. This will also create some disturbance in the data, which can have adverse effects.

The bottom line is that we should be careful when applying a log transform. We have implemented a transformer in src.transforms.target_transformations.py as LogTransformer with just one parameter, add_one, which adds one before the transform and subtracts one after the inverse. The key logic in Python is as simple as applying an np.log1p or np.log function in the transform and reversing it with np.expm1 or np.exp, respectively:

# Transform
np.log1p(y) if self.add_one else np.log(y)
# Inverse Transform
np.expm1(y) if self.add_one else np.exp(y)y_log = log_transformer.fit_transform(y)

All we have done is wrap this into a nice and easy-to-use transformer. Let’s see how we can use it:

from src.transforms.target_transformations import LogTransformer
log_transformer = LogTransformer(add_one=True)
y_log = log_transformer.fit_transform(y)

y_log is the log-transformed time series. We can call inverse_transform to get the original time series back.

Box-Cox transformation

The log transform, although effective and common, is very strong. But the log is not the only monotonic transform that we can use. There are many other transforms, such as y2, $\frac{1}{y}$, $\frac{1}{\sqrt{y}}$ and so on, which are collectively part of the family of power transforms. One set of transforms that is very famous and widely used in this family is the Box-Cox transformations:

$$y(\lambda) = \frac{y^\lambda - 1}{\lambda}, \quad \text{if } \lambda \neq 0$$

and , $\log(y), \quad \text{if } \lambda = 0$

Reference check:

The original research paper by Box and Cox is cited in the References section under reference 3.

Intuitively, we can see that the Box-Cox transformation is a generalized logarithm transform. The log transform is just a special case of the Box-Cox transformation (when $\lambda = 0$). At different values of $\lambda$, it approximates other transforms such as y2 when $\lambda = 2$, $\frac{1}{\sqrt{y}}$ when $\lambda = -0.5$, $\sqrt{y}$ when $\lambda = 0.5$, and so on. When $\lambda = 1$, there is no major transformation.

A lot of the disadvantages that we mentioned for log transforms apply here as well, but the degree to which those effects are there varies, and we have a parameter, $\lambda$, to help us decide on the right level of those effects. Like log transforms, Box-Cox transformations also only use strictly positive data. The same addition of a constant to offset the data distribution has to be done here as well. The flip side of the parameter is that there is one more hyperparameter to tune.

There are a few automated methods to find the optimum $\lambda$ for any data distribution. One of them is by minimizing the log-likelihood of the data distribution, assuming normality. So, essentially, what we will be doing is finding the optimal $\lambda$ that makes the data distribution most normal. This optimization is already implemented in popular implementations such as the boxcox function in the scipy.special module in scipy.

Another way to find the optimal $\lambda$ is to use Guerrero’s method, which is typically suited for a time series. In this method, instead of trying to conform the data distribution to a normal distribution, we try to minimize the variability of the time series across different sub-series in the time series that are homogenous. The definition of this sub-series is slightly subjective but, usually, we can safely assume the sub-series as the seasonal length. Therefore, what we will be trying to minimize is the variability of the time series across different seasonality cycles.

Reference check:

The research paper proposing Guerrero’s method is cited in the References section under reference 4.

There are stark differences in the way both these optimization methods work and we need to be careful when using them. If our main concern is to remove the heteroscedastic behavior of the time series, Guerrero’s method is what we can use.

We have made a transformer available in src.transforms.target_transformations.py called BoxCoxTransformer. There are a few parameters and properties that we need to be aware of:

The core logic implemented in the Transformer is as follows:

## Fit Process
# Add one if needed
y = self._add_one(y)
# Find optimum box cox lamda if optimization is Guerrero
self.boxcox_lambda = self._optimize_lambda(y)
## Transform Process
boxcox(y.values, lmbda=self.boxcox_lambda)
## Inverse Transform
self._subtract_one(inv_boxcox(y.values, self.boxcox_lambda))

Now, let’s see how we can use it:

from src.transforms.target_transformations import BoxCoxTransformer
boxcox_transformer = BoxCoxTransformer()
y_boxcox = boxcox _transformer.fit_transform(y)

y_boxcox will contain the Box-Cox transformed time series. To get back to the original time series, we can use the inverse_transform function.

Both Box-Cox and Log Transform can be used for correcting heteroscedasticity. But, as mentioned before, log transform is a strong transformation and Box-Cox gives us another lever to tweak and tune the transformation to suit our data. We can look at Box-Cox as a flexible log transform, which can be tuned to make the right transformation for our data. Do check out the notebook, where you can see and play around with these different transformations and get a feel of what it will do to your data.

When we approach the forecasting problem at scale, we will have hundreds, thousands, or millions of time series that we will need to analyze before forecasting. In such scenarios, an AutoML approach is needed to be practical.

AutoML approach to target transformation

So far, we have discussed many ways to make a series more stationary (we are using the word stationary here in the non-mathematical sense), such as detrending, deseasonalizing, differencing, and monotonic transformations. We’ve also looked at statistical tests to check whether trends, seasonality, and so on are present in a time series. So, the natural next step is to put it all together to carry out these transforms in an automated way while choosing good defaults wherever possible. This is exactly what we did and implemented an AutoStationaryTransformer in src.transforms.target_transformations.

The following flow chart explains the logic of this in an automated way:

Figure 7.5 – Flow chart for AutoStationaryTransformer

Figure 7.6: Flow chart for AutoStationaryTransformer

We have excluded differencing from this implementation for two reasons:

Now, let’s see what parameters we can use to tweak AutoStationaryTransformer:

Let’s apply this AutoStationaryTransformer to a synthetic time series and see how well it works (full code in the associated notebook):

from src.transforms.target_transformations import AutoStationaryTransformer
auto_stationary = AutoStationaryTransformer(seasonal_period=25)
y_stat = auto_stationary.fit_transform(y_final)

Fig 7.7: AutoStationaryTransformer—Before and after

We can see that the AutoStationaryTransformer has deseasonalized and de-trended the time series. In this particular example, Detrending, Deseasonalizing, and Box-Cox Transformation were applied by the AutoStationaryTransformer.

Now, let’s apply this automatic transformation to the dataset we have been working with:

train_df = pd.read_parquet(preprocessed/"selected_blocks_train_missing_imputed_feature_engg.parquet")
transformer_pipelines = {}
for _id in tqdm(train_df["LCLid"].unique()):
    #Initialize the AutoStationaryTransformer with a seasonality period of 48*7
    auto_stationary = AutoStationaryTransformer(seasonal_period=48*7)
    #Creating the timeseries with datetime index
    y = train_df.loc[train_df["LCLid"]==_id, ["energy_consumption","timestamp"]].set_index("timestamp")
    #Fitting and transforming the train
    y_stat = auto_stationary.fit_transform(y, freq="30min")
    # Setting the transformed series back to the dataframe
    train_df.loc[train_df["LCLid"]==_id, "energy_consumption"] = y_stat.values
    #Saving the pipeline
    transformer_pipelines[_id] = auto_stationary

The code to execute this is split into two notebooks called 02-Dealing_with_Non-Stationarity.ipynb and 02a-Dealing_with_Non-Stationarity-Train+Val.ipynb in the Chapter06 folder. The former does the auto-stationary transformation on the train data, while the latter does it on train and validation data combined. This is to simulate how we would predict for validation data (by just using train data for training) and for test data (where we use the train and validation data for training).

This process is slightly time-consuming. I suggest that you run the notebook, grab lunch or a snack, and come back. Once it’s done, the 02-Dealing_with_Non-Stationarity.ipynb notebook will save a couple of files:

The 02a-Dealing_with_Non-Stationarity-Train+Val.ipynb notebook also saves the corresponding files for the train and validation datasets.

The dataset we are working on has almost negligible trends and is pretty stationary throughout. The impact of these transformations will be more evident in time series with strong trends and heteroscedasticity.

Best practice:

This kind of explicit detrending and deseasonalizing before modeling can also be seen as a form of boosting. This should be considered as just another alternative to modeling all of this together. There can be situations where letting the model learn from end to end in a data-driven manner performs better than injecting these strong inductive biases using explicit detrending and deseasonalization and vice versa. Cross-validated test scores should always have the last word.

Congratulations on making it through a heavy chapter full of new concepts, some statistics, and mathematics. From the point of view of applying ML models for time series, the concepts in this chapter will be really helpful in taking your models to the next level.

Summary

After getting down to a practical level in the previous chapter, we stayed there and plowed on to review concepts such as stationarity and how to deal with such non-stationary time series. We learned about techniques we can use to explicitly handle non-stationary time series, such as differencing, detrending, deseasonalizing, and so on. To put this all together, we saw an automatic way of transforming the target, learned how to use the implementation provided, and applied it to our dataset. Now that we have the necessary skills to effectively transform a time series into an ML dataset, in the next chapter, we will start applying a few ML models to the dataset using the features we’ve created.

References

The following are the references for this chapter:

  1. Bayazit, M. and Önöz, B. (2007), To prewhiten or not to prewhiten in trend analysis?, Hydrological Sciences Journal, 52:4, 611–624. https://doi.org/10.1623/hysj.52.4.611.
  2. White, H. (1980), A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity. Econometrica Vol. 48, No. 4 (May 1980), pp. 817–838 (22 pages). https://doi.org/10.2307/1912934.
  3. Box, G. E. P. and Cox, D. R. (1964), An analysis of transformations. Journal of the Royal Statistical Society, Series B, 26, 211–252. http://www.ime.usp.br/~abe/lista/pdfQWaCMboK68.pdf.
  4. Guerrero, Victor M. (1993), Time-series analysis supported by power transformations. Journal of Forecasting, Volume 12, Issue 1, 37–48. https://onlinelibrary.wiley.com/doi/10.1002/for.3980120104.

Further reading

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