#  VQ-AR: Vector Quantized Autoregressive Probabilistic Time Series
Forecasting

Kashif Rasul
<EMAIL_ADDRESS>
&Young-Jin Park
NAVER CLOVA & NAVER AI Lab
<EMAIL_ADDRESS>
&Max Nihlén Ramström
NAVER CLOVA
<EMAIL_ADDRESS>&Kyung-Min Kim
NAVER CLOVA & NAVER AI Lab
<EMAIL_ADDRESS>
###### Abstract

Time series models aim for accurate predictions of the future given the past,
where the forecasts are used for important downstream tasks like business
decision making. In practice, deep learning based time series models come in
many forms, but at a high level learn some continuous representation of the
past and use it to output point or probabilistic forecasts. In this paper, we
introduce a novel autoregressive architecture, VQ-AR, which instead learns a
_discrete_ set of representations that are used to predict the future.
Extensive empirical comparison with other competitive deep learning models
shows that surprisingly such a discrete set of representations gives state-of-
the-art or equivalent results on a wide variety of time series datasets. We
also highlight the shortcomings of this approach, explore its zero-shot
generalization capabilities, and present an ablation study on the number of
representations. The full source code of the method will be available at the
time of publication with the hope that researchers can further investigate
this important but overlooked inductive bias for the time series domain.

##  1 Introduction

Time series forecasting is an important business and scientific machine
learning problem which is typically used to support decision making for down-
stream tasks. Classical methods like in [22] model each time series in a
dataset individually using hand-crafted features and are used extensively as
they providing strong baselines.

Recently deep learning based methods offer an alternative approach by
utilizing a shared model trained over the whole dataset of time series,
without the need for explicit feature engineering. These methods are
competitive, especially when dealing with a large amount of data, against the
classical methods, but can also suffer from overfitting or learn spurious
correlations due to their flexibility.

Often the forecasts are point valued, however it is much more useful for the
decision making process to also incorporate the inherent (aleatoric)
uncertainty present in the observations via _probabilistic forecasts_. In the
univariate setting this is usually done by learning the parameters of some
chosen distribution, or via quantile regression or non-parametric and semi-
parametric models which explicitly learn the conditional quantiles of the next
time step given the past [39]. Other approaches include using Conformal
Predictions [49, 54]. In the multivariate setting [50] however, one needs to
resort to some approximation of the full multivariate distribution [47] due to
computational tractability or learn the conditional distribution via
Normalizing Flows [46, 10], Energy Based Methods [45], or Generative
Adversarial Networks (GANs) [30].

As in most supervised deep learning methods, these models learn some
representation of the history of a time series in order to best forecast,
which can be thought of as a kind of self-supervised learning task (not unlike
GPT [7] or Bert [11] in the Natural Language Processing (NLP) setting). Thus
in order to provide better historical representations we have a large
collection of architectural choices we can make for the problem or dataset at
hand.

The core insight of this work is the realization that there might only be a
discrete collection of representations of time series histories which could
potentially be used for the prediction task. The success of modern NLP on the
back of discrete token representations leads us to consider if something
similar would also work in the time series domain. Thus, in this paper we take
inspiration from the Vector Quantized-Variational AutoEncoder (VQ-VAE) [51]
model to learn a set of _latent_ representations in order to best forecast in
an autoregressive fashion.

Thus, the main contribution of this paper is a novel autoregressive time
series forecasting method which is trained end-to-end and provides excellent
inductive bias for the probabilistic forecasting task. We highlight the
model’s performance and properties via extensive experimentation and show that
it achieves comparable results against ensembles of tree-based methods without
extensive feature engineering.

##  2 Background

We introduce some notation and go over some background material needed for the
model first.

In this work we assume the univariate forecasting setup of DeepAR [48],
however the method also applies to the multivariate setting. Formally we
assume we are given a training dataset of $D\geq 1$ time series
$\mathcal{D}_{\mathrm{train}}=\\{x_{1:T_{i}}^{i}\\}$ where
$i\in\\{1,\ldots,D\\}$ and at each time $t$ we have $x_{t}^{i}\in\mathbb{R}$
or in $\mathbb{N}$. The problem will require us to predict $P\geq 1$ steps
into the future and thus will come with a back-testing test set of these $D$
time series denoted by
$\mathcal{D}_{\mathrm{test}}=\\{x_{T_{i}+1:T_{i}+P}^{i}\\}$. Note that even
though we are denoting the time index $t$ by an incremental counter here, in
reality each $t$ has a date-time associated with it, which increments
regularly based on the frequency of the dataset. Often we will have all the
$T_{i}$ be the same date-time for the time series of a particular dataset we
use.

###  2.1 Probabilistic Forecasting

In the probabilistic forecasting problem, we wish to learn the unknown future
distribution of the data given its past. Rather than considering the whole
history of each time series $i$ in a dataset (which might not be of the same
size) we can instead consider some fixed context window of size $C\geq 1$ of
our choosing and try to learn this potentially complex unknown distribution of
the future values given the context window sized past denoted by

$p_{\mathcal{X}}(x_{T_{i}+1:T_{i}+P}^{i}|x_{T_{i}+1-C:T_{i}}^{i}).$ (1)

Thus, if we denote the parameters of our deep learning model by $\theta$, we
can approximate (1) by an autoregressive model which we can write via the
chain-rule of probability as

$p_{\mathcal{X}}(x_{T_{i}+1:T_{i}+P}^{i}|x_{T_{i}+1-C:T_{i}}^{i};\theta)=\Pi_{t=T_{i}+1}^{T_{i}+P}p_{\mathcal{X}}(x_{t}^{i}|x_{T_{i}+1-C:t-1}^{i};\theta).$

And as we can see at a high level the probabilistic forecasting problem
reduces down to learning some representation of the context window past
together with a distribution model of the next time step(s) given this
representation. For example, the DeepAR model uses a Recurrent Neural Network
(RNN) (like the LSTM [20] or GRU [8]) to represent the context window history
together with associated covariates, denoted by real-valued $\mathbb{R}^{F}$
sized vectors $\mathbf{c}^{i}_{1:T_{i}+P}$ (known for all times), as

$\mathbf{h}_{t}=\mathrm{RNN}(\mathtt{concat}(x^{i}_{t-1},\mathbf{c}^{i}_{t}),\mathbf{h}_{t-1};\theta),$

with $\mathbf{h}_{t}\in\mathbb{R}^{H}$ and $\mathbf{h}_{0}=\vec{0}$, to learn
$p_{\mathcal{X}}(x^{i}_{t}|\mathbf{h}_{t};\theta).$ The distribution head in
DeepAR learns the parameters of some chosen distribution, e.g. the mean and
variance of a Gaussian and we can maximize the log-likelihood of the
distribution with respect to the ground-truth $x^{i}_{t}$ via SGD for all $i$
and $t$ in $\mathcal{D}_{\mathrm{train}}$.

By incorporating different architectural inductive biases to learn the
historic representation together with different distribution heads,
researchers have come up with a zoo of models to solve the probabilistic
forecasting problem in the deep learning domain, for not only the univariate
and multivariate or regular and irregular time series setting, but also for
spatial-temporal problems as well.

###  2.2 VQ-VAE

Leaving time series forecasting aside for the moment we can consider the fact
that latent representations can potentially also be discrete. Thus instead of
learning a continuous valued distributional representation as in Variational
AutoEncoders (VAE) [28], the VQ-VAE model learns a fixed size or discrete set
of representations of the data via an encoder-decoder bottleneck.

Formally the model consists of an encoder that maps its input onto a fixed
sized set of latent variables and a decoder that reconstructs the input from
one of these fixed number of latents. Both the encoder and decoder use this
shared codebook of vectors. Thus if we denote the encoding representation by
$\mathbf{h}^{\mathrm{enc}}\in\mathbb{R}^{E}$, then this vector is quantized
based on its distance to the prototype vectors in the codebook
$\\{\mathbf{z}_{1},\ldots,\mathbf{z}_{J}\\}$ such that the
$\mathbf{h}^{\mathrm{enc}}$ is replaced by the closest prototype vector:

$\mathrm{Quantize}(\mathbf{h}^{\mathrm{enc}}):=\mathbf{z}_{n}\quad\mathrm{where}\quad
n=\mathop{\mathrm{arg\,min}}_{j}\|\mathbf{h}^{\mathrm{enc}}-\mathbf{z}_{j}\|_{2},$
(2)

for a codebook size $J\geq 1$ of our choosing.

The mappings are learned by back-propagating the gradient of a reconstruction
error via the decoder to the encoder using the Straight-Through gradient
estimator [19, 4]. Apart from this loss the VQ-VAE model also has two terms
that encourages the alignment of the vector space of the codebook with the
output of the encoder, namely the codebook loss which pushes the selected
codebook entry $\mathbf{z}$ closer to the encoder representation and the
commitment loss which encourages the output of the encoder to stay close to
the chosen codebook vector to stop it switching too frequently from one code
vector to another. The “stop-gradient” or “detach” operator
$\mathtt{sg}(\cdot)$ blocks the gradients from flowing into its argument and
$\beta$ is the hyperparameter which penalizes the code vector fluctuating,
then the two extra entries of the loss are given by

