
Introduction
I came across a new and promising Python Library for Time Series – Sktime. It provides a plethora of Time Series Functionalities like Transformations, Forecasting algorithms, the Composition of Forecasters, Model Validation, Pipelining the entire flow, and many more. In this article, we explore some features the library provides, the most important one being how to make a Machine Learning model – a Light GBM fit for time series forecasting.
Abstract
When it comes to Time Series Forecasting, ARIMA and its variants dominate the domain(simple yet powerful methods). However, having a strong personal liking for Ensemble Tree Models it is always tempting to use them for forecasting too! We explore ways in which we can formulate a Time Series forecasting problem into a Supervised Learning Task and solve it using Sktime.
Forecasting as Supervised Learning
Single-Step Forecasting
We take the previous k time-step values as our regressors and output the k+1 value. This is a One-Step Ahead forecaster, which can be better understood from the below diagram. In the diagram, we auto-regress on the previous four values and forecast the fifth one at each point in time. If you are familiar with Sklearn APIs, we can fit our regression model by using fit(X,y), where X and y are as shown below.

Assuming we use a GBM Regressor, the below diagram fits well to illustrate the concept!

But the real-world problems involve forecasting the future multiple steps ahead. How do we tackle that?
There are multiple ways to do that, for further reading, please refer to Reference[1]. We explore one of the strategies, namely the Recursive variant.
Multi-Step Forecasting
We achieve this by iterating over multiple-steps(which is known as Forecasting Horizon) and using the forecasted value as input for forecasting the next value. This is known as Recursive Strategy and despite some of its limitations(error sensitivity), it works well in real-world settings. The below diagram illustrates the idea.

A pseudo-code for the same:
for each_time_step in forecasting_horizon:
Use forecast from (each_time_step - 1) step in the Input.
forecast()

Implementation Using Sktime
Let’s start by installing Sktime and importing the libraries!
! pip install sktime==0.4.3
import pandas as pd
import numpy as np
import seaborn as sns
import warnings
import itertools
import numpy as np
import matplotlib.pyplot as plt
import Lightgbm as lgb
from pylab import rcParams
rcParams['figure.figsize'] = 18, 8
warnings.filterwarnings("ignore")
from statsmodels.tsa.seasonal import seasonal_decompose
#-----------------------Imports from Sktime-------------------------
from sktime.forecasting.base import ForecastingHorizon
from sktime.transformers.single_series.detrend import Deseasonalizer, Detrender
from sktime.forecasting.trend import PolynomialTrendForecaster
from sktime.forecasting.model_selection import (
temporal_train_test_split,
)
from sktime.utils.plotting import plot_series
from sktime.forecasting.compose import (
TransformedTargetForecaster,
ReducedRegressionForecaster
)
The Sales data that we will work on can be found on Kaggle. Assuming the data is in the project root, we load the data.
#------------------------Reading Data from CSV----------------------
df = pd.read_csv("train.csv",parse_dates=['date'])
df.head()

For our implementation and demonstration, we select the Data from "store 1" and aggregate the data at a monthly level.
#--------------------------Data Preparation-------------------------
#For the sake of demonstration, we will train our model on monthly aggregated Sales data of a particular store
# Select sales for Store 1 Only.
store1_agg = df.loc[df['store']==1].groupby(['date'])['sales'].sum()
store1_agg.index = pd.to_datetime(store1_agg.index)
#Aggregate the Data on a Monthly basis.
store1_agg_monthly = store1_agg.resample('M').sum()
#--------------------Visulaize Data on a Time Plot------------------
sns.lineplot(
data=store1_agg_monthly,
)
plt.title("Store-1 Sales Data aggreagted at Month Level")
plt.show()

Time Series Decomposition
The above plot gives us the idea that the time series is seasonal and heteroskedastic(the variations increase with time for a Season). We need to use a Multiplicative Decomposition in such a case, to examine its Seasonal, Trend components and the residuals. The Seasonality is Annual, so we decompose with a Seasonality hyper-parameter set to 12. We notice a Linear(approximately) Trend and infer the Seasonal Element from the plot below. The maximum sales occur during the Mid-Year.
#----------------------Time Series Decomposition--------------------
seasonal_decompose(store1_agg_monthly,model="multiplicative",period=12).plot()
plt.show()

