Kernel-based Time-varying Regression - Part I

Kernel-based time-varying regression (KTR) is a time series model to address

1. time-varying regression coefficients

2. complex seasonality pattern

The full details of the model structure with an application in marketing media mix modeling can be found in Ng, Wang and Dai (2021). The core of KTR is the use of latent variables to define a smooth time varying representation of model coefficients, which bears a similar ideas in Kernel Smoothsing. 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 the computation speed, and allows for handling of high dimensional data and detecting small variances.

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

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)

/Users/towinazure/edwinnglabs/orbit/orbit/utils/dataset.py:112: UserWarning: Parsing dates in DD/MM/YYYY format when dayfirst=False (the default) was specified. This may lead to inconsistently parsed dates! Specify a format to ensure consistent parsing.
  df["date"] = pd.date_range(

	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)

2023-01-20 18:01:49 - orbit - INFO - Optimizing (PyStan) with algorithm: LBFGS.
2023-01-20 18:01:49 - orbit - INFO - Using SVI (Pyro) with steps: 301, samples: 100, learning rate: 0.1, learning_rate_total_decay: 1.0 and particles: 100.
2023-01-20 18:01:49 - orbit - INFO - step    0 loss = -1946.6, scale = 0.093118
INFO:orbit:step    0 loss = -1946.6, scale = 0.093118
2023-01-20 18:01:53 - orbit - INFO - step  100 loss = -3126.6, scale = 0.0099929
INFO:orbit:step  100 loss = -3126.6, scale = 0.0099929
2023-01-20 18:01:56 - orbit - INFO - step  200 loss = -3146, scale = 0.0099024
INFO:orbit:step  200 loss = -3146, scale = 0.0099024
2023-01-20 18:02:00 - orbit - INFO - step  300 loss = -3147.7, scale = 0.0098069
INFO:orbit:step  300 loss = -3147.7, scale = 0.0098069

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

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.43009	9.56284	9.69732
1	2000-01-02	9.42618	9.56036	9.69355
2	2000-01-03	9.42676	9.56115	9.69353
3	2000-01-04	9.42838	9.56113	9.69383
4	2000-01-05	9.42540	9.56119	9.69603

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)

2023-01-20 18:02:04 - orbit - INFO - Optimizing (PyStan) with algorithm: LBFGS.
INFO:orbit:Optimizing (PyStan) with algorithm: LBFGS.
2023-01-20 18:02:05 - orbit - INFO - Using SVI (Pyro) with steps: 301, samples: 100, learning rate: 0.1, learning_rate_total_decay: 1.0 and particles: 100.
INFO:orbit:Using SVI (Pyro) with steps: 301, samples: 100, learning rate: 0.1, learning_rate_total_decay: 1.0 and particles: 100.
2023-01-20 18:02:05 - orbit - INFO - step    0 loss = -2190.9, scale = 0.093667
INFO:orbit:step    0 loss = -2190.9, scale = 0.093667
2023-01-20 18:02:08 - orbit - INFO - step  100 loss = -4339.3, scale = 0.007011
INFO:orbit:step  100 loss = -4339.3, scale = 0.007011
2023-01-20 18:02:11 - orbit - INFO - step  200 loss = -4362.5, scale = 0.0070783
INFO:orbit:step  200 loss = -4362.5, scale = 0.0070783
2023-01-20 18:02:14 - orbit - INFO - step  300 loss = -4296.8, scale = 0.0072582
INFO:orbit:step  300 loss = -4296.8, scale = 0.0072582

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

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.49836	9.57769	9.65776	9.46910	9.54843	9.62850	0.00000	0.00000	0.00000	-0.02756	-0.02756	-0.02756	0.05683	0.05683	0.05683
1	2000-01-02	9.44648	9.52624	9.60597	9.46715	9.54692	9.62665	0.00000	0.00000	0.00000	-0.07838	-0.07838	-0.07838	0.05770	0.05770	0.05770
2	2000-01-03	9.50748	9.58766	9.66673	9.46750	9.54769	9.62675	0.00000	0.00000	0.00000	-0.01856	-0.01856	-0.01856	0.05854	0.05854	0.05854
3	2000-01-04	9.56713	9.64639	9.72560	9.46843	9.54769	9.62690	0.00000	0.00000	0.00000	0.03936	0.03936	0.03936	0.05934	0.05934	0.05934
4	2000-01-05	9.55023	9.63102	9.71150	9.46704	9.54783	9.62831	0.00000	0.00000	0.00000	0.02310	0.02310	0.02310	0.06009	0.06009	0.06009

Tips: there is an additional arg seasonality_fs_order to control the number of orders in fourier series terms we want to approximate the seasonality. In general, they cannot violate the condition that \(2 \times \text{fourier series order} < \text{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