$\|\mathtt{sg}(\mathbf{h}^{\mathrm{enc}})-\mathbf{z}\|^{2}_{2}+\beta\|\mathtt{sg}(\mathbf{z})-\mathbf{h}^{\mathrm{enc}}\|^{2}_{2}.$
(3)

Since the optimal code would be the average of the representations, [51]
presents an exponential moving average update scheme of the latents instead of
the codebook loss (first term of (3)) during training. The SoundStream [55]
paper proposes to initialize the codebook by the k-means centroids of the
representations from the first batch. Additionally, as proposed in Jukebox
[12], we deploy the heuristic where we replace codebook vectors that have an
exponential moving average cluster size less than some threshold $Q$ by a
random vector from the batch.

##  3 VQ-AR Method

We motivate this method with the observation that time series data in a
dataset visually looks similar with repeated patterns, and thus we ask if
instead of a continuous representation of the histories as learned by most
deep learning models, could a discrete set of representations suffice? With
this in mind, we introduce the VQ-AR model which incorporates an encoder-
decoder RNN together with a Vector Quantizer (VQ) bottleneck.

Unlike in the DeepAR setting, where we learn the temporal representation of
each time series’ context window to predict the distribution of the next step,
we instead first encode it to a quantized representation $\mathbf{z}_{n}$
given by (2) using an encoding-RNN’s state denoted by

$\mathbf{h}_{t}^{\mathrm{enc}}=\mathrm{RNN}_{\mathrm{enc}}(\mathtt{concat}(x^{i}_{t-1},\mathbf{c}^{i}_{t}),\mathbf{h}_{t-1}^{\mathrm{enc}};\theta),$

with $\mathbf{h}_{t}^{\mathrm{enc}}\in\mathbb{R}^{E}$ and
$\mathbf{h}^{\mathrm{enc}}_{0}=\vec{0}$. Subsequently we model the
distribution of the next time step
$p_{\mathcal{X}}(x^{i}_{t}|\mathbf{z}_{n};\theta)$ using a decoding-RNN

$\mathbf{h}_{t}^{\mathrm{dec}}=\mathrm{RNN}_{\mathrm{dec}}(\mathrm{Quantize}(\mathbf{h}^{\mathrm{enc}}_{t}),\mathbf{h}_{t-1}^{\mathrm{dec}};\theta),$

which _only_ takes the Vector-Quantized latent $\mathbf{z}_{n}$ (2) as input,
with $\mathbf{h}_{t}^{\mathrm{dec}}\in\mathbb{R}^{H}$ and
$\mathbf{h}^{\mathrm{dec}}_{0}=\vec{0}$. See Figure 1 for a schematic of the
model.

Figure 1: VQ-AR model schematic at time point $t-1$ for time series $i$.
During training the model takes as input the target $x^{i}_{t-1}$ and
covariates $\mathbf{c}^{i}_{t}$ and outputs the parameters of some chosen
distribution to minimize the negative log-likelihood and Vector Quantizer loss
(5); and during inference (dashed) it allows for sampling of the next time
point $\hat{x}^{i}_{t}$, in an autoregressive fashion.

Because the _prior_ $p(\mathbf{z}_{j})=1/J$ is a uniform categorical
distribution and the _posterior_

$q(\mathbf{z}_{n}|x^{i}_{1:t-1},\mathbf{c}^{i}_{2:t};\theta)=\begin{cases}1&\text{if
}n=\mathop{\mathrm{arg\,min}}_{j}\|\mathbf{h}_{t}^{\mathrm{enc}}-\mathbf{z}_{j}\|_{2}\\\
0&\text{otherwise},\end{cases}$ (4)

is a categorical distribution parametrized by one-hot probability vectors, the
KL-divergence term in the Evidence Lower Bound (ELBO) [6] is constant ($\log
J$), and since the ELBO’s expectation term is over a deterministic
distribution (4), we can simply train the whole model by maximizing the log-
likelihood term of the next time step: $\log
p_{\mathcal{X}}(x^{i}_{t}|\mathbf{z}_{n},\mathbf{h}_{t-1}^{\mathrm{dec}};\theta)$.

Thus, similar to the DeepAR model we minimize the negative log-likelihood of
some chosen parametric distribution and train the whole model end-to-end
together with the VQ loss (3). For example, to model Gaussian emissions the
distribution head consist of a linear layer which outputs the mean
$\mu_{t}\in\mathbb{R}$ and standard deviation $\sigma_{t}\in\mathbb{R}^{+}$,
making sure that the standard deviation is always positive via the
$\mathtt{softplus}(\cdot)$ non-linearity. These parameters are used to
instantiate the distribution $\mathcal{N}(\mu_{t},\sigma_{t})$ from which we
calculate the negative log-likelihood: $-\log\ell_{\mathcal{N}}(x^{i}_{t})$.
The full loss of the model for time point $t-1$ in the training window and
series $i$ is thus denoted by

$\mathcal{L}^{i}_{t-1}(\theta):=-\log\ell_{\mathcal{N}}(x^{i}_{t})+\|\mathtt{sg}(\mathbf{h}^{\mathrm{enc}}_{t})-\mathbf{z}_{n}\|^{2}_{2}+\beta\|\mathtt{sg}(\mathbf{z}_{n})-\mathbf{h}^{\mathrm{enc}}_{t}\|^{2}_{2}.$
(5)

###  3.1 Training

We construct batches $\mathcal{B}$ of context and subsequent prediction window
of total size $T+1$, by sampling a random time series (and corresponding
covariates) from $\mathcal{D}_{\mathrm{train}}$ and then sample this window
randomly within. We will then minimize the mini-batch loss, which is the mean
of all the individual losses (5) denoted by

$\mathcal{L}(\theta):=\frac{1}{|\mathcal{B}|T}\sum_{x^{i}_{1:T+1}\in\mathcal{B}}\sum_{t=1}^{T}\mathcal{L}_{t}^{i}(\theta),$

via the Adam [27] SGD optimizer at each step of training with an appropriate
batch size $|\mathcal{B}|$ and learning rate.

###  3.2 Inference

At inference time we go over the last context sized window of each time series
in the training dataset $\mathcal{D}_{\mathrm{train}}$ and then start
forecasting for the $P$ prediction steps by sequentially sampling from the
distribution head and autoregressively passing this value back into the model
(together with the covariates which are known for all times). In fact, we can
repeat this process $S$ number of times (e.g. $S=100$) to then report any
empirical probability interval of interest as well as probabilistic metrics,
like Continuous Ranked Probability Score (CRPS) [15, 37], with respect to the
ground-truth values in $\mathcal{D}_{\mathrm{test}}$. Point forecasting
metrics can also be evaluated with respect to the empirical median or mean
from the $S$ samples at each time step of the prediction.

Note, unlike in some generative modeling use cases, we _do not_ sample with a
reduced temperature distribution to obtain higher quality samples, nor do we
Beam-search or do Nucleus sampling as in NLP settings during inference.

###  3.3 Covariates

As mentioned in section 2.1, we can incorporate covariates for each time point
of the time series (as well as for all future times) by creating time features
given the granularity of the time series in question. For example, for daily
data we could consider the day of the week, week of the month, month of the
year, etc., features. Although not strictly necessary, but certainly helpful,
we can take inspiration from classical methods and construct time series
covariates, such as lagged features, moving averages, difference to previous
time points, exponential moving averages, and moving percentiles, etc., of
$x^{i}$ (and/or of temporal covariates) depending on the frequency of the time
series in question.

Where applicable holidays or other recurring events can also be turned into
temporal features. A running “age” counter can also serve as an indication of
the length of a time series. Apart from these temporal features, one can embed
the categorical identity $i\in\\{1,\ldots,D\\}$ of each time series into a
vector via a trainable embedding layer whose time independent output is copied
over for all time points being considered. If the time series are grouped then
their categorical tree structure can also be represented via similar
embeddings.

Unless explicitly stated, for all our experiments we only use time, age,
lagged features, and categorical identities as covariates.

###  3.4 Scaling

Entities in time series data can potentially have arbitrary magnitudes and so
to be able to train a shared deep learning model with respect to inputs of any
scale we can utilize the following heuristic: for each context window we can
calculate the mean value of the data in the interval, specifically
$\nu^{i}=\sum_{t=1}^{T}x^{i}_{t}/T$ (if it is not zero or $1$ otherwise) and
divide the time series by this before using it in the model. The outputs from
the distribution heads are then multiplied by this $\nu^{i}$ when sampling or
we rescale in the parameter space of the distribution (e.g. when we require
integer samples from a Poisson or Negative-Binomial head).

We omit stating this heuristic in all our derivations and treat it as an
implementation detail.

##  4 Experiments