Sktime – Data Splitting
We go for a train-test split next using Sktime‘s specialized function – temporal_train_test_split(20% of the data is set for validation). This function is tailored for Sequential Data Splitting as it maintains the ordering of the data without shuffling across time.
#--------------------Time Series Train-Test split-------------------
store1_agg_monthly.index = store1_agg_monthly.index.to_period('M')
y_train, y_test = temporal_train_test_split(store1_agg_monthly, test_size=0.2)
Sktime – Trend Forecaster and Detrending
What if we just want to forecast the Trend for analysis purposes or in order to use it as one of the steps in Composition Forecasting?
Sktime provides PolynomialTrendForecaster for this purpose, which can be used to fit a polynomial trend. The forecast it generates can be in turn combined with forecasts from other components to generate a final forecast. To illustrate its effects, we can use Sktime’s Detrender which takes the Polynomial forecaster as input in order to detrend the Time Series. It returns the residuals after removing the trend generated by the Trend Forecaster. The trend generated can be seen in the yellow line in the below plot. As we see the residuals consist mostly of the Seasonal component.
#--------------------------Detrender-----------------------------
#degree=1 for Linear
forecaster = PolynomialTrendForecaster(degree=1)
transformer = Detrender(forecaster=forecaster)
#Get the residuals after fitting a linear trend
y_resid = transformer.fit_transform(y_train)
# Internally, the Detrender uses the in-sample predictions
# of the PolynomialTrendForecaster
forecaster = PolynomialTrendForecaster(degree=1)
fh_ins = -np.arange(len(y_train)) # in-sample forecasting horizon
y_pred = forecaster.fit(y_train).predict(fh=fh_ins)
plot_series(y_train, y_pred, y_resid, labels=["y_train", "fitted linear trend", "residuals"]);

Sktime – Deseasonalizer
Similarly, in order to extract the seasonal component from the time series, we can use Sktime’s Deseasonalizer. This removes the seasonal component leaving behind the trend and other cyclic components in the time series. We extract the Annual Seasonality in this step(period=12).
#--------------------------Deseasonalizer---------------------------
#Multiplicative Deseasonalizer, period = 12(for Monthly Data)
deseasonalizer = Deseasonalizer(model="multiplicative", sp=12)
plot_series(deseasonalizer.fit_transform(y_train))
seasonal = deseasonalizer.fit_transform(y_train)

