Kernel-based Time-varying Regression - Part I¶
To address many of the time-series modeling problems at Uber, for example:
time-varying regression coefficients
complex seasonality pattern
This work proposes a kernel-based time-varying regression (KTR). The full details of model structure with an example application in marketing data science are found in Ng, Wang and Dai (2021). The core of KTR is that latent variables along with their own time-point location to approximate the coefficient curves varied over time. The core of KTR is the use of latent variables to define a smooth time varying representation of model coefficients. These smooth representations are Kernel Smooths. The KTR approach sharply reduces the number of parameters compared to typical dynamic linear models such as Harvey (1989) and Durbin and Koopman (2002). The reduced number of parameters improves computation speed; this allows for KTR tp practically handle larger numbers of regressors and detect small variance as explained in Ng (2021).
To topics covered here in Part I, are 1. KTR model structure 2. syntax to initialize, fit and predict a model 3. fit a model with complex seasonality 4. visualization of prediction and decomposed components
[1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import orbit
from orbit.models import KTR
from orbit.diagnostics.plot import plot_predicted_data, plot_predicted_components
from orbit.utils.dataset import load_electricity_demand
%matplotlib inline
pd.set_option('display.float_format', lambda x: '%.5f' % x)
[2]:
print(orbit.__version__)
1.1.3dev
Model Structure¶
This section gives the mathematical structure of the KTR model. In short, it considers a time-series (\(y_t\)) as the linear combination of three parts which are the local-trend (\(l_t\)), seasonality (\(s_t\)), and regression (\(r_t\)) terms. Mathematically,
where the \(\epsilon_t\) comprise a stationary random error process.
In KTR, the distinction between the local-trend, seasonality, and regressors while useful is semi-arbitrary and the time-series can also be considered as
where \(\beta_t\) is a \(P\)-dimensional vector of coefficients that vary over time (i.e., \(\beta_i\) is almost certainly different from \(\beta_j\) for \(i \neq j\)) and \(X_t\) \(P\)-dimensional covariate vector (i.e., the \(t\)th row of \(X\), the design matrix).
To reduce the total number of parameters in the model (potentially \(P \times T\)) the \(\beta_t\) are parameterized with a weighted sum of \(J\) local latent variables (\(b_1,..b_j,..b_J\)). That is
where - coefficient matrix \(B\) has size \(T \times P\) with rows equal to the \(\beta_t\). - knot matrix \(b\) with size \(P\times J\); each entry is a latent variable \(b_{p, j}\). The \(b_j\) can be viewed as the “knots” from the perspective of spline regression and \(j\) is a time index such that \(t_j \in [1, \cdots, T]\). - kernel matrix \(K\) with size \(T\times J\) where the \(i\)th row and \(j\)th element can be viewed as the normalized weight \(k(t_j, t) / \sum_{j=1}^{J} k(t_j, t)\)
For the level/trend,
It can also be viewed as a dynamic intercept (where the regressor is a vector of ones).
For the seasonality,
We use Fourier series to handle the seasonality; i.e., sin waves with varying periods are used for the columns of \(X_{ \text{seas}}\).
The details for the additional regressors are given in Part II, as they are not used in this tutorial. Note this includes different choices of kernel function (which determines the kernel matrix \(K\)) and prior for matrix \(b\).
Data¶
To illustrate the usage of KTR, consider the daily series of electricity demand in Turkey from the 2000 - 2008.
[3]:
# from 2000-01-01 to 2008-12-31
df = load_electricity_demand()
date_col = 'date'
response_col = 'electricity'
df[response_col] = np.log(df[response_col])
print(df.shape)
df.head()
(3288, 2)
[3]:
date | electricity | |
---|---|---|
0 | 2000-01-01 | 9.43760 |
1 | 2000-01-02 | 9.50130 |
2 | 2000-01-03 | 9.63565 |
3 | 2000-01-04 | 9.65392 |
4 | 2000-01-05 | 9.66089 |
[4]:
print(f'starts with {df[date_col].min()}\nends with {df[date_col].max()}\nshape: {df.shape}')
starts with 2000-01-01 00:00:00
ends with 2008-12-31 00:00:00
shape: (3288, 2)
Train / Test Split¶
Split the data into a training set and test set for model validation.
[5]:
test_size=365
train_df=df[:-test_size]
test_df=df[-test_size:]
A Quick Start on KTR¶
Here the Similar to other model types in Orbit, KTR follows sklearn model API style. First an instance of the Orbit class KTR
is created. Second fit and predict methods are called for that instance. Note that unlike version <=1.0.15
, the fitting API arg are within the function; thus, KTR
is called directly.
[6]:
ktr = KTR(
response_col=response_col,
date_col=date_col,
seed=2021,
estimator='pyro-svi',
# bootstrap sampling to capture uncertainties
n_bootstrap_draws=1e4,
# pyro training config
num_steps=301,
message=100,
)
[7]:
ktr.fit(train_df)
INFO:orbit:Optimizing(PyStan) with algorithm:LBFGS .
INFO:orbit:Using SVI(Pyro) with steps:301 , samples:100 , learning rate:0.1, learning_rate_total_decay:1.0 and particles:100.
INFO:root:Guessed max_plate_nesting = 1
INFO:orbit:step 0 loss = -1946.6, scale = 0.092812
INFO:orbit:step 100 loss = -3132.3, scale = 0.0099682
INFO:orbit:step 200 loss = -3140.4, scale = 0.0098539
INFO:orbit:step 300 loss = -3149.5, scale = 0.0097989
[7]:
<orbit.forecaster.svi.SVIForecaster at 0x143c8e110>
We can take a look how the level is fitted with the data.
[8]:
predicted_df = ktr.predict(df=df)
predicted_df.head()
[8]:
date | prediction_5 | prediction | prediction_95 | |
---|---|---|---|---|
0 | 2000-01-01 | 9.46082 | 9.59296 | 9.72753 |
1 | 2000-01-02 | 9.45634 | 9.59054 | 9.72242 |
2 | 2000-01-03 | 9.45696 | 9.59080 | 9.72278 |
3 | 2000-01-04 | 9.45798 | 9.59035 | 9.72344 |
4 | 2000-01-05 | 9.45498 | 9.59037 | 9.72499 |
One can use .get_posterior_samples()
to extract the samples for all sampling parameters.
[9]:
ktr.get_posterior_samples().keys()
[9]:
dict_keys(['lev_knot', 'lev', 'yhat', 'obs_scale'])
[10]:
_ = plot_predicted_data(training_actual_df=train_df, predicted_df=predicted_df,
date_col=date_col, actual_col=response_col,
test_actual_df=test_df, markersize=20, lw=.5)

It can also be helpful to see the trend knot locations and levels. This is done with the plot_lev_knots
function.
[11]:
_ = ktr.plot_lev_knots()

Fitting with Complex Seasonality¶
The previous model fit is not satisfactory as there is clear seasonality in the electrical demand time-series that is not accounted for. In this modelling example the electrical demand data is fit with a dual seasonality for weekly and yearly patterns. Since the data is daily, the seasonality periods are 7 and 365.25. These are added into the KTR object as a list through the seasonality
arg. Otherwise the process is the same as the previous example.
[12]:
ktr_with_seas = KTR(
response_col=response_col,
date_col=date_col,
seed=2021,
seasonality=[7, 365.25],
estimator='pyro-svi',
n_bootstrap_draws=1e4,
# pyro training config
num_steps=301,
message=100,
)
[13]:
ktr_with_seas.fit(train_df)
INFO:orbit:Optimizing(PyStan) with algorithm:LBFGS .
INFO:orbit:Using SVI(Pyro) with steps:301 , samples:100 , learning rate:0.1, learning_rate_total_decay:1.0 and particles:100.
INFO:root:Guessed max_plate_nesting = 1
INFO:orbit:step 0 loss = -2190.9, scale = 0.093359
INFO:orbit:step 100 loss = -4301.4, scale = 0.0069362
INFO:orbit:step 200 loss = -4296.9, scale = 0.0071429
INFO:orbit:step 300 loss = -4368.3, scale = 0.0069215
[13]:
<orbit.forecaster.svi.SVIForecaster at 0x144404550>
[14]:
predicted_df = ktr_with_seas.predict(df=df, decompose=True)
[15]:
predicted_df.head(5)
[15]:
date | prediction_5 | prediction | prediction_95 | trend_5 | trend | trend_95 | regression_5 | regression | regression_95 | seasonality_7_5 | seasonality_7 | seasonality_7_95 | seasonality_365.25_5 | seasonality_365.25 | seasonality_365.25_95 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2000-01-01 | 9.54501 | 9.62752 | 9.71168 | 9.51590 | 9.59841 | 9.68258 | 0.00000 | 0.00000 | 0.00000 | -0.02766 | -0.02766 | -0.02766 | 0.05676 | 0.05676 | 0.05676 |
1 | 2000-01-02 | 9.49290 | 9.57605 | 9.65914 | 9.51370 | 9.59685 | 9.67994 | 0.00000 | 0.00000 | 0.00000 | -0.07844 | -0.07844 | -0.07844 | 0.05764 | 0.05764 | 0.05764 |
2 | 2000-01-03 | 9.55393 | 9.63741 | 9.72024 | 9.51391 | 9.59739 | 9.68022 | 0.00000 | 0.00000 | 0.00000 | -0.01846 | -0.01846 | -0.01846 | 0.05848 | 0.05848 | 0.05848 |
3 | 2000-01-04 | 9.61321 | 9.69599 | 9.77808 | 9.51450 | 9.59727 | 9.67937 | 0.00000 | 0.00000 | 0.00000 | 0.03944 | 0.03944 | 0.03944 | 0.05927 | 0.05927 | 0.05927 |
4 | 2000-01-05 | 9.59604 | 9.68016 | 9.76398 | 9.51296 | 9.59708 | 9.68090 | 0.00000 | 0.00000 | 0.00000 | 0.02305 | 0.02305 | 0.02305 | 0.06003 | 0.06003 | 0.06003 |
Tips: there is an additional arg seasonality_fs_order
to control the number of orders we want to approximate the seasonality. In general, they cannot violate the condition that 2 * seasonality_fs_order
< seasonality
since each order represents adding a pair of sine and cosine regressors.
More Diagnostic and Visualization¶
Here are a few more diagnostic and visualization. The fit is decomposed into components, the local trend and both periods of seasonality.
[16]:
_ = plot_predicted_data(training_actual_df=train_df, predicted_df=predicted_df,
date_col=date_col, actual_col=response_col,
test_actual_df=test_df, markersize=10, lw=.5)

[17]:
_ = plot_predicted_components(predicted_df=predicted_df, date_col=date_col, plot_components=['trend', 'seasonality_7', 'seasonality_365.25'])

References¶
Ng, Wang and Dai (2021) Bayesian Time Varying Coefficient Model with Applications to Marketing Mix Modeling, arXiv preprint arXiv:2106.03322
Hastie, Trevor and Tibshirani, Robert. (1990), Generalized Additive Models, New York: Chapman and Hall.
Wood, S. N. (2006), Generalized Additive Models: an introduction with R, Boca Raton: Chapman & Hall/CRC
Harvey, C. A. (1989). Forecasting, Structural Time Series and the Kalman Filter, Cambridge University Press.
Durbin, J., Koopman, S. J.. (2001). Time Series Analysis by State Space Methods, Oxford Statistical Science Series