We test the performance of VQ-AR for the forecasting task in terms of
performance and also highlight the limitations of the model and generalization
properties of the representations in this section.

For our experiments we use Exchange [31], Solar [31],
Elec.111https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014,
Traffic222https://archive.ics.uci.edu/ml/datasets/PEMS-SF,
Taxi333https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page,
Wikipedia444https://github.com/mbohlkeschneider/gluon-
ts/tree/mv_release/datasets open datasets, preprocessed exactly as in [47], as
well as M5 [35] and proprietary E-Commerce dataset, with their properties
listed in Table 5 in the appendix B. As can be noted in the table, we do not
need to normalize scales for the Traffic dataset. From the names of the
datasets, we see that we cover a number of time series domains including
finance, weather, energy, logistics, page-views, as well as E-commerce.

As indicated above, for evaluation we employ the CRPS metric which measures
the compatibility of a cumulative distribution function (CDF) $F$ with the
ground-truth observation $x$ as

$\mathrm{CRPS}(F,x)=\int_{\mathbb{R}}(F(y)-\mathbb{I}\\{x\leq
y\\})^{2}\,\mathrm{d}y,$

where $\mathbb{I}\\{x\leq y\\}$ is the indicator function which is one if
$x\leq y$ and zero otherwise. CRPS is a _proper scoring function_ , meaning
CRPS attains its minimum when the predictive distribution $F$ and the data
distribution are equal. Employing the empirical CDF of $F$, i.e.
$\hat{F}(y)=\frac{1}{S}\sum_{s=1}^{S}\mathbb{I}\\{x^{(s)}\leq y\\}$ using $S$
samples $x^{(s)}\sim F$ as a natural approximation of the predictive CDF, CRPS
can be directly computed from sampled predictions from the model at each time
point [26]. The final metric is averaged over all the prediction time steps
and time series in a dataset.

###  4.1 Forecasting

We will compare VQ-AR with the following deep learning baseline
_probabilistic_ univariate models

  * •

DeepAR [48]: an RNN based probabilistic model which learns the parameters of
some chosen distribution for the next time point;

  * •

MQCNN [53]: a Convolutional Neural Network model which outputs chosen
quantiles of the forecast upon which we regress the ground truth via Quantile
loss;

  * •

SQF-RNN [14]: an RNN based non-parametric method which models the quantiles
via linear splines and also regresses the Quantile loss;

  * •

IQN-RNN [16]: combines an RNN model with an Implicit Quantile Network (IQN)
[9] head to learn the distribution similar to SQF-RNN;

  * •

LSF [18]: a method that transform a point-estimator, coming from for example
gradient boosting methods, into a probabilistic one;

as well as the classical ETS [23] which is an exponential smoothing method
using weighted averages of past observations with exponentially decaying
weights as the observations get older together with Gaussian _additive_ errors
(E) modeling trend (T) and seasonality (S) effects separately.

We evaluate the model on the datasets detailed above and follow the
recommendations of the M4 competition [34] regarding forecasting performance
metrics. Thus, we also report the mean scale interval score [15]
(MSIS555http://www.unic.ac.cy/test/wp-
content/uploads/sites/2/2018/09/M4-Competitors-Guide.pdf) for a 95% prediction
interval, the 50-th and 90-th quantile percentile loss (QL50 and QL90
respectively), as well as the CRPS score. The point-forecasting performance of
models is measured by the normalized root mean square error (NRMSE), the mean
absolute scaled error (MASE) [24], and the symmetric mean absolute percentage
error (sMAPE) [33]. For pointwise metrics, we use sampled _medians_ with the
exception of NRMSE, where we take the _mean_ over our prediction samples. The
results of our extensive experiments are detailed in Table 1. Note that the
VQ-AR is essentially compressing the time series histories to discrete tokens
which could have applications in edge-computing, by trading off performance
with space.

Table 1: Comparison metrics using different methods: SQF-RNN with 50 knots, DeepAR and VQ-AR with Student-T (-t), Negative Binomial (-nb) or IQN (-iqn) emission heads, ETS, MQCNN, and IQN-RNN on the datasets. Dataset | Method | CRPS | QL50 | QL90 | MSIS | NRMSE | sMAPE | MASE
---|---|---|---|---|---|---|---|---
Exchange | SQF-RNN-50 | 0.010 | 0.013 | 0.006 | 14.15 | 0.020 | 0.013 | 1.800
DeepAR-t | 0.012 | 0.016 | 0.007 | 69.29 | 0.022 | 0.030 | 9.980
ETS | 0.008 | 0.010 | 0.005 | 15.89 | 0.015 | 0.011 | 1.517
IQN-RNN | 0.007 | 0.010 | 0.004 | 17.37 | 0.014 | 0.013 | 3.041
MQCNN | 0.015 | 0.016 | 0.011 | 60.04 | 0.026 | 0.045 | 5.440
VQ-AR-t | 0.010 | 0.013 | 0.007 | 18.10 | 0.019 | 0.015 | 2.658
Solar | SQF-RNN-50 | 0.330 | 0.431 | 0.175 | 5.65 | 0.929 | 1.342 | 1.004
DeepAR-t | 0.418 | 0.543 | 0.254 | 7.33 | 1.072 | 1.393 | 1.275
ETS | 0.646 | 0.661 | 0.383 | 18.55 | 1.112 | 1.546 | 1.938
IQN-RNN | 0.373 | 0.491 | 0.165 | 5.99 | 1.037 | 1.356 | 1.15
MQCNN | 0.928 | 0.960 | 1.535 | 73.58 | 1.920 | 1.838 | 2.248
VQ-AR-iqn | 0.320 | 0.414 | 0.174 | 5.64 | 0.885 | 1.346 | 0.969
Electricity | SQF-RNN-50 | 0.078 | 0.097 | 0.044 | 8.66 | 0.632 | 0.144 | 1.051
DeepAR-t | 0.062 | 0.078 | 0.046 | 6.79 | 0.687 | 0.117 | 0.849
ETS | 0.076 | 0.100 | 0.050 | 9.99 | 0.838 | 0.156 | 1.247
IQN-RNN | 0.060 | 0.074 | 0.040 | 8.74 | 0.543 | 0.138 | 0.897
MQCNN | 0.129 | 0.148 | 0.132 | 30.54 | 1.230 | 0.240 | 2.000
VQ-AR-t | 0.054 | 0.068 | 0.036 | 5.88 | 0.653 | 0.107 | 0.717
Traffic | SQF-RNN-50 | 0.153 | 0.186 | 0.117 | 8.40 | 0.401 | 0.243 | 0.76
DeepAR-t | 0.172 | 0.216 | 0.117 | 8.02 | 0.472 | 0.244 | 0.89
ETS | 0.373 | 0.386 | 0.287 | 17.67 | 0.647 | 0.489 | 1.543
IQN-RNN | 0.139 | 0.168 | 0.117 | 7.11 | 0.433 | 0.171 | 0.656
MQCNN | 1.220 | 0.563 | 2.005 | 116.69 | 0.723 | 0.636 | 2.712
VQ-AR-t | 0.138 | 0.164 | 0.113 | 7.79 | 0.409 | 0.185 | 0.641
Taxi | SQF-RNN-50 | 0.286 | 0.362 | 0.188 | 5.53 | 0.570 | 0.609 | 0.741
DeepAR-nb | 0.299 | 0.379 | 0.203 | 5.44 | 0.610 | 0.582 | 0.771
ETS | 1.059 | 1.297 | 0.617 | 12.24 | 2.147 | 1.159 | 1.552
IQN-RNN | 0.295 | 0.370 | 0.201 | 6.51 | 0.583 | 0.629 | 0.758
MQCNN | 1.262 | 1.451 | 0.488 | 48.61 | 2.645 | 0.912 | 3.041
VQ-AR-nb | 0.286 | 0.362 | 0.193 | 5.43 | 0.572 | 0.570 | 0.741
Wikipedia | SQF-RNN-50 | 0.283 | 0.328 | 0.321 | 23.71 | 2.24 | 0.261 | 1.44
DeepAR-nb | 0.321 | 0.383 | 0.361 | 26.48 | 2.354 | 0.327 | 1.852
DeepAR-t | 0.235 | 0.27 | 0.267 | 23.77 | 2.15 | 0.219 | 1.295
ETS | 0.788 | 0.440 | 0.836 | 61.68 | 3.261 | 0.301 | 2.214
IQN-RNN | 0.221 | 0.254 | 0.251 | 21.78 | 2.102 | 0.193 | 1.214
MQCNN | 0.398 | 0.453 | 0.327 | 38.79 | 2.202 | 0.379 | 2.336
VQ-AR-iqn | 0.231 | 0.266 | 0.252 | 22.09 | 2.106 | 0.208 | 1.261
M5 | SQF-RNN-50 | 0.558 | 0.708 | 0.466 | 8.26 | 1.634 | 1.556 | 0.912
DeepAR-nb | 0.539 | 0.679 | 0.469 | 8.01 | 1.547 | 1.550 | 0.915
ETS | 0.838 | 1.051 | 0.696 | 23.67 | 2.560 | 1.560 | 3.161
IQN-RNN | 0.539 | 0.677 | 0.475 | 8.26 | 1.511 | 1.602 | 0.898
MQCNN | 0.574 | 0.725 | 0.497 | 8.04 | 1.775 | 1.616 | 0.921
VQ-AR-nb | 0.527 | 0.694 | 0.457 | 7.20 | 1.628 | 1.535 | 0.895
E-Commerce | SQF-RNN-50 | 0.545 | 0.627 | 0.648 | 24.82 | 8.492 | 1.653 | 1.246
DeepAR-nb | 0.531 | 0.612 | 0.638 | 38.12 | 8.160 | 1.409 | 1.071
ETS | 2.605 | 2.991 | 2.457 | 117.50 | 31.236 | 1.791 | 10.270
IQN-RNN | 0.534 | 0.627 | 0.590 | 32.92 | 7.588 | 1.650 | 1.076
MQCNN | 0.649 | 0.751 | 0.712 | 27.56 | 9.336 | 1.677 | 1.325
VQ-AR-iqn | 0.483 | 0.562 | 0.555 | 22.33 | 6.702 | 1.629 | 1.074