Sktime – Forecasting using Recursive Strategy
To implement a recursive strategy(as discussed above) needs to take care of various boundary conditions and the implementation might not be a robust one. This implementation is abstracted by Sktime’s RecursiveRegressionForecaster or a more generalized ReducedRegressionForecaster, which takes a Sklearn API compatible Regressor and reduces it into a Multi-Step Ahead Forecaster based on Recursive strategy. Here we auto-regress on the previous four historic values.
Some important hyperparameters for ReducedRegressionForecaster:
- _windowlength: The number of (immediate)previous historic values to consider as regressors.
- strategy: Multi-Step Forecast strategy, in this case, its "recursive".
regressor = lgb.LGBMRegressor()
forecaster = ReducedRegressionForecaster(
regressor=regressor, window_length=4, strategy="recursive" #hyper-paramter to set recursive strategy
)
Sktime -Meta Forecasting
Finally, we can use a Meta-Forecaster which can combine the effects of multiple forecasters to produce the final forecast. To exemplify its use, we take the following case. We use three different forecasters and combine them:
- Seasonality Forecaster(Using Deseasonalizer)
- Trend Forecaster(Using Detrender)
- AutoRegressor, based on immediate previous values(Using LightGBM and ReducedRegressionForecaster)
We use Sktime’s TransformedTargetForecaster for this purpose which abstracts the meta forecasting techniques. Internally it not only fits multiple forecasters as a pipeline but also generates the forecasts and combines them(the seasonal, the trend, and the autoregressive components) into a final one. The pipeline can be found in the below code snippet. All the sequential steps conveniently pipelined using TransformedTargetForecaster.
#----------------------------Create Pipeline--------------------
def get_transformed_target_forecaster(alpha,params):
#Initialize Light GBM Regressor
regressor = lgb.LGBMRegressor(alpha = alpha,**params)
#-----------------------Forecaster Pipeline-----------------
#1.Separate the Seasonal Component.
#2.Fit a forecaster for the trend.
#3.Fit a Autoregressor to the resdiual(autoregressing on four historic values).
forecaster = TransformedTargetForecaster(
[
("deseasonalise", Deseasonalizer(model="multiplicative", sp=12)),
("detrend", Detrender(forecaster=PolynomialTrendForecaster(degree=1))),
(
# Recursive strategy for Multi-Step Ahead Forecast.
# Auto Regress on four previous values
"forecast",
ReducedRegressionForecaster(
regressor=regressor, window_length=4, strategy="recursive",
),
),
]
)
return forecaster
Uncertainty Estimation
Real-life forecasts involve uncertainty estimations mostly in the form of Confidence/Prediction Intervals. We use LightGBM’s quantile loss to get regressed values based on quantiles(10th, median, and 90th being popular ones). A blog on the same can be found here.
#-------------------Fitting an Auto Regressive Light-GBM------------
#Setting Quantile Regression Hyper-parameter.
params = {
'objective':'quantile'
}
#A 10 percent and 90 percent prediction interval(0.1,0.9 respectively).
quantiles = [.1, .5, .9] #Hyper-parameter "alpha" in Light GBM
#Capture forecasts for 10th/median/90th quantile, respectively.
forecasts = []
#Iterate for each quantile.
for alpha in quantiles:
forecaster = get_transformed_target_forecaster(alpha,params)
#Initialize ForecastingHorizon class to specify the horizon of forecast
fh = ForecastingHorizon(y_test.index, is_relative=False)
#Fit on Training data.
forecaster.fit(y_train)
#Forecast the values.
y_pred = forecaster.predict(fh)
#List of forecasts made for each quantile.
y_pred.index.name="date"
y_pred.name=f"predicted_sales_q_{alpha}"
forecasts.append(y_pred)
#Append the actual data for plotting.
store1_agg_monthly.index.name = "date"
store1_agg_monthly.name = "original"
forecasts.append(store1_agg_monthly)
Finally, we plot the results-
#-------------------Final Plotting of Forecasts------------------
plot_data = pd.melt(pd.concat(forecasts,axis=1).reset_index(), id_vars=['date'],
value_vars=['predicted_sales_q_0.1', 'predicted_sales_q_0.5',
'predicted_sales_q_0.9','original'])
plot_data['date'] = pd.to_datetime(plot_data['date'].astype(str).to_numpy())
plot_data['if_original'] = plot_data['variable'].apply(
lambda r:'original' if r=='original' else 'predicted'
)
sns.lineplot(data = plot_data,
x='date',
y='value',
hue='if_original',
style="if_original",
markers=['o','o'],
)
plt.title("Final Forecast")
plt.show()

Conclusion
The library seems to be under active development so it’s good to use the latest releases(and the corresponding functions/modules). Overall, this seems to be a very promising Time Series package and worthwhile to explore. To download the notebook please refer to this Kaggle Kernel. Thanks for the read!
References
- Bontempi, Gianluca & Ben Taieb, Souhaib & Le Borgne, Yann-Aël. ( 2013). Machine Learning Strategies for Time Series Forecasting – https://www.researchgate.net/publication/236941795_Machine_Learning_Strategies_for_Time_Series_Forecasting.
- Markus Löning, Anthony Bagnall, Sajaysurya Ganesh, Viktor Kazakov, Jason Lines, Franz Király (2019): "sktime: A Unified Interface for Machine Learning with Time Series"
- LightGBM-Quantile loss – https://towardsdatascience.com/lightgbm-for-quantile-regression-4288d0bb23fd