Kernel-based Time-varying Regression - Part I

To address many of the time-series modeling problems at Uber, for example:

1. time-varying regression coefficients

2. 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

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)

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,

\[y_t = l_t + s_t + r_t +\epsilon_t, ~ t = 1,\cdots, T,\]

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

\[y_t = X_t^T\beta_t + \epsilon_t, ~ t = 1,\cdots, T,\]

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

\[B = K b^T\]

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,

\[l_t = \beta_{t, \text{lev}}\]

It can also be viewed as a dynamic intercept (where the regressor is a vector of ones).

For the seasonality,

\[B_{\text{seas}} = K_{\text{seas}}b^T_{\text{seas}}\]

\[s_t = X_{t, \text{seas}}\beta_{t, \text{seas}}\]

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.

# 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)

	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

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.

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.

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,
)

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

<orbit.forecaster.svi.SVIForecaster at 0x143c8e110>

We can take a look how the level is fitted with the data.

predicted_df = ktr.predict(df=df)
predicted_df.head()

	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.

ktr.get_posterior_samples().keys()

dict_keys(['lev_knot', 'lev', 'yhat', 'obs_scale'])

_ = 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.

_ = 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.

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,
)

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

<orbit.forecaster.svi.SVIForecaster at 0x144404550>

predicted_df = ktr_with_seas.predict(df=df, decompose=True)

predicted_df.head(5)

	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.

_ = 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)

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

References

1. Ng, Wang and Dai (2021) Bayesian Time Varying Coefficient Model with Applications to Marketing Mix Modeling, arXiv preprint arXiv:2106.03322

2. Hastie, Trevor and Tibshirani, Robert. (1990), Generalized Additive Models, New York: Chapman and Hall.

3. Wood, S. N. (2006), Generalized Additive Models: an introduction with R, Boca Raton: Chapman & Hall/CRC

4. Harvey, C. A. (1989). Forecasting, Structural Time Series and the Kalman Filter, Cambridge University Press.

5. Durbin, J., Koopman, S. J.. (2001). Time Series Analysis by State Space Methods, Oxford Statistical Science Series