The M5 competition was a twin competition measuring both point forecasting
accuracy [35] as well as probabilistic uncertainty [36] and the top solutions
were dominated by tree-based methods. In fact, the winning solution to the
accuracy benchmark is based on an ensemble of 220 gradient boosting tree (GBT)
models while the winner of the uncertainty challenge is an ensemble of 126
GBTs with hundreds of hand-crafted features. In Table 2 we provide a
comparison between tree-based methods and neural methods. As detailed in [25]
it is still an open problem to improve deep-learning models with respect to
tree-based methods and as can be seen, the VQ-VR method provides competitive
metrics via a simple _single_ model without any extensive feature engineering.

Table 2: Weighted Quantile Loss (lower is better) for M5 test set predictions using different methods for specified quantiles and their average. Quantile |  M5 winner |  M5 runner-up | DeepAR-t | DeepAR-nb | MQCNN | LSF | VQ-AR-nb
---|---|---|---|---|---|---|---
0.005 | 0.010 | 0.042 | 0.023 | 0.016 | 0.010 | 0.012 | 0.013
0.025 | 0.050 | 0.086 | 0.075 | 0.061 | 0.051 | 0.054 | 0.057
0.165 | 0.312 | 0.337 | 0.319 | 0.319 | 0.320 | 0.316 | 0.316
0.25 | 0.444 | 0.461 | 0.458 | 0.450 | 0.451 | 0.451 | 0.445
0.50 | 0.698 | 0.690 | 0.766 | 0.705 | 0.712 | 0.724 | 0.697
0.75 | 0.687 | 0.722 | 0.821 | 0.699 | 0.709 | 0.728 | 0.690
0.835 | 0.585 | 0.598 | 0.754 | 0.600 | 0.609 | 0.607 | 0.591
0.975 | 0.182 | 0.196 | 0.430 | 0.202 | 0.213 | 0.218 | 0.189
0.995 | 0.055 | 0.076 | 0.310 | 0.074 | 0.084 | 0.086 | 0.063
Mean | 0.336 | 0.357 | 0.440 | 0.347 | 0.351 | 0.355 | 0.340

###  4.2 Robustness

Table 3: Comparison of test set performance of IQN-RNN and VQ-AR on Solar
dataset with different levels of Gaussian noise added to the context window
input at the start of inference. The best metric for each noise level $l$ is
highlighted in bold.

Method | Level | CRPS | MSIS | sMAPE | MASE
---|---|---|---|---|---
IQN-RNN | 0.0 | 0.373 | 5.99 | 1.356 | 1.115
0.2 | 0.397 | 8.22 | 1.358 | 1.182
0.4 | 0.431 | 11.97 | 1.364 | 1.234
0.6 | 0.479 | 17.49 | 1.370 | 1.324
0.8 | 0.536 | 22.71 | 1.382 | 1.453
1.0 | 0.610 | 29.42 | 1.399 | 1.612
VQ-AR-iqn | 0.0 | 0.320 | 5.64 | 1.346 | 0.969
0.2 | 0.388 | 8.95 | 1.362 | 1.144
0.4 | 0.466 | 12.76 | 1.384 | 1.340
0.6 | 0.540 | 15.25 | 1.410 | 1.564
0.8 | 0.623 | 17.66 | 1.429 | 1.801
1.0 | 0.712 | 21.18 | 1.451 | 2.042

We check the robustness of VQ-AR with respect to the noise added to the
context window inputs (but not covariates) when forecasting. Our hypothesis is
that due to the quantization nature of the method, the added noise (up to a
certain level) will correspond to the _same_ codebook vector and thus this
method would be robust to this, as opposed to for example DeepAR. However, too
much noise could also cause the encoder to flip to nearby codebook vectors,
and cause the predictions to be completely off since there is no semantic
meaning in nearby codebook vectors that is explicitly built into VQ-AR.

For some noise level $l=\\{0.2,0.4,\ldots,1.0\\}$ we calculate the standard
deviation of the context window target $\sigma_{\textrm{context}}$ and _add_
to each time point noise sampled from
$\mathcal{N}(0,l\times\sigma_{\textrm{context}})$ and then pass this together
with the covariates to the appropriate RNN.

Table 3 details the metrics for different levels of noise. As we can see with
a small amount of noise the VQ-AR method is more robust than IQN-RNN, however
as we increase the noise the selected codebook vectors become too different
leading to worse performance than with an RNN based method with smooth
representations. The MSIS metric which measures the upper 95% and lower 5%
probability interval however is more robust.

###  4.3 Generalization

The VQ-AR method is essentially an information bottleneck that serves as a
form of regularization. Thus our hypothesis is that the learned cookbook
representations $\\{\mathbf{z}_{j}\\}$ would need to be somewhat universal
time series representations in comparison to methods that can potentially
learn continuous representations without such a constraint. In order to test
this we measure the zero-shot prediction performance of VQ-AR against DeepAR
by first training these two models on the Electricity dataset (without any
categorical time series covariates) and then testing the trained model’s
performance on the test set of Solar and Traffic.

Table 4: Comparison of test set performance of DeepAR and VQ-AR trained on the
Electricity dataset and then tested on Solar and Traffic in a zero-shot
setting. The best metric for each test dataset is given in bold.

Method | Test set | CRPS | MSIS | sMAPE | MASE
---|---|---|---|---|---
DeepAR-t | Solar | 0.522 | 16.60 | 1.441 | 1.432
Traffic | 0.267 | 12.91 | 0.393 | 1.270
VQ-AR-t | Solar | 0.511 | 14.69 | 1.435 | 1.413
Traffic | 0.214 | 10.52 | 0.300 | 0.972

Table 4 highlights the prediction metrics on unseen datasets and as can be
noted, the representations learned as discrete entities tend to perform better
in the zero-shot setting especially since the datasets come from dissimilar
domains.

###  4.4 Ablation

For our ablation study, we wish to test the performance of the model with
respect to the number of codebook vectors $J$. In this experiment for the M5
and E-Commerce datasets we train and test the performance of VQ-AR using
$J=\\{2,4,\ldots,512\\}$. We record the metrics in Figure 2 and to our
surprise, we see that we get decent forecasts using only _two_ codebook
vectors for all the datasets considered (and not only for the ones in the
figure). More codebook vectors tend to give better performance. As
$J\to\infty$ we approach the DeepAR method albeit with more RNN layers.

(a)

(b)

Figure 2: Test set CRPS of VQ-AR-iqn on (a) M5 and (b) E-Commerce datasets
with different number of codebook vectors ($J$). Surprisingly, the decoding
RNN can still give valid predictions given some sequence of just _two_
codebook vectors.

##  5 Related Work

The use of Vector Quantization (VQ) in time series _forecasting_ has not been
explored, as far as we are aware. However, generative models like VQ-VAE have
been used extensively for sequential modeling problems like for videos [43]
and audio/speech [12, 55, 2] generation. Typically these models are trained in
a two-stage process where a VQ-VAE model is first trained and then a
sequential generative model is trained on top of the codebook vectors. The VQ-
Wav2Vec [2] uses a causal convolution to first learn a representation of an
audio signal some steps into the future for further downstream task. This
model is not autoregressive, takes bounded audio inputs only and is not able
to be used as a forecaster which would need to account for the inductive bias
of scale and additional covariates.

The SOM-VAE [13] model is related in that it uses VQ to learn interpretable
discrete representation of sequential data for _clustering_. Also in the paper
[42], the authors discretize the _inputs_ to time series models via different
binning transformations, as a way to mitigate the scaling issue discussed in
section 3.4, and investigate performance given this input transformation on
different architectures.

Neural forecasting methods [5] are related with the caveat that all of them,
as far as we are aware, can be thought of as some component that learns a
representation of the context window in a continuous fashion together with an
emission component, be it point forecasting or probabilistic. CNN [3] based
methods can typically be replaced by RNNs or causal Transformers [52] and most
works try to introduce some time series or problem specific inductive bias
into the architectural choices [32] or emission models [44] or inputs [42]
rather than explicitly on the representations themselves.

##  6 Summary and Discussion

We have presented VQ-AR a conceptually simple and novel yet powerful
autoregressive time series forecasting method based on learning discrete
historic representations of time series histories in order to forecast,
trained in an end-to-end fashion. We have shown not only does it perform
better or similar to the competitive methods considered, but is also able to
achieve good performance in comparison to ensembles of hundreds of tree-based
models which typically dominate forecasting benchmarks, all via a single
model. The architectural principle applied in this method can be adapted for
multivariate forecasting and the underlying RNN based encoder-decoder can be
replaced by a Transformer if so desired.

For future work, we would like to test the inductive bias introduced in this
work for the multivariate setting and also incorporate it within the
Transformer architecture especially for very long time series sequences. As
mentioned in the appendix A, the gradient estimator considered can also be
replaced by its smooth-relaxation variant or by methods with less bias to help
improve the VQ training.

Finally, we have also highlighted a critique of VQ-AR especially in terms of
its robustness to noisy inputs, and in the future, we would like to
investigate how to mitigate this potential issue.

## Acknowledgments

We wish to acknowledge and thank the authors and contributors of the many open
source libraries that were used in this work, in particular: GluonTS [1],
NumPy [17], Pandas [38], Matplotlib [21], Jupyter [29] and PyTorch [40].

## References

  * [1] Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, Lorenzo Stella, Ali Caner Türkmen, and Yuyang Wang.  GluonTS: Probabilistic and Neural Time Series Modeling in Python.  Journal of Machine Learning Research, 21(116):1–6, 2020.
  * [2] Alexei Baevski, Steffen Schneider, and Michael Auli.  vq-wav2vec: Self-supervised learning of discrete speech representations.  In International Conference on Learning Representations, 2020.
  * [3] Shaojie Bai, J. Zico Kolter, and Vladlen Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling, 2018.
  * [4] Yoshua Bengio, Nicholas Léonard, and Aaron C. Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation.  arXiv preprint arXiv:1308.3432, 2013.
  * [5] Konstantinos Benidis, Syama Sundar Rangapuram, Valentin Flunkert, Bernie Wang, Danielle Maddix, Caner Turkmen, Jan Gasthaus, Michael Bohlke-Schneider, David Salinas, Lorenzo Stella, Laurent Callot, and Tim Januschowski.  Neural forecasting: Introduction and literature overview, 2020.
  * [6] Christopher M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2006.
  * [7] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel Ziegler, Jeffrey Wu, Clemens Winter, Chris Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language models are few-shot learners.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 1877–1901. Curran Associates, Inc., 2020.
  * [8] Kyunghyun Cho, Bart van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder–decoder approaches.  In Proceedings of SSST-8, Eighth Workshop on Syntax, Semantics and Structure in Statistical Translation, pages 103–111, Doha, Qatar, October 2014. Association for Computational Linguistics.
  * [9] Will Dabney, Georg Ostrovski, David Silver, and Remi Munos.  Implicit Quantile Networks for Distributional Reinforcement Learning.  In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1096–1105, Stockholmsmässan, Stockholm Sweden, 2018. PMLR.
  * [10] Emmanuel de Bézenac, Syama Sundar Rangapuram, Konstantinos Benidis, Michael Bohlke-Schneider, Richard Kurle, Lorenzo Stella, Hilaf Hasson, Patrick Gallinari, and Tim Januschowski.  Normalizing kalman filters for multivariate time series analysis.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 2995–3007. Curran Associates, Inc., 2020.
  * [11] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota, June 2019. Association for Computational Linguistics.
  * [12] Prafulla Dhariwal, Heewoo Jun, Christine Payne, Jong Wook Kim, Alec Radford, and Ilya Sutskever.  Jukebox: A generative model for music.  CoRR, abs/2005.00341, 2020.
  * [13] Vincent Fortuin, Matthias Hüser, Francesco Locatello, Heiko Strathmann, and Gunnar Rätsch.  Som-vae: Interpretable discrete representation learning on time series.  In ICLR (Poster). OpenReview.net, 2019.
  * [14] Jan Gasthaus, Konstantinos Benidis, Yuyang Wang, Syama Sundar Rangapuram, David Salinas, Valentin Flunkert, and Tim Januschowski.  Probabilistic Forecasting with Spline Quantile Function RNNs.  In Kamalika Chaudhuri and Masashi Sugiyama, editors, Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics AISTATS 2019, Naha, Okinawa, Japan, volume 89 of Proceedings of Machine Learning Research, pages 1901–1910. PMLR, 2019.
  * [15] Tilmann Gneiting and Adrian E Raftery.  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378, 2007.
  * [16] Adèle Gouttes, Kashif Rasul, Mateusz Koren, Johannes Stephan, and Tofigh Naghibi.  Probabilistic Time Series Forecasting with Implicit Quantile Networks.  In Time Series Workshop @ ICML 2021, 2021.
  * [17] Charles R. Harris, K. Jarrod Millman, St’efan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fern’andez del R’ıo, Mark Wiebe, Pearu Peterson, Pierre G’erard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [18] Hilaf Hasson, Bernie Wang, Tim Januschowski, and Jan Gasthaus.  Probabilistic forecasting: A level-set approach.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [19] Geoffrey Hinton, Nitish Srivastava, Kevin Swersky, Tijmen Tieleman, and Abdel-rahman Mohamed.  Neural networks for machine learning: Lecture 15b deep autoencoders.  https://www.cs.toronto.edu/~hinton/coursera_lectures.html, 2012\.
  * [20] S. Hochreiter and J. Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, November 1997.
  * [21] J. D. Hunter.  Matplotlib: A 2D graphics environment.  Computing in Science & Engineering, 9(3):90–95, 2007.
  * [22] R.J. Hyndman and G. Athanasopoulos.  Forecasting: Principles and practice.  OTexts, 2021.
  * [23] Rob J. Hyndman and Yeasmin Khandakar.  Automatic time series forecasting: The forecast package for R.  J. Stat. Soft., 27(3):1–22, 2008.
  * [24] Rob J. Hyndman and Anne B. Koehler.  Another look at measures of forecast accuracy.  International Journal of Forecasting, 22(4):679–688, 2006.
  * [25] Tim Januschowski, Yuyang Wang, Kari Torkkola, Timo Erkkilä, Hilaf Hasson, and Jan Gasthaus.  Forecasting with trees.  International Journal of Forecasting, 2021.
  * [26] Alexander Jordan, Fabian Krüger, and Sebastian Lerch.  Evaluating Probabilistic Forecasts with scoringRules.  Journal of Statistical Software, Articles, 90(12):1–37, 2019.
  * [27] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [28] Diederik P. Kingma and Max Welling.  An Introduction to Variational Autoencoders.  Foundations and Trends in Machine Learning, 12(4):307–392, 2019\.
  * [29] Thomas Kluyver, Benjamin Ragan-Kelley, Fernando Pérez, Brian Granger, Matthias Bussonnier, Jonathan Frederic, Kyle Kelley, Jessica Hamrick, Jason Grout, Sylvain Corlay, Paul Ivanov, Damián Avila, Safia Abdalla, and Carol Willing.  Jupyter notebooks – a publishing format for reproducible computational workflows.  In F. Loizides and B. Schmidt, editors, Positioning and Power in Academic Publishing: Players, Agents and Agendas, pages 87–90. IOS Press, 2016\.
  * [30] Alireza Koochali, Andreas Dengel, and Sheraz Ahmed.  If you like it, gan it—probabilistic multivariate times series forecast with gan.  Engineering Proceedings, 5(1), 2021.
  * [31] Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu.  Modeling Long- and Short-Term Temporal Patterns with Deep Neural Networks.  In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, SIGIR ’18, pages 95–104, New York, NY, USA, 2018. ACM.
  * [32] Bryan Lim, Sercan Ö. Arık, Nicolas Loeff, and Tomas Pfister.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  International Journal of Forecasting, 37(4):1748–1764, 2021.
  * [33] Spyros Makridakis.  Accuracy measures: theoretical and practical concerns.  International Journal of Forecasting, 9(4):527–529, 1993.
  * [34] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos.  The M4 Competition: 100,000 time series and 61 forecasting methods.  International Journal of Forecasting, 36(1):54–74, 2020.  M4 Competition.
  * [35] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos.  The m5 competition: Background, organization, and implementation.  International Journal of Forecasting, 2021.
  * [36] Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos, Zhi Chen, Anil Gaba, Ilia Tsetlin, and Robert L. Winkler.  The m5 uncertainty competition: Results, findings and conclusions.  International Journal of Forecasting, 2021.
  * [37] James E. Matheson and Robert L. Winkler.  Scoring Rules for Continuous Probability Distributions.  Management Science, 22(10):1087–1096, 1976.
  * [38] The Pandas development team.  pandas-dev/pandas: Pandas, February 2020.
  * [39] Youngsuk Park, Danielle Maddix, François-Xavier Aubet, Kelvin Kan, Jan Gasthaus, and Yuyang Wang.  Learning quantile functions without quantile crossing for distribution-free time series forecasting.  In Gustau Camps-Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 8127–8150. PMLR, 28–30 Mar 2022.
  * [40] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8026–8037. Curran Associates, Inc., 2019.
  * [41] Max B. Paulus, Chris J. Maddison, and Andreas Krause.  Rao-blackwellizing the straight-through gumbel-softmax gradient estimator.  In International Conference on Learning Representations, 2021.
  * [42] Stephan Rabanser, Tim Januschowski, Valentin Flunkert, David Salinas, and Jan Gasthaus.  The effectiveness of discretization in forecasting: An empirical study on neural time series models.  In 6th Workshop on Mining and Learning from Time Series, 2020.
  * [43] Ruslan Rakhimov, Denis Volkhonskiy, Alexey Artemov, Denis Zorin, and Evgeny Burnaev.  Latent video transformer.  In Giovanni Maria Farinella, Petia Radeva, José Braz, and Kadi Bouatouch, editors, VISIGRAPP (5: VISAPP), pages 101–112. SCITEPRESS, 2021\.
  * [44] Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski.  Deep state space models for time series forecasting.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 7796–7805. Curran Associates, Inc., 2018.
  * [45] Kashif Rasul, Calvin Seward, Ingmar Schuster, and Roland Vollgraf.  Autoregressive denoising diffusion models for multivariate probabilistic time series forecasting.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 8857–8868. PMLR, 18–24 Jul 2021.
  * [46] Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf.  Multivariate Probabilistic Time Series Forecasting via Conditioned Normalizing Flows.  In International Conference on Learning Representations 2021 (Conference Track), 2021.
  * [47] David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, and Jan Gasthaus.  High-dimensional multivariate forecasting with low-rank Gaussian Copula Processes.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 6824–6834. Curran Associates, Inc., 2019.
  * [48] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  International Journal of Forecasting, 2019.
  * [49] Kamilė Stankevičiūtė, Ahmed Alaa, and Mihaela van der Schaar.  Conformal time-series forecasting.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021.
  * [50] Ruey S. Tsay.  Multivariate Time Series Analysis: With R and Financial Applications.  Wiley Series in Probability and Statistics. Wiley, 2014.
  * [51] Aaron van den Oord, Oriol Vinyals, and Koray Kavukcuoglu.  Neural discrete representation learning.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [52] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is All you Need.  In I. Guyon, U.V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 5998–6008. Curran Associates, Inc., 2017.
  * [53] Ruofeng Wen, Kari Torkkola, and Balakrishnan Narayanaswamy.  A multi-horizon quantile recurrent forecaster.  In Vitaly Kuznetsov, Oren Anava, Scott Yang, and Azadeh Khaleghi, editors, NIPS 2017 Time Series Workshop, 2017.
  * [54] Chen Xu and Yao Xie.  Conformal prediction interval for dynamic time-series.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 11559–11569. PMLR, 18–24 Jul 2021.
  * [55] Neil Zeghidour, Alejandro Luebs, Ahmed Omran, Jan Skoglund, and Marco Tagliasacchi.  Soundstream: An end-to-end neural audio codec.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2021.

##  Appendix A Gradient Estimators

Since the sampling and quantization in (4) are not differentiable the VQ-VAE
paper applies a Straight-Through gradient estimator which propagates the
gradients with respect to the codebook vectors through to the output of the
encoding network. If the codebook vector and encoding representation are close
then this is an adequate approximation but there is no bound on the bias of
the gradients with this scheme and thus it falls under the biased-low variance
taxonomy [41] of gradient estimators as shown in Figure 3.

Figure 3: A taxonomy of different gradient estimator methods, adapted from
[41], ordered chronologically with the method used in this paper marked with a
tick.

##  Appendix B Experiment Details

Table 5: Number of time series, domain, frequency, total training time steps and prediction length properties of the training datasets used in the experiments. Dataset | $D$ | Dom. | Freq. | Time step | Pred. len.
---|---|---|---|---|---
Exchange | $8$ | $\mathbb{R}^{\geq 0}$ | day | $6,071$ | $30$
Solar | $137$ | $\mathbb{R}^{\geq 0}$ | hour | $7,009$ | $24$
Elec. | $320$ | $\mathbb{R}^{\geq 0}$ | hour | $15,782$ | $24$
Traffic | $862$ | $(0,1)$ | hour | $14,036$ | $24$
Taxi | $1,214$ | $\mathbb{N}^{\geq 0}$ | 30-min | $1,488$ | $24$
E-Commerce | $2,731$ | $\mathbb{N}^{\geq 0}$ | day | $887$ | $28$
Wikipedia | $9,535$ | $\mathbb{N}^{\geq 0}$ | day | $762$ | $30$
M5 | $30,490$ | $\mathbb{N}^{\geq 0}$ | day | $1,941$ | $28$

To help reduce confusion we also collect the different hyperparameters in a
single location with their description and default values listed in Table 6.

Table 6: Glossary of the main hyperparameters, their description and default values unless specified. Parameter | Description | Default
---|---|---
$D$ | Number of time series in a dataset |
$P$ | Prediction horizon length |
$C$ | Context window length | $P\times 6$
$F$ | Size of covariate vector |
$E$ |  Encoding RNN’s hidden size and size of codebook vectors | 64
$H$ | Decoding RNN’s hidden size | 40
$J$ | Number of codebook vectors | 128
$\beta$ | Commitment cost | $0.25$
$Q$ | Threshold to replace codebook vector | 2
$S$ | Number of samples | 100

For our experiments we selected the default options for the models as
implemented in the GluonTS [1] library. In particular for VQ-AR, we modified
the base implementation of DeepAR and thus kept the parameters of the models
to be similar where applicable. We trained each model for max_epochs=50 with
the batch_size=256 and a default learning_rate=1e-3. For the evaluation we
used quantiles=[0.1, 0.2, ..., 0.9] to calculate the CRPS metric.


11institutetext: Jagiellonian University, Faculty of mathematics and computer
science, Prof. Stanisława Łojasiewicza 6, 30-348 Kraków, Poland,
contact email: 11email<EMAIL_ADDRESS>
# Exploring shareability networks of probabilistic ride-pooling problems

Michał Bujak 11 Rafał Kucharski 11

###### Abstract

Travellers sharing rides in ride-pooling systems form various kinds of
networks. While the notions of the so-called shareability graphs, has been in
the core of many ride-pooling algorithms, so far they have not been explicitly
analysed. Here, we introduce and examine four kinds of networks resulting from
ride-pooling problems.

We use 147 NYC taxi requests from 2016 pooled into attractive shared-rides
with our utility based ExMAS algorithm and explore resulting shareability
networks. To cover unknown individual properties of pooling travellers, we run
1000 replications of probabilistic pooling process, resulting in richer
representation of weighted graphs instrumental to reveal complex network
structures.

Our findings reveal substantial differences between network structures and
topologies. Properties of ride-pooling networks may be further analysed to
better understand and improve ride-pooling systems.

###### keywords:

ride-pooling, shareability graph, matching graph, structural properties

##  1 Introduction

Ride-pooling, where travellers share a vehicle to reach their destinations, is
becoming an increasingly popular urban mobility alternative. Thanks to
pooling, co-travellers can reduce their travel costs, mobility platforms (like
Uber and Lyft) can more efficiently utilise their fleets and cities can reduce
congestion.

Naturally, travellers sharing their rides form a graph structures, already
identified in the seminal works of Santi et al. (2014) and Alonso-Mora et al.
(2017). While the so-called shareability networks are central to solve the
ride-pooling problems they have not been hitherto explicitly analysed.

Ride-pooling, from the networks perspective, can be seen as linking two or
more nodes. Node is a traveller, or more precisely her/his single trip request
(tuple of origin, destination and request time). The link is formed at the two
levels. First, when all the co-travellers agree to travel together (in our
case this means that the ride is attractive for all the co-travellers).
Second, when they actually travel together (when each traveller is uniquely
assigned to a given ride). This two interpretations of links, yields two
different network structures (second one being a subgraph of the first one).
We will refer to the former representation as _shareability_ , i.e. potential
to travel together, and the latter as _matching_ , i.e. actually travelling
together. Moreover, the links can be formed twofold: first, when travellers
(nodes) are linked directly, and second when we introduce a new kind of node:
the ride, which can be seen as a vehicle trip during which co-travellers
travel together. Linking co-travellers with a ride results in a bipartite
network. We will denote the former as _simple networks_ and latter as
_bipartite networks_. This two levels of twofold distinctions yield four kinds
of networks that we explore in this paper.

Moreover, to represent the probabilistic nature of underlying pooling
decisions made by travellers, we propose the probabilistic version of pooling
algorithm. By adding the random noise to the utility formulas of ride-pooling
algorithm we introduce the behavioural non-determinism to the pooling problem.
When replicated multiple times we obtain rich and meaningful structures of
weighted graphs which well estimate probabilistic graphs. With those we can
observe both the stable matches with strong links, as well as components where
pooling compositions can vary substantially between replications.

Such weighted shareability graphs, consistent with the actual travel
behaviour, can be useful not only for a deeper understand of the still
challenging ride-pooling problems, but also to identify communities in ride-
pooling networks, improve attractiveness of pooling services or control the
virus spreading in ride-pooling networks (like in Kucharski et al. (2021)).

###  1.1 Background

The potential of ride-pooling to improve urban mobility has triggered a
bursting stream of research (Agatz et al. (2011); Chan and Shaheen (2012);
Furuhata et al. (2013); Tachet et al. (2017); Bischoff et al. (2017)), which
led to an algorithmic breakthrough in the seminal works of Santi et al. (2014)
and Alonso-Mora et al. (2017). First, in Santi et al. (2014), the so-called
shareability network was introduced along with a methodology to match
travellers in pairs. Further exploited in Alonso-Mora et al. (2017) with the
complete algorithm to efficiently match incoming requests with available
vehicles by sequentially adding new co-travellers into vehicles with empty
seats and adjusting their routes. By applying cutoffs on maximal detours and
delays the pooled rides are acceptable for travellers and by minimising the
service costs (vehicle hours) in the matching the solution becomes optimal for
the supply side (platform and/or drivers). In parallel, behavioural research
has made advancements in better understanding users’ so-called willingness to
share (Kang et al. (2021); Lavieri and Bhat (2019); Alonso-González et al.
(2021); Lazarus et al. (2021)). Such, utility-driven approaches, not only put
the traveller in the centre but also allow to incorporate inherent non-
determinism of travellers’ decision processes.

Hitherto shareability and matching structures have been only introduced. We
formalise and analyse different methods of presenting their topologies.
Furthermore, we introduce random noise, hence the novel non-deterministic
approach to ride-pooling. Based on replications, we introduce weighted graphs
which estimate the probabilistic graphs. We analyse their evolution and
stability.

##  2 Methodology

###  2.1 Ride-pooling algorithm

ExMAS (Kucharski and Cats, 2020) is an offline algorithm which matches the
travellers into a shared ride if it is more attractive than travelling alone.
More precisely, travellers are assigned to the shared rides if individual
utility of the shared ride exceeds a non-shared ride utility for each of the
co-travellers. Value of the utility depends on the individual- and system-
dependent variables with the following formula, given as the difference
between shared ride and non-shared ride:

$U=U^{s}-U^{ns}=\beta^{c}\lambda
l+\beta^{t}(t-\beta^{s}(\hat{t}+\beta^{d}\hat{t}^{p}))+\epsilon>0,$ (1)

where $U$, $U^{s}$, $U^{ns}$ denote respectively utility gain due to sharing,
utility of shared ride and utility of non-shared ride. $\lambda$ stands for
discount for sharing a ride and is controlled by the system operator.
$\beta^{c}$, $\beta^{t}$, $\beta^{s}$ $\beta^{d}$ are the behavioural
parameters: cost sensitivity, value of time, sharing discomfort and delay
sensitivity, respectively. $t$ and $\hat{t}$ stand for travel time of non-
shared and shared ride, respectively, $\hat{t}^{p}$ is a delay associated with
pooling and $\epsilon$ is a random term.

We say that a group of travellers can share a ride, if for each of them the
individually perceived utility of the ride is positive. Those rides contribute
to the set of attractive (feasible) rides. ExMAS identifies such rides
hierarchically, first exploring pairwise shareable trips, which are then
extended to triples, quadruples, etc. until no cliques can be found in
shareability networks. On the set of identified attractive shared rides we
solve the assignment problem in which travellers are uniquely assigned to a
rides, formulated as a matching on a bipartite graph aimed to minimise the
total distance travelled. From now on, we will refer to the set of feasible,
attractive rides identified with ExMAS as _shareability_ and the optimal
solution as a _matching_. Full details of ExMAS algorithm are available in
Kucharski and Cats (2020).

Hitherto, the deterministic version of ExMAS has been used, with
$\epsilon=0$.Here, we decided to adopt a probabilistic model in which $U$
follows a normal distribution (probit model). Formally, probability that
traveller $i$ wants to participate in ride $j$ is described as follows,
$\mathbf{P}(t_{i}\in r_{j})=\Phi_{0,\sigma}(U_{s}-U_{ns})$, where
$\Phi_{0,\sigma}$ stands for CDF of the normal distribution with mean $0$ and
standard deviation $\sigma$. For the sake of simplicity, we omit in $U$
additional indexes, however keeping in mind that $U$, $U_{ns}$ and $U_{s}$ are
individually calculated for each traveller for each ride. $\Phi_{0,\sigma}$ is
a strictly increasing bijection, hence, we can write:

$\displaystyle\mathbf{P}(t_{i}\in r_{j})=\Phi_{0,\sigma}(U_{s}-U_{ns})>v$
$\displaystyle U_{s}-U_{ns}>\Phi_{0,\sigma}^{-1}(v).$

Aforementioned inequalities yield a simple way to draw values to estimate the
distributions of non-deterministic ride-pooling. Each replication yields
simulated values of the random term, hence giving explicit values of utility
per each set of travellers.

###  2.2 Construction of networks

Both finding the shareability and the matching can be represented with
networks. The most straightforward representation is a network, where
travellers constitute nodes, linked when they share a ride (row one in tab.
1). In such _shareability_ network (first column in tab.1), travellers are
linked if they find it attractive, whereas in the _matching_ network (second
column) if the ride is part of the optimal solution. While the latter is the
actual realisation of ride-pooling, i.e. linked travellers actually travel
together, the former reflects the possibility to travel together (which may
not be realised due to other constrains).

Such network representation is, however, incomplete. From the shareability
graph, we cannot conclude if, for example, a triangle means that the three
travellers can share a ride together or whether they are just pairwise
shareable (Figure 1(a)). In order to avoid this information loss, an
alternative representation if often proposed. Both a shareability and matching
network can be formulated as bipartite graphs, where nodes are not only
travellers but also the rides (vehicles). Traveller is connected to a ride in
the bipartite shareability network if for her/him the ride is attractive
(Figure 1(c)), in case of bipartite matching network travellers are linked
only to rides being part of the optimal solution (Figure 1(d)). This leads to
the four distinct network structures, as summarised in the Table 1 and
illustrated with Figure 1.

network | nodes | links
---|---|---
shareability | travellers | attractive pooling (between travellers)
matching |  | optimal matching (between travellers)
bipartite shareability | travellers | attractive pooling (between travellers and rides)
bipartite matching | and rides | optimal matching (between travellers and rides)
Table 1:  Four networks of ride pooling.

(a) Shareability network

(b) Matching network

(c) Bipartite shareability

(d) Bipartite matching

Figure 1: Four representations of networks in ride-pooling on a single
replication from the case of 147 pooling travellers in NYC. The shareability
network (a) is highly connected, many nodes (travellers) are linked with
multiple other potential co-travellers. However, when travellers are uniquely
assigned to vehicles in the matching network (b) the structure becomes
dramatically more sparse. The bipartite forms of those two networks have
similar structures. Travellers (top) are linked with rides (bottom). Among
multiple potentially feasible links (c) only few remain the part of actual
matching (d).

##  3 Experiment

We experiment with trips actually requested within half an hour in New York
City on January 2016. Each of 147 requests links origin with destination at a
given time (as illustrated on Figure 2). We assumed the fare $\beta^{c}$ is
$1.5$ $/km, sharing discount $\lambda$ is $30\%$, value of time
$\beta^{t}=9$$/h and willingness to share $\beta_{s}=1.3$ (most of which based
either on empirical observations or actual system values). We set pooling
parameters to be constant among all the replications (deterministic part) and
assume the random term $\epsilon$ of utility $U$ to follow standard normal
distribution.

Figure 2: Demand dataset for experiments: 147 trip requests from Jan 2016 in
Manhattan. Green dots are origins, orange destinations.

It is important to note that even though we can explicitly calculate utility
of a ride for any group of travellers hence probability of corresponding
links, there is no easy way to determine probability in the matching. Thus, to
observe different realisations of the random term $\epsilon$, we have
performed $1000$ replications. For each, we have constructed networks such as
shown in the Figure 1.

###  3.1 Results

First, we look at the properties of ride-pooling networks and their stability
across the replications. Aggregated results are presented in the Table 2. Some
of the networks’ properties remain constant by their definition. For example,
the number of nodes in the simple networks is equal to the number of
travellers. Similarly, the number of links in bipartite matching. However,
investigating different structural characteristics reveals interesting
patterns. Relatively big standard deviation in the number of nodes and edges
in the bipartite shareability suggests that there is indeed a significant
structural variance introduced by the random noise. Low deviation in the
bipartite shareability network of the average degree of a ride suggests that
introduced noise has a subsidiary impact on the average occupancy of vehicles
among feasible rides. For the same type of network, average degree of
traveller nodes has a significant deviation. The fact implies that there is a
substantial group of rides which balance on the verge of attractiveness.
Hence, minor disturbance introduced by the noise can have a significant impact
on the number of feasible rides for a traveller. Furthermore, relating to the
variance of average degree of nodes in the simple shareability network, we may
conclude that noise is pronounced both in potential sharing between two
travellers and for the higher degree rides. Relative size of the greatest
component shows significant variance for the bipartite shareability network.
Considering the fact that the maximal value is $1$, while the average is
$0.96$, we can conclude that often the greatest component becomes much
smaller. That may suggest forming at least two big components which are not
always connected. Further investigation could prove communities presence
within the components.

| nodes | edges | degree | deg.rides | deg.travellers | greatest component
---|---|---|---|---|---|---
shareability | 147 (0) | 470 (47) | 6,394 (0,638) |  |  | 0,853 (0,027)
matching | 147 (0) | 61 (6) | 0,831 (0,082) |  |  | 0,030 (0,004)
bipartite shareability | 1164 (166) | 2306 (497) | 3,922 (0,287) | 2,249 (0,121) | 15,692 (3,384) | 0,961 (0,11)
bipartite matching | 250 (3) | 147 (0) | 1,175 (0,15) | 1 (0) | 1,425 (0,043) | 0,022 (0,002)

Table 2: General properties and their variability in $1000$ replications (std
in brackets) for four introduced graph structures.

We can now introduce a weighted graphs which aggregate realisations of the
probabilistic graph. For the simple networks, nodes remain constant. Hence, we
link two nodes if for any of the replications they were connected. The weight
of the link is exactly the number of occurrences of the link in the
experiment. We can introduce bipartite networks in nearly the same way. The
only difference is that for rides nodes, we must take union of all rides which
occurred in any replication. The introduced graph structures are presented in
the Figure 3. Width of the link proportionally depends on the weight.

(a) Shareability network

(b) Matching network

(c) Bipartite shareability

(d) Bipartite matching

Figure 3: Weighted shareability networks. After several replications of the
stochastic matching process the network wiring evolves. The potential
shareability network (a) is becoming even more connected, similarly to the
matching, where some travellers are stably matched with high probability,
while others build a broad set of potential co-travellers (b). Similar trends
can be seen in the bipartite versions of the networks (c and d).

Intrigued by the structure of graphs in the Figure 3, we investigated the
evolution of bipartite shareability networks. Figure 4 presents weighted
graphs constructed after a certain part of the simulation. We have found that
even in the latter replications, new ride nodes appeared. Those nodes usually
correspond to high order rides, where travellers are matched thanks to
positive realisations of the random noise. Considering the fact that the
solution is optimal with respect to saved vehicle hours, attractive high order
rides are likely to be assigned in the final matching. Nonetheless, we can see
that there are certain rides which remain attractive and assigned in nearly
every replication.

(a) 1st step

(b) 5 rep.

(c) 10 rep.

(d) 100 rep.

(e) 900 rep.

Figure 4: Evolution of the bipartite matching between travellers (left) and
rides (right) (panel d on figs. 1 and 3). The number of rides (on the right)
as well as number of connections increases. Some highly attractive relations
become stronger, while others appear rarely, with low probability.

To gain more understanding of evolution in the Figure 4, we investigated
properties of the (weakly) growing weighted matching graph in the non-
bipartite form. Figure 5 resents how the average node degree and relative size
of the greatest component grow. Those properties have a direct, physical
interpretation and application. Average degree (blue) of the node reflects the
average number of pairwise distinct co-travellers that a traveller meets
within the $1000$ realisations. Size of the greatest component (red) gives an
upper boundary for spreading in the network starting in one point (to obtain
exact results we should translate our partial results to temporal graph).

Figure 5:  Evolution of graph properties along the replications. As we explore
the probabilistic nature of ride-pooling, the graph properties evolve. The
greatest component of the matching graph (network B in the Figure 3) quickly
reaches stable value and only slightly evolves after 900 days of evolution.
Average degree also increases, yet slower, gradually increasing until the end
of the simulation.

##  4 Conclusions

We have presented four ways to represent the networks resulting from ride-
pooling problems. In the first type of the networks, each node corresponds to
a traveller. Second representation incorporates also rides as nodes, enabling
an exhaustive representation for the feasibility network.

Considering individual preferences not captured by a deterministic model, we
introduce a random noise assigned to each traveller. As a result, within each
replication we got different realisations of the random variables. We analyse
the introduced variability in the context of structural properties of
networks. Weighted graphs representing the estimated probability of occurrence
of a link yield new form of networks in ride-pooling. Tracking evolution of
the graphs gave us a better understating of the stability of the underlying
structures. In future works, we plan to focus on the communities and analysing
spreading patterns in the network.

#### Acknowledgements

This research was funded by National Science Centre in Poland program OPUS 19
(Grant Number 2020/37/B/HS4/01847).

## References

  * Agatz et al. (2011) Niels Agatz, Alan L Erera, Martin WP Savelsbergh, and Xing Wang.  Dynamic ride-sharing: A simulation study in metro atlanta.  _Procedia-Social and Behavioral Sciences_ , 17:532–550, 2011.
  * Alonso-González et al. (2021) María J Alonso-González, Oded Cats, Niels van Oort, Sascha Hoogendoorn-Lanser, and Serge Hoogendoorn.  What are the determinants of the willingness to share rides in pooled on-demand services?  _Transportation_ , 48(4):1733–1765, 2021.
  * Alonso-Mora et al. (2017) Javier Alonso-Mora, Samitha Samaranayake, Alex Wallar, Emilio Frazzoli, and Daniela Rus.  On-demand high-capacity ride-sharing via dynamic trip-vehicle assignment.  _Proceedings of the National Academy of Sciences_ , 114(3):462–467, 2017.
  * Bischoff et al. (2017) Joschka Bischoff, Michal Maciejewski, and Kai Nagel.  City-wide shared taxis: A simulation study in berlin.  In _2017 IEEE 20th international conference on intelligent transportation systems (ITSC)_ , pages 275–280. IEEE, 2017.
  * Chan and Shaheen (2012) Nelson D Chan and Susan A Shaheen.  Ridesharing in north america: Past, present, and future.  _Transport reviews_ , 32(1):93–112, 2012.
  * Furuhata et al. (2013) Masabumi Furuhata, Maged Dessouky, Fernando Ordóñez, Marc-Etienne Brunet, Xiaoqing Wang, and Sven Koenig.  Ridesharing: The state-of-the-art and future directions.  _Transportation Research Part B: Methodological_ , 57:28–46, 2013.
  * Kang et al. (2021) Shuqing Kang, Aupal Mondal, Aarti C Bhat, and Chandra R Bhat.  Pooled versus private ride-hailing: A joint revealed and stated preference analysis recognizing psycho-social factors.  _Transportation Research Part C: Emerging Technologies_ , 124:102906, 2021.
  * Kucharski and Cats (2020) Rafał Kucharski and Oded Cats.  Exact matching of attractive shared rides (exmas) for system-wide strategic evaluations.  _Transportation Research Part B: Methodological_ , 139:285–310, 2020.
  * Kucharski et al. (2021) Rafał Kucharski, Oded Cats, and Julian Sienkiewicz.  Modelling virus spreading in ride-pooling networks.  _Scientific Reports_ , 11:1–11, 2021.
  * Lavieri and Bhat (2019) Patrícia S Lavieri and Chandra R Bhat.  Modeling individuals’ willingness to share trips with strangers in an autonomous vehicle future.  _Transportation research part A: policy and practice_ , 124:242–261, 2019.
  * Lazarus et al. (2021) Jessica R Lazarus, Juan D Caicedo, Alexandre M Bayen, and Susan A Shaheen.  To pool or not to pool? understanding opportunities, challenges, and equity considerations to expanding the market for pooling.  _Transportation Research Part A: Policy and Practice_ , 148:199–222, 2021.
  * Santi et al. (2014) Paolo Santi, Giovanni Resta, Michael Szell, Stanislav Sobolevsky, Steven H Strogatz, and Carlo Ratti.  Quantifying the benefits of vehicle pooling with shareability networks.  _Proceedings of the National Academy of Sciences_ , 111(37):13290–13294, 2014.
  * Tachet et al. (2017) R Tachet, O Sagarra, P Santi, G Resta, M Szell, S H Strogatz, and C Ratti.  Scaling Law of Urban Ride Sharing.  _Nature Publishing Group_ , pages 1–6, 2017.  10.1038/srep42868.  URL http://dx.doi.org/10.1038/srep42868.