# An overview of time series point and interval forecasting based on similarity of trajectories, with an experimental study on traffic flow forecasting

İlker Arslan, Can Hakan Dağıdır, Ümit Işlak

<figure><img src="image_1.png" /></figure>

To cite this version:

İlker Arslan, Can Hakan Dağıdır, Ümit Işlak. An overview of time series point and interval forecasting based on similarity of trajectories, with an experimental study on traffic flow forecasting. 2024. hal-04670731

HAL Id: hal-04670731

https://hal.science/hal-04670731v1

Preprint submitted on 12 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An overview of time series point and interval forecasting based on similarity of trajectories, with an experimental study on traffic flow forecasting

İlker Arslan*

Can Hakan Dağıdır†

Ümit Işlak‡

August 12, 2024

## Abstract

The purpose of this paper is to give an overview of the time series forecasting problem based on similarity of trajectories. Various methodologies are introduced and studied, and detailed discussions on hyperparameter optimization, outlier handling and distance measures are provided. The suggested new approaches involve variations in both the selection of similar trajectories and assembling the candidate forecasts. After forming a general framework, an experimental study is conducted to compare the methods that use similar trajectories along with some other standard models (such as ARIMA and Random Forest) from the literature. Lastly, the forecasting setting is extended to interval forecasts, and the prediction intervals resulting from the similar trajectories approach are compared with the existing models from the literature, such as historical simulation and quantile regression. Throughout the paper, the experimentations and comparisons are conducted via the time series of traffic flow from the California PEMS dataset.

## 1 Introduction

A time series is a list of data points that are indexed by time. Examples of time series include the stock price of a given company, hourly electricity demand, retail sales of a product and the daily temperature at a given city. The problems related to time series forecasting has become of extreme importance in recent decades due to their wide applications. We refer to the texts [12] and [15] for two general references on time series analysis. The purpose of this manuscript is to give an overview of a certain time series forecasting methodology in detail, namely, the use of similar trajectories.

The experiments below will be based on the well studied problem of traffic flow forecasting. The data set we use is the traffic flow data of California PEMS data, and the relevant details are explained in Section 2.1. Our forecasting approach follows the lines of [17] in which they use similar trajectories to obtain a point forecast for the next time step. The method in cited work can be summarized as follows. In order to provide a point (or, interval) forecast for time $T + 1$ when you are at time $T \in \mathbb{N}$, first look for (windows of) trajectories in the past that are similar to the recent observations up to time $T$. Then after obtaining some candidates from these similar trajectories, ensemble them via an appropriate function to reach at a final forecast. Since the process requires the selection and assembly of nearest neighbors, the method is sometimes known to be the K-NN, K referring to the number of similar trajectories that are used in obtaining the forecast. Focusing on the traffic flow case, note that one has seasonal patterns such as daily, weekly and monthly

<footnote>*MEF University, Mechanical Engineering Dept., Istanbul, Turkey. email: arslanil@mef.edu.tr</footnote>

<footnote>†Boğaziçi University, Mathematics Dept., Istanbul, Turkey. email: can.dagidir1@boun.edu.tr</footnote>

<footnote>‡Boğaziçi University, Mathematics Dept., Istanbul, Turkey. email: umit.islak1@boun.edu.tr</footnote>
seasonality, and these need to be taken into account in selection of the nearest neighbors, and
relevant discussions will be included below.

We refer to the papers [17], [40], [44] for the use of similarity of trajectories in traffic flow. Among
these, as noted earlier, [17] should be emphasized in our case since our setup will primarily use the
same reasoning behind it. Of course the idea here is quite natural and has emerged in various other
fields other than the traffic flow. For a general discussion of the topic centered around economics,
but also applicable to other fields, see [8]. The papers [9] and [10] are again on similarity based
approaches for predictions, with more technical notes included. Trajectory similarity is used in
various other fields as well, see [38] for prediction of remaining life time, [24] for location prediction
and [47] for predicting daily tourist arrival as examples. The literature is too vast to list here, but
let us note that most of the work in terms of time forecasting is within the title K-NN.

Although our main interest here is not the particular case of traffic flow forecasting, we mention
some pointers to the literature in this area too for the sake of completeness. One classical approach
in traffic flow forecasting is the use of econometric models, such as ARIMA and its variations. See,
for example [5], [22] and [23], for relevant work. Both machine learning models and deep learning
models are extensively used lately. Let us point out the papers [23] on SVR, [19] on ANN and
Gradient Boosting on [37], as some exemplary work. Also note that these can be combined with
ARIMA and related ones in an additive decomposition manner ([25], [45]). Regarding the deep
learning models, we mention a few work based on LSTM: [20], [35], and [43]. Lastly note that [3],
[23] and [27] provide comparisons on different aspects of traffic flow forecasting.

Besides the point forecast problem in time series, we will also be interested in interval forecasts
below. Although these intervals are not studied in the traffic flow case, they are widely used in
several fields in order to supplement the point forecasts. As an example, electricity price forecasting
is a popular topic in which prediction intervals are widely used, and obtaining such intervals in this
case is quite challenging due to the spiky behavior of the corresponding time series. See [28] and
[29] for relevant discussions. There are various methods for obtaining prediction intervals in the
literature, for example, the historical method, bootstrapping and distributional prediction intervals
[29]. To the best of our knowledge, the attempts using a direct use of similar trajectories in this
manuscript are new, and we were not able to find such references neither in the field of traffic flow
nor in other ones.

Overall, our experiments showed that the point forecasting methods based on similarity of tra-
jectories are competitive with various models used in the literature. We will provide the comparison
results between AR, ARIMA, Random Forest and a basic version of the similarity method in Section
3.2. Our experimentations actually included various other models such as gradient boosting, ANN
and SVR, but we will not be include a comprehensive list, since our main purpose is to experiment,
discuss and compare the methodologies within the similarity of trajectories framework. As will be
clear from the upcoming sections, the similarity approach has an extensive flexibility in it, and
the theoretical classification of problems that can be handled with a certain given method can be
interesting. Due to the presence of several options in design of the algorithms here, the computa-
tional cost may or may not be satisfying. The proposed methodologies in general are grounded in
the premise that historical patterns sharing similar characteristics can provide valuable insights for
estimating future trajectories. Similarity models benefit from the observations from both distant
history and outlier (or minority) segments of historical data. With this ability, similarity models
can utilize the diversity in the historical data without sacrificing the overall accuracy. In addi-
tion to the accuracy and robustness of our forecasts, analyzing similarity models may enhance our
understanding of the underlying system dynamics.

The contributions of the current work can be summarized as follows:

1. The basics of point and interval forecasting based on similarity of trajectories are explained in a simple setup in detail in Section 3. We hope that this exemplary case gives a clear
description of the method to a general audience.

2. Various point forecast mechanisms based on similarity of trajectories are discussed. In particular, the multi-hourly model and the local regression, that are to be described at length, approaches provide the best results among our experiments. The designs are flexible by nature and seem to be open to further improvement. Also the forecasts obtained via the direct use of similar trajectories help in terms of the interpretability.

3. Results corresponding to various distances are collected and listed together, and are compared with the methodologies suggested in the current manuscript. Some models outside of our approach are also included to make the position of the present attempts clear.

4. Experimentations for multi-step forecasts via the use of similarity are done, and the results seem to be promising, especially under the setting of seasonal similarity.

5. Similarity of trajectories approach is used to obtain prediction intervals in time series framework, complementing the point forecasts. Comparisons for this case with other standard methods are included as well. A new algorithm for obtaining prediction intervals is introduced which is a mixture of the so-called historical simulation and the similarity of trajectories methods.

The rest of the paper is organized as follows. The next section is on a description of the traffic flow data set we use, and it further contains the relevant background on the distances of interest, the forecast evaluation metrics, and sample quantiles. In Section 3, we stop for a case analysis based on the weighted $ℓ^2$ distance. Section 4 is devoted to discussions on selection and handling of the candidate trajectories that are to be used for forecasting purposes. The point forecast methods handled in this paper are listed in Section 5. The results of our experiments for point forecasts are also given in this part. Section 6 includes the approaches for interval forecasts based on similarity, and also provides the experimentation results for the prediction intervals. We conclude the paper with a short discussion in Section 7.

# 2 Background

## 2.1 Problem Description and Data set

Given the first $t$ instances $X_1, X_2, \dots, X_{t-1}, X_t$ of a time series $X$, the general problem we are interested in is to provide a forecast for $X_{t+1}$. Using similarity of trajectories, this problem can be approached as follows. First, the time series data is transformed into a trajectory-based representation. (Trajectories can be thought as finite subsequences corresponding to consecutive indices of sequence of time series. The trajectories to be considered here will be of the same length.) Secondly, similar trajectories are identified using $K$-NN. Then, the succeeding observations of similar trajectories, which are to be called *candidates*, are used for forecasting.

In our experiments, we use real world traffic data from California Freeway Performance Measure Systems (PeMS). In particular, the experiments are done with the flow data for Station 312757. The data contains five-minute intervals of observations such as total traffic flow and average speed. The use of 15-minute long intervals is advocated in the literature since using 15 minutes instead of 5 smoothens out the noise while keeping the data meaningful [23]. In our case, we achieve this by combining 3-tuples of non-overlapping five-minute intervals by taking their sum.
<figure><img src="image_3.png" /><figcaption>Figure 1: Weekly traffic flow data</figcaption></figure>

During this study we use the data of 20 months duration from October 2020 to June 2022. The data split is designed to contain at least a year long observations in each split to capture seasonality.

We split the data into four parts: tune reference, tune query, test reference and test query. Tune splits are used for optimizing the hyperparameters and test splits are used for measuring the performance. Tune split, containing tune reference and tune query, ranges from October 2020 to February 2022 and test split ranges from February 2021 to June 2022. Due to the nature of similarity based models, query and reference splits are introduced. Query split contains the observations that are desired to be estimated. For each observation in the query split, a reference set, consisting of past observations, is defined. Reference splits contain past observations ranging along at least 12 and at most 16 months long time interval. Tune query split ranges from October 5, 2021 to February 5, 2022 and for each observation in the query split, a reference split is defined to range from October 5, 2020 to the date (and time) of that observation. Test query split ranges from February 2022 to June 2022 with reference split starting from February 2021, again to the date (and time) of that observation. Here we emphasize that the size of reference sets differ for each observation. To be more concrete, for two observations at November 10, 2021, 10:30 AM and January 3, 2022, 6:00 AM, two different reference sets are defined. The former reference set ranges from October 5, 2020 to November 10, 2021, 10:30 AM and the latter reference set again starts from October 5, 2020 but ends in January 3, 2022, 6:00 AM.

The following figure demonstrates the data split done in our study. Note that the significant change in traffic flow around October 2020 is due to the relaxation of Covid-19 measures.
<figure><img src="image_4.png" /><figcaption>Figure 2: Data split</figcaption></figure>

The missing data in our study is handled in two different ways. First, when the length of consecutive missing data is less than one hour, the corresponding flow data from previous week is substituted. Secondly, for longer consecutive periods (which actually only occurs in February 2022) we use the average of the same days and hours of previous three weeks.

## 2.2 Distances

In this subsection we define certain distance functions that can be used for measuring the distance between two trajectories. Before moving further let us provide some references that contain more details about the distances discussed here, and several other ones. [32] compares various distances, and discuss their use in the clustering framework. [34] gives an overview of distances such as Frechet distance, the dynamic time warping, longest common subsequences. The paper contains four experimental studies for comparison purposes. [33] and [36] are two other references that discuss various aspects of measures for similarity of trajectories.

Now, let us provide the definitions of some common distances between trajectories, and discuss them briefly. Below $\mathbf{x} = (x_1, \dots, x_L)$ and $\mathbf{y} = (y_1, \dots, y_L)$ are considered to be $L \in \mathbb{N}$ instances of a real valued time series.

1. $\ell^p$ distance and related. For $p \ge 1$, $p$-distance between $\mathbf{x}$ and $\mathbf{y}$ is defined to be

$$d_p(\mathbf{x}, \mathbf{y}) = \left( \sum_{i=1}^{L} |x_i - y_i|^p \right)^{1/p}.$$

$\|\mathbf{x} - \mathbf{y}\|_p$ is sometimes written for $d_p(\mathbf{x}, \mathbf{y})$. In this notation the norm of $\mathbf{x}$ is $\|\mathbf{x}\|_p$.

Two special cases should be noted. When $p = 2$, $\ell^p$ distance is known as the **Euclidean distance** between $\mathbf{x}$ and $\mathbf{y}$, perhaps the most well-known distance in the literature. Some authors
also consider the **squared Euclidean distance**

$$d_2^2 = \sum_{i=1}^{L} |x_i - y_i|^2$$

in order to avoid dealing with the square roots, but that will not be an issue for us.

The second special case of the $l^p$ distance we shall mention is the **Manhattan distance** corresponding to $p = 1$:

$$d_M(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{L} |x_i - y_i|.$$

Manhattan distance is also known as $l^1$ distance and the taxicab distance. Computations involving Manhattan distance are faster than the Euclidean distance due to absence of taking square roots.

We lastly note the weighted version of the Euclidean distance. Given weights $w_i \in (0, 1)$ adding up to 1, the corresponding **weighted Euclidean distance** is defined by

$$d_w(\mathbf{x}, \mathbf{y}) = \left( \sum_{i=1}^{L} w_i |x_i - y_i|^2 \right)^{1/2}.$$

Our preliminary analysis of the similarity based forecasting methods in Section 3 will be using the weighted Euclidean distance, where the weights will be selected in a linear way according to the recentness of observations [17]. In present notations, considering $\mathbf{x}$ and $\mathbf{y}$ as the compared trajectories, these weights can be written as $w_i = \frac{i}{L(L+1)/2}$.

**2. sup distance** is defined by

$$d_{\infty}(\mathbf{x}, \mathbf{y}) = \max_{i=1,...,L} |x_i - y_i|.$$

Other names for the sup distance include the maximum distance and Chebyshev distance. Since the sup distance looks for uniformly similar patterns, it can be difficult to find similar trajectories in our setting below.

In case one is interested in only the similarity of certain coordinates of $\mathbf{x}$, $\mathbf{y}$, one may consider the maximum in previous definition over a subset $S$ of $\{1, ..., L\}$. Below we will focus on a specific case which is to be called the head-tail distance. Letting $\ell_1, \ell_2$ be two positive integers for which $\ell_1 + \ell_2 \le L - 1$, the $(\ell_1-\ell_2)$ **head-tail distance**, based on some underlying distance $d$ (such as Euclidean), between the given two windows is defined by

$$d_{HT}(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{\ell_1} d(x_i, y_i) + \sum_{i=L-\ell_2+1}^{L} d(x_i, y_i).$$

Note that one may further prefer to choose the distances in these two summations differently.

**3. Cosine distance** is defined to be

$$d_{\cos}(\mathbf{x}, \mathbf{y}) = 1 - \frac{\sum_{i=1}^{L} x_i y_i}{\|\mathbf{x}\|_2 \|\mathbf{y}\|_2}.$$

Writing $\theta$ for the angle between $\mathbf{x}$ and $\mathbf{y}$, $d_{\cos}(\mathbf{x}, \mathbf{y}) = 1 - \cos \theta$. The cosine distance is known to work effectively in higher dimensions, and it is used in areas such as natural language processing and text mining.

A closely related distance is **Pearson correlation** defined by

$$d_{\text{Pearson}}(\mathbf{x}, \mathbf{y}) = 1 - \frac{\sum_{i=1}^{L} (x_i - \bar{x})(y_i - \bar{y})}{\left(\sum_{i=1}^{L} (x_i - \bar{x})^2\right)^{1/2} \left(\sum_{i=1}^{L} (y_i - \bar{y})^2\right)^{1/2}}$$
Note that correlation distance is not a metric. When $\bar{x} = \bar{y} = 0$, $\|x\|_2 = \|y\|_2 = 1$, $d_{\text{Pearson}}(\mathbf{x}, \mathbf{y})$ simplifies to $1 - \sum_{i=1}^{L} x_i y_i$. Then under the same assumptions, we have $d_{\text{Pearson}}(\mathbf{x}, \mathbf{y}) = d_{\cos}(\mathbf{x}, \mathbf{y})$, and also

$$d_2^2(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{L} (x_i - y_i)^2 = \sum_{i=1}^{L} x_i^2 + \sum_{i=1}^{L} y_i^2 - 2 \sum_{i=1}^{L} x_i y_i = 2d_{\text{Pearson}}(\mathbf{x}, \mathbf{y}).$$

Let us note that another distance similar to Pearson correlation can be defined by using Spearman correlation. The main difference from Pearson correlation is that the variables of interest in this case are rank-ordered.

**4. Canberra distance** is defined to be

$$d_{\text{Canberra}}(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{L} \frac{|x_i - y_i|}{|x_i| + |y_i|}.$$

This distance is a weighted version of the $l^1$ distance, where the weights help robustness to outliers. It is used in various areas including clustering and classification, and has several variations which we will not discuss here.

**5. Others.** There is vast list of other similarity measures in the literature that can be adapted to measure similarity of trajectories. These include but are not limited to mean character difference, Frechet distance, coefficient of divergence, matching subsequence counts, dynamic time warping, longest common subsequence ([2], [18], [32], [31], [34]). In our study we had some preliminary experiments with the longest common subsequences approach. The length of the longest common subsequences is a classical measure of similarity for discrete sequences. It is used, for example, in molecular biology for DNA comparison, computer science for binary sequence comparison, among others. See [4], [30] and [31] for various applications of the longest common subsequences as well as algorithmic discussions. Recent theoretical advances on the common subsequences problem can be traced in [13]. The theoretical background can also be carried over to a score function setting which can be used to treat continuous data. See [11] for relevant discussions and references.

Following the discussion in [34], the variant of the longest common subsequences for comparison of continuous trajectories can be defined as a dynamic programming algorithm as follows. Letting $\mathbf{x}$ and $\mathbf{y}$ be as before, the length of the longest common subsequence $LC(\mathbf{x}, \mathbf{y})$ is defined by

$$LC(\mathbf{x}, \mathbf{y}) = \begin{cases} 0, & \text{if } \mathbf{x} \text{ or } \mathbf{y} \text{ is empty,} \\ 1 + LC(\mathbf{x}[1:k-1], \mathbf{y}[1:m-1]), & \text{if } |x_k - y_m| < \epsilon \text{ and } |k-m| \le \delta, \\ \max\{LC(\mathbf{x}[1:k-1], \mathbf{y}), LC(\mathbf{x}, \mathbf{y}[1:m-1])\}, & \text{otherwise.} \end{cases}$$

Here, the values of the $\epsilon$ and $\delta$ are to be chosen by the researcher. The choice of $\delta = L - 1$ yields the natural analogue of the longest common subsequences in a continuous setup. Our basic experiments with the longest common subsequences yielded poor performance, and we therefore do not list them below. A variant of this subsequence approach is the dynamic time warping for which relevant discussions can be found in [34] and [2]. We hope to get back to the longest common subsequence, dynamic time warping and their variations in more detail in a subsequent work.

## 2.3 Evaluation metrics

In this subsection we briefly describe the evaluation metrics for point and interval forecasts that are used in this manuscript.

**Point forecasts.** Below we will be using the mean absolute error (MAE) and the mean absolute percentage error (MAPE) for the evaluation of our point forecasts. One may refer to [15] and [16]
for general discussions on the point forecast evaluation metrics. MAE corresponding to the forecasts $\hat{y}_1, \dots, \hat{y}_m$ for the realized sequence $y_1, \dots, y_m$ is defined by

$$MAE = \frac{1}{m} \sum_{i=1}^{m} |\hat{y}_i - y_i|.$$

The mean absolute percentage error (MAPE) is defined by

$$MAPE = \frac{1}{m} \sum_{i=1}^{m} \left| \frac{\hat{y}_i - y_i}{y_i} \right|.$$

Some researchers in the traffic literature use other related metrics such as $MAE_{100}$ for the traffic flow prediction, which basically ignores the time instances $t$ for which $y_t < 100$. For example, see [23]. We do not consider a similar strategy here since minimal flow in our case is not too small, and making use of such a metric does not make much difference than using the standard MAE.

Once point forecasts are obtained via two distinct models, one may also like to know whether one method is significantly different than (better or worse) the other one. For this purpose we will be using the Diebold-Mariano test, see [6] and [7] for details.

**Interval forecasts.** The two criteria used for comparison of performance of prediction interval mechanisms below are (i) Unconditional coverage, and (ii) Winkler score. Let us now briefly go over these.

* **Unconditional coverage.** The simplest metric one may consider for a prediction interval is the unconditional coverage, which is just the proportion of samples in test data covered by our prediction intervals. Formally, given the prediction intervals $[\hat{L}_t, \hat{U}_t]$, $t = 1, \dots, S$, with the corresponding actual values $y_1, \dots, y_S$, the unconditional coverage is defined to be

$$UC = \frac{1}{S} \sum_{t=1}^{S} 1(y_t \in [\hat{L}_t, \hat{U}_t]),$$

where $1$ is the characteristic function giving $1$ if the input is true, and $0$ otherwise. For $\alpha \in (0, 1)$, a $(1 - \alpha)100\%$ prediction interval construction will be considered to work well if the UC is close to $(1 - \alpha)100\%$.

Note that there is also the conditional coverage of prediction intervals which also penalizes the clustering of 0's and 1's (for the sequence $1(y_t \in [\hat{L}_t, \hat{U}_t])$). Such clusterings may occur in presence of data with significant outliers. The electricity price forecasting with its spiky behavior is one such example, see [29]. We will not be reporting on conditional coverage below.

* **Winkler score.** Besides the obvious coverage probabilities, one other natural criteria to consider is the length of prediction interval. One naturally wishes to have the prediction intervals as short as possible (keeping the UC around the desired level). For this purpose we will be using Winkler score and penalize the longer prediction intervals [41]. Before giving the definition for the Winkler score of a prediction interval, let us note that the constructions of our prediction intervals are symmetric around the point forecasts. Now, for $\alpha \in (0, 1)$, for the time series $y_t$, $t \ge 1$, and for the corresponding prediction intervals $[\hat{L}_t, \hat{U}_t]$, we follow [29] and define the corresponding Winkler score by setting

$$W_t = \begin{cases} \hat{U}_t - \hat{L}_t & \text{if } y_t \in [\hat{L}_t, \hat{U}_t], \\ \hat{U}_t - \hat{L}_t + \frac{2}{\alpha}(\hat{L}_t - y_t) & \text{if } y_t < \hat{L}_t \\ \hat{U}_t - \hat{L}_t + \frac{2}{\alpha}(y_t - \hat{U}_t) & \text{if } y_t > \hat{U}_t. \end{cases}$$
As is clear from the definition, Winkler score penalizes the cases where the actual value of the time series is outside the prediction interval, and in this case, it is proportional to the distance between the actual value and the closest point on the prediction interval. We may then define the average Winkler score over the time interval 1, ..., S to be

$$ \bar{W} = \frac{1}{S} \sum_{t=1}^{S} W_t. $$

It is also clear from these definitions that a lower Winkler score is an indicator of a better prediction interval performance.

## 2.4 Sample Quantile

For our prediction interval models, we need to estimate upper and lower quantiles from a sample. There are multiple ways to define the corresponding sample quantiles, and in this section we fix one that is to be used throughout the paper. When the sample size is small (e.g. less than 50), a precise calculation of the quantile would be more desirable. Due to this, the (precise) sample quantile is defined as follows. Assume we have an ordered sample {$a_1, a_2, ..., a_N$} of size N, such that $a_i \le a_{i+1}$ for each i. Then i-th element of the sample $a_i$, corresponds to the $\frac{i}{N+1}$-th quantile. With that reasoning we define the q-th sample quantile in the following way. For $q \in [0, 1]$,

$$ a(q) = \begin{cases} a_{q(N+1)}, & \text{if } q(N+1) \in \mathbb{Z}^{\ge 0} \\ c_1 a_{\lfloor q(N+1) \rfloor} + c_2 a_{\lfloor q(N+1) \rfloor}, & \text{otherwise} \end{cases} $$

where $c_1 = 1 - (q(N+1) - \lfloor q(N+1) \rfloor)$ and $c_2 = q(N+1) - \lfloor q(N+1) \rfloor$.

## 3 Analysis for weighted Euclidean case

This section contains a detailed analysis of the general methodology in this paper in a specific setting. Here, following [17], we specialize on the case where the similar trajectories are obtained via a weighted variant of $\ell^2$ distance as explained in Section 2.2 (In particular, the weights will be $w_i = i/(L(L+1)/2)$, $i = 1, ..., L$). In later sections, there will be various approaches in use, but they share certain common characteristics. Focusing on a special case will retain us repeating several issues, and we hope that the study in this section will clearly expose the ideas which are involved in general.

Let us start by recalling the general reasoning of the similarity of trajectories approach. First, we generate trajectories based on lagged observations of traffic flows. Then we identify the nearest neighbors for each trajectory and save the corresponding flow values as candidates. These candidates should closely resemble the flow that we aim to predict, but as usual they may contain outliers, which are to be handled before obtaining point and interval forecasts via the candidate sets. After the outlier removal is done, the remaining candidates are used to arrive at the final forecasts. The following diagram provides a broad description of the methods we use throughout this paper.
<figure><img src="image_5.png" /><figcaption>Figure 3: Flowchart showing the general similarity of trajectories methodology</figcaption></figure>

## 3.1 Obtaining the Nearest Neighbors

The selection of the nearest neighbors is one of the most critical parts of historical similarity methods. In order to have promising results, it is crucial that the nearest neighbors have a common and meaningful pattern. Furthermore, the number of nearest neighbors used must be sufficient enough to provide a good view. Keeping these in mind, we must select optimal values for the three primary hyperparameters:

* Window size, *L*, for defining the trajectories,

* Distance function, *d*, for measuring similarity,

* Number of nearest neighbors, *K*, for obtaining as many neighbors as required in order to have satisfactory relationships.

Recall that we fixed the distance measure to a weighted Euclidean distance throughout this section. The time windows that are to be considered in this section will be ranging from 30 minutes to 5 hours (i.e., 2 to 20 past observations) resulting from 15-minute intervals in our data set. We save up to 200 nearest neighbors for forecasting purposes. Further hyperparameter selection discussions will be given in the sequel.

Now, before proceeding further, let us include an independent note on the distribution of nearest neighbors. One may expect to find seasonal or other distinguishable patterns in historical search, but our experimentations do not reveal any visible pattern at all. In particular, in Figure 4, we present an exemplary distribution of the 10 nearest neighbors for 100 specific consecutive trajectories (To find the nearest neighbors in this specific case, we set the trajectory length to 16 and use the weighted Euclidean distance). Although the most of the nearest neighbors were from more recent observations, we observed that there were several candidates from as far back as a year ago. This may imply that extending the interval of the data may be beneficial for our models.
<figure><img src="image_6.png" /><figcaption>Figure 4: A sample distribution of the nearest neighbors</figcaption></figure>

## 3.2 Benchmarks

In this section we divert from the main discussion in order to describe our benchmark models. We have chosen four methods for this purpose, namely, Naive (or, Random Walk), Autoregressive model (AR), Autoregressive Integrated Moving Average (ARIMA) and Random Forest.

The naive approach simply uses the current flow value as the forecast for the next time step ($X_t = X_{t-1}$). This method sets the baseline performance results and any proposed method is expected to outperform it. A slightly more complicated version of the standard naive method is the seasonal naive, in which the forecast we give is based on an averaging with respect to certain previous time instances (such as previous days or weeks), taking seasonality into account. Since the results of seasonal naive were not close to the ones that are to be described below, we will only provide the ones for standard naive approach.

AR and ARIMA models are econometric models that are widely used in time series analysis and forecasting. AR models use previous values of the target variable to forecast future values. In particular, if the series of interest is $X_t$, then the AR($p$) model can explicitly be written as
$$X_t = \sum_{i=1}^{p} \alpha_i X_{t-i} + \epsilon_t,$$
where $\epsilon_t$ is white noise and $\alpha_1, \dots, \alpha_p$ are the parameters of the model. On the other hand, ARIMA and SARIMA models combine autoregressive and moving average components to account for trends and seasonal patterns in the data. These models are capable of capturing data patterns and dependencies, such as day-of-week effects and hourly variations, and have the potential to provide relatively accurate predictions even with a limited data. As a result, AR and ARIMA serve as robust and simple benchmarks methods for our research. We refer to [12] and [14] for further details.

Random Forest (RF) is widely recognized for its effectiveness in capturing nonlinear relationships within high-dimensional data. In short-term traffic flow forecasting literature, Random Forest has demonstrated competitive performance as an ensemble learning method [46]. Given the resemblance in methodological approaches between our proposed model and Random Forest, together with the demonstrated efficacy of Random Forest in previous applications, we selected it as the primary benchmark here.

The benchmark results we obtained are summarized in the following table.
<table><caption>Table 1: Benchmark results for point forecast</caption><thead><tr><th>Model</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td>Naive</td><td>54.74 / 56.73</td><td>5.35 / 5.36</td></tr><tr><td>AR</td><td>50.34 / 52.62</td><td>4.92 / 4.94</td></tr><tr><td>ARIMA</td><td>49.34 / 51.83</td><td>4.82 / 4.83</td></tr><tr><td>SARIMA</td><td>42.09 / 46.08</td><td>4.18 / 4.39</td></tr><tr><td>RF</td><td>43.58 / 46.85</td><td>4.35 / 4.40</td></tr><tr><td>Seasonal RF</td><td>42.91 / <strong>45.99</strong></td><td>4.31 / <strong>4.32</strong></td></tr></tbody></table>

Note that in Table 1 and below, the result on left-hand-side shows the tune performance and the one on right-hand-side the test performance. Now we continue with some comments on the results of benchmark models, and on some notes relevant to their implementation. First let us mention the setup for ARIMA and AR models which are seen to improve the naive model. For the AR model, we identified the optimal value of *p* to be 9, yielding a mean absolute error (MAE) of 52.62. The ARIMA model, with optimal values of *p* = 9, *d* = 0, and *q* = 3, exhibited superior performance with an MAE of 51.83. Finally, the SARIMA model demonstrated the best performance among the three benchmarks, with optimal hyperparameters *p* = 9, *d* = 0, *q* = 3, *P* = 3, *D* = 0, *Q* = 2, and seasonal frequency *S* = 96, achieving an MAE of 46.08. The results highlight the importance of seasonality.

Random Forest models consist of several hyperparameters. We refer to *ranger* R package [42] for naming and detailed description. In our models, the “mtry” hyperparameter that determines the number of features randomly sampled at each split, ranges from 4 to 11. Meanwhile the “min.node.size” hyperparameter, determining the minimal node size to split at, varies from 2 to 20. There are various studies regarding the “num.trees” hyperparameter that determines the number of trees in each forest. These studies imply that higher values simply yield better performance. However, our experiments indicate that the improvement in performance becomes insignificant after the number of trees exceed 150. Therefore, for our experimentation, the maximum value for “num.trees” in the grid search is set to 200. Also, the models trained with minimizing mean absolute error performed better, compared to the ones that minimize mean square error.

We had two main approaches for the random forest models. The first one only used lagged observations as the features. We tried various lengths for window size and stopped at 30 since no further improvement was observed in the accuracy. The first approach yielded an MAE of 47.01, performing better than non-seasonal autoregressive benchmarks but worse than SARIMA, implying a motivation for a random forest model with seasonal parameters (SRF). In the second approach, we also included observations from one day and one week before. Again we conducted a similar grid search for window sizes and other hyper-parameters. SRF performed best among other benchmarks with MAE of 45.85.

Lastly, the train/test split of benchmark models differ from the splits of similarity models. We keep the size of test split the same, 4 months, and use 12 months of observations for training.<fnref n="1" />

## 3.3 Point Forecasts: Similarity of Trajectories

There are various issues to be considered in providing point forecasts via similarity of trajectories. In this section, focusing on the aforementioned weighted Euclidean distance used in this section, we give a general overview of these, by focusing on hyperparameter tuning, outlier handling and multi-step forecasting.

<footnote n="1">¹As an additional note, using larger sets for training resulted in worse test scores.</footnote>
**Tuning Hyperparameters.** Hyperparameter optimization is a crucial step in forecasting problems. Recall in our case that the hyperparameters are the distance function *d*, the number of neighbors selected *K* and the window size *L*. Since the distance function is already fixed, we are left with adjusting *K* and *L*. Similar tuning is required to be done for the forecasting mechanisms described in Section 5, but not all the details will be presented there.

In our treatment below in this section, the optimal values of window size *L* is searched from 2 to 20. Then for each *L*, average of *K* candidates from nearest neighbors is used as the forecast, where *K* ranges from 5 to 200. The following heat map summarizes our experiment results for this range of *K* and *L*.

<figure><img src="image_7.png" /><figcaption>Figure 5: Initial Grid Searches for *K* and *L*. Left figure shows a wider search. Right figure focuses on the better performing hyperparameters.</figcaption></figure>

In Figure 5, the left hand side shows a general view, and then with zooming in on the right hand side we see a nicely shaped region containing the optimal hyperparameters. This region suggests us to make a second grid search with more frequent values for *K* between 10 and 50, and their pairings with $9 \le L \le 17$. Note that we were conducting this initial experimentation on the tuning split, so the results differ from the final results given in Section 5.

<figure><img src="image_8.png" /><figcaption>Figure 6: Best *L* values</figcaption></figure>

In Figure 6, we observe that $L = 14$ provides the best results and that there is a steep increase in error after $L > 15$ for all values of *K*. Then, while searching for the optimal *K* value, a distinguishable tick pattern can be observed in Figure 7. These observations suggest finalizing the hyperparameter search, and set $L = 14$ and $K = 25$ as the optimal values. While analysing the other aspects of the similarity models, we will mainly be using the optimal hyperparameters. However, we will include results with other values of *K* and *L* whenever it contributes to the discussion.
<figure><img src="image_9.png" /><figcaption>Figure 7: Best *K* values</figcaption></figure>

The following table shows the result of the similarity of trajectories method discussed in this section along with the previously listed benchmark results. Recall that this basic similarity approach was based on using the weighted Euclidean distance and the arithmetic mean of the candidates as forecast.

<table><caption>Table 2: Preliminary comparison of similarity approach with the best benchmark results</caption><thead><tr><th>Model</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td>Naive</td><td>54.74 / 56.73</td><td>5.35 / 5.36</td></tr><tr><td>AR (9)</td><td>50.34 / 52.62</td><td>4.92 / 4.94</td></tr><tr><td>ARIMA (9,0,3)</td><td>49.34 / 51.83</td><td>4.82 / 4.83</td></tr><tr><td>SARIMA (9,0,3), (3,0,2)</td><td>42.09 / 46.08</td><td>4.18 / 4.39</td></tr><tr><td>RF</td><td>43.58 / 46.85</td><td>4.35 / 4.40</td></tr><tr><td>Seasonal RF</td><td>42.91 / <b>45.99</b></td><td>4.31 / <b>4.32</b></td></tr><tr><td>Similarity (L = 14, K = 25)</td><td>45.46 / 47.84</td><td>4.44 / 4.50</td></tr></tbody></table>

Table 2 shows that the similarity of trajectories model outperform autoregressive models (without seasonal parameters) but it performs worse than both Random Forest models and SARIMA. This leads to a motivation for using seasonality in similarity models.

**Handling Outliers.** Extreme values among candidates may lead to biased predictions. Just like any other method, removing these outliers can improve the performance of forecasting based on similarity of trajectories. In our specific case based on traffic data, it is observed that the flow values of candidates for certain observations have an interesting distribution. Figure 8 below shows one such example.
<figure><img src="image_10.png" /><figcaption>Figure 8: Exemplary distribution of candidate flow values. Red line shows the observed flow.</figcaption></figure>

In Figure 8, observe that there is a cluster of outliers having very low value compared to other candidates. We have observed numerous examples sharing similar candidate distribution and addressing this characteristic may be helpful while handling the outliers. [17] uses winsorization for smoothening the extreme values, and they observe that this process improves their forecasting results. The details of winsorization are discussed in Section 4.3. Let us note here that this method is computationally appealing and computational efficiency should be kept as a priority as this procedure is an extra step.

The definition of winsorization deals only with the largest and smallest values, but as seen in Figure 8, there can be many outliers that need to be taken into account. With this in mind, we experimented with various other approaches towards the outlier values. In particular, it is observed that using the z-score to capture outliers is slightly more successful than winsorization in our setting. Besides this, a straightforward approach which we call “Symmetric Tail Removal” is described in Section 4.3, again providing small improvements in results.

Below is a comparison of outlier removal methods with respect to the weighted Euclidean distance, where optimal hyperparameters are used. Let us note that since some of the outliers are removed, the optimal value of *K* may be subject to change. More specifically, here *L* is fixed to be 14 and a grid search among various values of *K* is conducted.

<table><caption>Table 3: Comparison of outlier methods</caption><thead><tr><th>Outlier Method</th><th>L</th><th>K</th><th>MAE</th></tr></thead><tbody><tr><td>None</td><td>14</td><td>25</td><td>45.46 / 47.84</td></tr><tr><td>Winsorization</td><td>14</td><td>21</td><td>45.31 / 47.71</td></tr><tr><td>Z-Score</td><td>14</td><td>27</td><td>45.21 / 47.54</td></tr><tr><td><strong>Sym-Tail-Remove</strong></td><td>14</td><td>25</td><td>45.10 / <strong>47.39</strong></td></tr></tbody></table>

The results in Table 3 suggest that the use of simple outlier removal methods may improve the performance slightly. Although we tried several distinct values for *K*, we have seen that the optimal value tends to be around 25. We observe that symmetric tail removal is the best method for our case. In the following discussions and the next sections we will only use this method.

**Multi-step forecasting.** For obtaining multi-step ahead forecasts, a similar strategy to the one-step ahead forecast is used, but the candidates are chosen according to the number of steps of interest. Since similarity models are non-parametric and they rely on use of patterns, one may hope
that a successful similarity model perform well even when the step-size is huge.

<table><caption>Table 4: Multi-step ahead forecasts</caption><thead><tr><th>Model</th><th>MAE (Step-2)</th><th>MAE (Step-3)</th><th>MAE (Step-5)</th></tr></thead><tbody><tr><td>Naive</td><td>78.160 / 79.889</td><td>98.890 / 101.991</td><td>143.850 / 146.936</td></tr><tr><td>AR</td><td>64.325 / 67.552</td><td>73.414 / 77.103</td><td>106.138 / 111.294</td></tr><tr><td>ARIMA</td><td>58.702 / 59.905</td><td>61.027 / 63.243</td><td>88.823 / 92.032</td></tr><tr><td>RF</td><td>53.975 / 56.480</td><td>59.056 / 62.658</td><td>70.866 / 76.639</td></tr><tr><td>SRF</td><td>53.571 / 54.836</td><td>58.260 / 60.568</td><td>69.863 / 73.694</td></tr><tr><td>Similarity</td><td>53.743 / 54.852</td><td>59.125 / 60.764</td><td>80.308 / 83.789</td></tr><tr><td>Seasonal Similarity</td><td>49.812 / 51.072</td><td>54.297 / 55.516</td><td>61.596 / 62.197</td></tr></tbody></table>

Table 4 demonstrates that the non-seasonal similarity based model performs well when the step size is 2. However, the performance rapidly falls as the step size increases. This could be an indication that the model under consideration misses important patterns from the data. On the other hand, the table contains the results from a “seasonal similarity” model, details of which to be discussed in Section 5.3.1, showing Seasonal Similarity model outperforms all of the benchmarks as the step size increases.

## 3.4 Prediction Intervals

Instead of giving a point forecast, providing a prediction interval may be desirable in many scenarios. To the best of our knowledge, the similarity of trajectories approach has not been explored in the literature in this aspect yet. The purpose of this section is to describe one particular strategy for such a construction. The idea is simple and can be described as follows: After obtaining the candidates, sort them by their value in ascending order and match each observation with the corresponding sample quantile. Afterwards, the required prediction interval is generated using the corresponding quantile forecasts. Note that we use the definition of the sample quantile given in Section 2.4.

Before proceeding further, let us mention the benchmark methods that are to be used in the prediction interval framework. A possible natural approach here is to pass from point forecasts to prediction intervals by making use of the sample standard deviation of the errors. A second possible option is training the models with quantile loss function. There are of course several other directions for obtaining prediction intervals, but our experiments for benchmarking purposes will be based on seven variations around these.

Among these seven, three of them use point forecast benchmarks to generate prediction intervals and the other four are Quantile Autoregression (QAR) [21], Seasonal Quantile Autoregression (SQAR), Quantile Random Forest (QRF) [26] and Seasonal Quantile Random Forest (SQRF) that use quantile loss function. Note that the prediction interval from ARIMA models are only obtained by using sample variance of the point forecast errors.<fnref n="2" />

The optimal hyperparameters for quantile autoregressive models are given in Table 5. QRF and SQRF models are trained simultaneously with their point forecast version using *ranger* R package [42]. The details regarding hyperparameters can be found on Page 12, and the R package follows [26] for obtaining quantile predictions.

We measure performance of these models by using the unconditional coverage and Winkler score defined in Section 2.3. Table 5 shows that the use of quantile loss function provides significantly

<footnote n="2">²Training ARIMA (or moving average models) with quantile loss function does not have a standard definition and we could not find any discussion in the literature. Also, calculation of moving average part using the regression error from a quantile forecast does not seem to be practical.</footnote>
better results. In particular, Seasonal Random Forest trained with quantile loss function. Note that (SD) represents the models that are obtained with adding standard deviation to point forecast models and the numbers near QAR and SQAR represent the optimal autoregressive parameters.

<table><caption>Table 5: Benchmark results for prediction interval forecasts</caption><thead><tr><th>Model</th><th>Unconditional Coverage</th><th>Winkler Score</th></tr></thead><tbody><tr><td>Naive (SD)</td><td>0.9665 / 0.9665</td><td>426.75 / 420.46</td></tr><tr><td>AR (SD)</td><td>0.9671 / 0.9633</td><td>391.38 / 389.16</td></tr><tr><td>ARIMA (SD)</td><td>0.9421 / 0.9393</td><td>378.12 / 373.05</td></tr><tr><td>QAR (9)</td><td>0.9655 / 0.9591</td><td>368.91 / 369.95</td></tr><tr><td>SQAR (8)(6)</td><td>0.9450 / 0.9386</td><td>337.50 / 366.81</td></tr><tr><td>QRF</td><td>0.9291 / 0.9548</td><td>341.08 / 321.89</td></tr><tr><td>SQRF</td><td>0.9349 / 0.9625</td><td>329.34 / 314.48</td></tr></tbody></table>

**Hyperparameter Tuning.** The optimal value of the window size *L* is expected to be similar to the point forecasting case. However, the value of the number of nearest neighbors *K* may require to be larger in order to capture the underlying distribution better. The initial grid search is done similar to the point forecasting case, with *L* between 2 and 20, and *K* between 10 and 200.

<figure><img src="image_11.png" /><figcaption>Figure 9: Unconditional coverage and Winkler score</figcaption></figure>

Figure 9 demonstrates that the unconditional coverage stabilizes for $K \ge 70$. Also the Winkler score does not appear to change significantly for $K \ge 50$, for all values of *L*. The best Winkler score is obtained for *L* = 9 and *K* = 60.
<table><caption>Table 6: Prediction interval results, similarity approach included.</caption><thead><tr><th>Model</th><th>Unconditional Coverage</th><th>Winkler Score</th></tr></thead><tbody><tr><td>Naive (SD)</td><td>0.9665 / 0.9665</td><td>426.75 / 420.46</td></tr><tr><td>AR (SD)</td><td>0.9671 / 0.9633</td><td>391.38 / 389.16</td></tr><tr><td>ARIMA (SD)</td><td>0.9421 / 0.9393</td><td>378.12 / 373.05</td></tr><tr><td>QAR</td><td>0.9655 / 0.9591</td><td>368.91 / 369.95</td></tr><tr><td>SQAR</td><td>0.9450 / 0.9386</td><td>337.50 / 366.81</td></tr><tr><td>QRF</td><td>0.9291 / 0.9548</td><td>341.08 / 321.89</td></tr><tr><td>SQRF</td><td>0.9349 / 0.9625</td><td>329.34 / 314.48</td></tr><tr><td>Similarity</td><td>0.9501 / 0.9451</td><td>333.85 / 329.50</td></tr></tbody></table>

Table 6 shows the comparison between the best similarity model and benchmark models. We see that quantile random forest models are better than the basic similarity approach of this section. A more detailed treatment will be given in Section 6.

# 4 Obtaining and Processing Candidates

## 4.1 The design of the query set and reference sets

Let $L \in \mathbb{Z}^+$ be the window size and $T \in \mathbb{Z}^+$ be the size of the traffic flow data $\{x_1, x_2, \dots, x_T\}$. In this section, by a target, we mean a candidate in our previous setup. Now we design the time series data as trajectories and their targets as follows:

$$ 
\begin{align*}
x_1^{\text{traj}} &= (x_1, x_2, \dots, x_L), & y_1^{\text{targ}} &= x_{L+1} \\
x_2^{\text{traj}} &= (x_2, x_3, \dots, x_{L+1}), & y_2^{\text{targ}} &= x_{L+2} \\
&\vdots & & \\
x_{T-L}^{\text{traj}} &= (x_{T-L}, x_{T-L+1}, \dots, x_{T-1}), & y_{T-L}^{\text{targ}} &= x_T
\end{align*}
$$

We define the query test set as $Q_{\text{test}} = \{[x_i^{\text{traj}}, y_i^{\text{targ}}] : s \le i \le T - L\}$ and the query tune set as $Q_{\text{tune}} = \{[x_i^{\text{traj}}, y_i^{\text{targ}}] : u \le i \le s - 1\}$, where $s$ and $u$ are chosen in our case to be natural numbers satisfying $T - L - s = s - 1 - u$ so that we have $|Q_{\text{tune}}| = |Q_{\text{test}}|$. Note that the data in $Q_{\text{tune}}$ and $Q_{\text{test}}$ correspond to disjoint sets of time intervals. We train models through $Q_{\text{tune}}$ by tuning hyperparameters/parameters and measure accuracies on the test set $Q_{\text{test}}$. Predictions are based on the previously observed trajectories. We pick one trajectory (i.e. $x_i^{\text{traj}}$) from $Q_{\text{tune}}$ and define a reference set for it. We also pick one trajectory from $Q_{\text{test}}$ and again define a reference set for it. For this purpose, we let $R_q^{\text{tune}} = \{[x_i^{\text{traj}}, y_i^{\text{targ}}] : 1 \le i \le q - 1\}$ a reference set for each $q \in \{u, u + 1, \dots, s - 1\}$ and $R_{q'}^{\text{test}} = \{[x_i^{\text{traj}}, y_i^{\text{targ}}] : w \le i \le q' - 1\}$ a reference set for each $q' \in \{s, s + 1, \dots, T - L\}$, where we choose $w$ satisfying $s - 1 = T - L - w$. Choosing $w$ a particular value subject to this condition guarantees the fact that we consider trajectories with dates going back to at most one year for both the trajectories in the query tune set and the query test set. In fact, putting $w = 1$ leads to no problem in practice. However, we preferred to see the results by considering trajectories through a history of equal lengths for both the query tune set and the query test set. Let $q$ be in $\{u, u + 1, \dots, s - 1\}$ and $q'$ be in $\{w, w + 1, \dots, T - L\}$ and let $Z_i = [x_i^{\text{traj}}, y_i^{\text{targ}}]$. The following figure demonstrates the construction of these sets.
<figure><img src="image_12.png" /></figure>

## 4.2 Obtaining the Nearest Trajectories

Below we let $[x_q^{\text{traj}}, y_q^{\text{targ}}] \in Q_{\text{tune}}$. In general, we are willing to forecast $y_q^{\text{targ}}$ and the prediction provided by the similarity model is denoted by $y_q^{\text{pred}}$. In this brief subsection we define the necessary notation in order to explain the selection of the nearest trajectories rigorously.

Let $K \in \mathbb{Z}^+$ be the number of closest trajectories drawn from $R_q^{\text{tune,traj}} = \{x_i^{\text{traj}} : 1 \le i \le q-1\}$ with respect to some distance function $d$. We set $D_q = \{d(x_q^{\text{traj}}, x) : x \in R_q^{\text{tune,traj}}\}$ and denote by $\{d_{q,s}\}_{s=1}^K$ the $K$ smallest numbers of $D_q$ and $\{x_s^{\text{q,near-traj}}\}_{s=1}^K \subset R_q^{\text{tune,traj}}$ the $K$ trajectories which are closest to the trajectory $x_q^{\text{traj}}$. So, $d_{q,s} = d(x_q^{\text{traj}}, x_s^{\text{q,near-traj}})$ for $s = 1, 2, \dots, K$. Here, without loss of generality we consider $d_{q,s} \le d_{q,s+1}$ for $1 \le s \le K-1$. Finally, that gives us a sub-collection of the reference set $R_q^{\text{traj}}$ defined by $C_q = \{[x_s^{\text{q,near-traj}}, y_s^{\text{q,near-targ}}] : 1 \le s \le K\}$, where $y_s^{\text{q,near-targ}}$'s are targetsof $x_s^{\text{q,near-traj}}$'s, equipped with the distances $\{d_{q,s}\}_{s=1}^K$ corresponding to $x_q^{\text{traj}}$.

## 4.3 Handling Outliers

Once the nearest neighbors are selected, outlier detection and removal may as usual help in obtaining better forecasts. Smoothing the possible outliers via winsorization is used in [17] in setting of similar trajectories. Below we briefly discuss the tail removal along with winsorization, although we experimented with other approaches as well.

**Winsorization:** We recall the definition of winsorization from [17]. Assuming the candidate values for our forecast are given by $\{f_1, f_2, \dots, f_k\}$ where $f_i \le f_{i+1}$ for each $i$, the corresponding winsorized candidates are defined by

$$f_i^{(w)} = \begin{cases} f_{i+1}, & \text{if } f_i = \min(f_j) \\ f_{i-1}, & \text{if } f_i = \max(f) \\ f_i, & f_i \notin \{\min(f), \max(f)\} \end{cases}$$

Simply, this procedure replaces the smallest and highest values of the candidates with the values to them. Another natural approach to outliers is the following.

**Tail Removal:** The method removes the smallest $r_1$ and the highest $r_2$ candidates. The values of $r_1$ and $r_2$ can be set in the following ways.

* (Constant) $r_1 = c_1$ and $r_2 = c_2$ for some properly chosen $c_1, c_2$.

* (Percentile) $r_1 = \lfloor \gamma_1 K \rfloor$, $r_2 = \lfloor \gamma_2 K \rfloor$ for some positive $\gamma_1, \gamma_2$ such that $\gamma_1 + \gamma_2 < 1$.

When $\gamma_1 = \gamma_2$, these two can be called as the symmetric constant and symmetric percentile approaches.

# 5 Methodologies for point forecasting

Recall from Section 4.2 how we obtain $C_q = \{[x_s^{\text{q,near-traj}}, y_s^{\text{q,near-targ}}] : 1 \le s \le K\}$ which is the collection of nearest trajectories and the distances $\{d_{q,s}\}_{s=1}^K$ with respect to some distance function
d. We use these for point-forecasting in various ways. A basic method for reaching $y_q^{pred}$ is just a simple averaging:

$$y_q^{pred} = \frac{1}{K} \sum_{s=1}^{K} y_s^{q, \text{near-targ}}.$$

One perspective in our approaches below is to provide forecasts using weighted means of targets $y_s^{q, \text{near-targ}}$ where the weights are obtained in specific ways that are to be described in Section 5.1. Another one is making use of local regression which is discussed in Section 5.2. Before proceeding further into these, let us compare the performance of using different distance functions and forecasting via just the arithmetic mean of the candidates.<fnref n="3" />

To compare different distance functions, first a grid search for the hyperparameters is conducted on tune set. Even values between 2 and 20 were possible candidates for hyperparameter $L$, and the optimal $K$ value is searched from the set {10, 15, 20, ..., 200}. For the seasonal filtered models, $\mathcal{R}$ hyperparameter ranges from 0 to 6.

<table><caption>Table 7: Accuracy comparison with respect to different distances. The bottom half of the table uses seasonal filtering.</caption><thead><tr><th>Distance</th><th>best hyperp.</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td>Canberra</td><td>L = 8, K = 25</td><td>46.80 / 48.69</td><td>4.58 / 4.60</td></tr><tr><td>Correlation</td><td>L = 20, K = 10</td><td>146.02 / 130.58</td><td>13.89 / 11.82</td></tr><tr><td>Euclidean</td><td>L = 8, K = 30</td><td>46.40 / 48.40</td><td>4.54 / 4.57</td></tr><tr><td>Head-Tail</td><td>L = 8, K = 75</td><td>48.16 / 49.45</td><td>4.68 / 4.65</td></tr><tr><td>Maximum</td><td>L = 8, K = 30</td><td>47.01 / 49.08</td><td>4.59 / 4.62</td></tr><tr><td>W. Euclidean</td><td>L = 14, K = 25</td><td>45.46 / <b>47.84</b></td><td>4.44 / <b>4.50</b></td></tr><tr><td>Canberra (s)</td><td>L = 6, K = 25, R = 3</td><td>45.66 / 46.45</td><td>4.64 / 4.40</td></tr><tr><td>Correlation (s)</td><td>L = 20, K = 20, R = 2</td><td>90.98 / 82.46</td><td>9.22 / 8.13</td></tr><tr><td>Euclidean (s)</td><td>L = 8, K = 20, R = 3</td><td>45.41 / 46.66</td><td>4.60 / 4.40</td></tr><tr><td>Head-Tail (s)</td><td>L = 8, K = 20, R = 1</td><td>46.97 / 47.33</td><td>4.80 / 4.49</td></tr><tr><td>Maximum (s)</td><td>L = 8, K = 20, R = 5</td><td>46.08 / 47.21</td><td>4.65 / 4.47</td></tr><tr><td>W. Euclidean (s)</td><td>L = 10, K = 25, R = 3</td><td>44.59 / <b>45.82</b></td><td>4.51 / <b>4.32</b></td></tr></tbody></table>

Table 7 shows that the model using the weighted Euclidean distance is clearly the best among listed. Also, one may observe that applying seasonal filters improves the performance regardless of the distance function used. In the following sections, unless otherwise is stated, the distance function used for obtaining the candidates will be the weighted Euclidean function.

## 5.1 Weights for closest trajectories

Model-1 In general, one would expect that if a trajectory from reference set is closer than another one, then its target should be weighted more. Instead of taking the arithmetic mean of corresponding targets we may consider specific weights and take the weighted average:

$$y_q^{pred} = \sum_{s=1}^{K} w_s y_s^{q, \text{near-targ}},$$

<footnote n="3">³In Table 7 (s) means seasonal filtering and $\mathcal{R}$ represents radius, both of which are to be described in Section 5.3.1.</footnote>
where $w_s > 0$ and $\sum_{s=1}^{K} w_s = 1$. We define the weights via a non-decreasing function $f$ :
$\{1, \dots, K\} \to \mathbb{R}$ by setting $w_s = \frac{f(K-s+1)}{S}$, where $S = \sum_{x=1}^{K} f(x)$. For instance, if $f = 1$, then the prediction is just the arithmetic mean of the targets. Our experiments involved the following choices of the functions: $f_1(x) = 1$, $f_2(x) = x$, $f_3(x) = \sqrt{x}$, $f_4(x) = \ln(1+x)$ and $f_5(x) = \ln^2(1+x)$.

Model-2 Instead of choosing weights uniformly, we may also customize them according to the distances:

$$y_q^{\text{pred}} = \sum_{s=1}^{K} w_{q,s} y_s^{q,\text{near-targ}},$$

where $w_{q,s} = \frac{g(d_{q,s})}{\sum_{s=1}^{K} g(d_{q,s})}$ with $g$ an arbitrary positive decreasing function on $[0, \infty)$. This is based on the intuition that the closer a trajectory is the greater weight of its target should be. The functions experimented consisted of the following: $g_1(x) = \frac{1}{x+0.01}$, $g_2(x) = \frac{1}{\sqrt{x+0.01}}$, $g_3(x) = \frac{1}{x\sqrt{x+0.01}}$ and $g_4(x) = \frac{1}{x^2+0.01}$.

Model-3 We may also determine the weights of $y_s^{q,\text{near-targ}}$'s using linear regression. The linear regression is used in the following form:

$$y = \sum_{s=1}^{K} w_s y_s + w_{K+1}.$$

We find the minimizers as weights $\tilde{w}_1, \dots, \tilde{w}_K, \tilde{w}_{K+1}$ and then consider:

$$y_q^{\text{pred}} = \left[ \sum_{s=1}^{K} \tilde{w}_s y_s^{q,\text{near-targ}} \right] + \tilde{w}_{K+1}$$

The training is processed using the tune query set. Then, we measure the trained model via test query set.

## 5.2 A local regression approach

We may apply linear regression to a local data in contrast to the one discussed in the previous section. Given $x_q^{\text{traj}}$, we consider $C_q = \{[x_s^{q,\text{near-traj}}, y_s^{q,\text{near-targ}}] : 1 \le s \le K\}$ and set up a linear regression model on $C_q$. If we denote by $x_s^{q,\text{near-traj}}(t)$ the $t^{th}$ component of the vector $x_s^{q,\text{near-traj}}$ for $s \in \{1, 2, \dots, K\}$, and $x_q^{\text{traj}}(t)$ the $t^{th}$ component of the vector $x_q^{\text{traj}}$, then we indeed minimize the following function

$$L(w_1, w_2, \dots, w_L, w_{L+1}) = \sum_{s=1}^{K} \left| \left( \sum_{t=1}^{L} w_t x_s^{q,\text{near-traj}}(t) \right) + w_{L+1} - y_s^{q,\text{near-targ}} \right|^2,$$

and use the prediction:

$$y_q^{\text{pred}} = \left( \sum_{t=1}^{L} \tilde{w}_t x_q^{\text{traj}}(t) \right) + \tilde{w}_{L+1},$$

where $\tilde{w}_1, \tilde{w}_2, \dots, \tilde{w}_L, \tilde{w}_{L+1}$ are the optimum weights.
## 5.3 Variations

### 5.3.1 Filtering & Seasonality

An immediate variation of the basic method can be given by first using a filter on the overall observed trajectories in the corresponding reference set, and then doing a search for the similar trajectories over the set of remaining candidate trajectories. One particular strategy that can be used is as follows. Start by fixing a trajectory size, say $L_1$, and a distance function $d_1$, and choose similar trajectories with respect to $d_1$. Calling the collection of chosen trajectories $R$, we may then select the nearest trajectories from $R$ with respect to another distance function $d_2$ (and with possibly a different window size $L_2$).

An important related example would be the case where $R$ holds a certain seasonal characteristics. In general, the set $R$ may consist of the trajectories corresponding to previous days, or weeks. This is implicitly done in [17] where the candidates are only selected from the same time of the previous days.

Extending their work, we define a seasonal filter with a radius hyperparameter, which we denote by $\mathcal{R}$. When $\mathcal{R} = 0$, the set $R$ below will only contain the observations from the same hour and minute of the day. However, as expected, this filter results with a set $R$ with very few trajectories. This is especially undesirable for making interval forecasts. To increase the number of possible candidates, $\mathcal{R}$ is introduced. This hyperparameter loosens the seasonality condition via appending $R$ by adding the consecutive observations of the observations in $R$. As an example, when $\mathcal{R}$ is 0, for an observation made at 15:00, set $R$ would only contain the observations made at 15:00 on previous days. However, setting $\mathcal{R}$ to 1 would add the observations made at 14:45 and 15:15 to set $R$. Our experiments show that such a modification improves the results significantly, especially in interval forecasting.

### 5.3.2 Hourly Model Forecasting

The underlying patterns in different parts of the day may vary from each other and the optimal hyperparameters may change with respect to these underlying patterns. Thus, using multiple models to make forecasts in different parts of the data could improve overall performance. There are numerous ways to split the data. Time of the day is a simple and reasonable way to split in case of forecasting traffic flow. Our experiments reveal that such partitioning and using different models for each hour improves the performance significantly. In particular, in our experiments below, for each hour of the day, hyperparameters of the best performing similarity models are selected on tune set.

Let us provide a brief discussion on the details of experimentation. We used 10-fold cross validation to choose the optimal hyperparameters. For simplicity we fixed the distance to be the weighted Euclidean, used seasonal filtering, and searched for optimal $L, K$ and $\mathcal{R}$ values from 2 to 20, 10 to 200, and 0 to 6, respectively. The grid can be extended to contain other distances. To give a motivation for using multi-hourly models, we provide an hour by hour comparison of some similarity models with different hyperparameters.
<figure><img src="image_13.png" /><figcaption>Figure 10: Hourly performance comparison of different similarity models (Point Forecast).</figcaption></figure>

Figure 10 displays an example of performance differences between two different similarity models in various hours of the day. Overall, the model with optimal hyperparameters (shown by red line) performs better. However, in some hours of the day the other model exhibits significantly better performance.

## 5.4 Results

Table 9 below shows the results of implementations of Model-1, Model-2, Model-3, the Local regression, and their seasonal versions, and lastly the Multi-hourly model, by running a grid search through hyperparameters *L* and *K* shown in Table 8 below. The Euclidean (E) and the weighted Euclidean (WE) distances were used for Model-1, Model-2, Model-3, the Local Regression and their seasonal versions Model-1-S, Model-2-S, Model-3-S and Local R.-S. The weighted Euclidean is the only one used for the multi-hourly model.<fnref n="4" />

<footnote n="4">⁴In the following tables, *R* = 0 except the multi-hourly model.</footnote>
<table><caption>Table 8: Hyperparameters used for grid search</caption><thead><tr><th>Model</th><th>Hyperparameter range</th></tr></thead><tbody><tr><td>Model-1</td><td>L = 2, 3, ..., 20; K = 10, 20, ..., 100; f<sub>1</sub>, f<sub>2</sub>, f<sub>3</sub>, f<sub>4</sub>, f<sub>5</sub>; E, WE</td></tr><tr><td>Model-2</td><td>L = 2, 3, ..., 20; K = 10, 20, ..., 100; g<sub>1</sub>, g<sub>2</sub>, g<sub>3</sub>, g<sub>4</sub>; E, WE</td></tr><tr><td>Model-3</td><td>L = 2, 3, ..., 19, 20; K = 10, 20, ..., 800; E, WE</td></tr><tr><td>Local R.</td><td>L = 2, 3, ..., 20; K = 400, 410, ..., 800; E, WE</td></tr><tr><td>Model-1-S</td><td>L = 2, 3, ..., 20; K = 10, 20, ..., 100; f<sub>1</sub>, f<sub>2</sub>, f<sub>3</sub>, f<sub>4</sub>, f<sub>5</sub>; E, WE</td></tr><tr><td>Model-2-S</td><td>L = 2, 3, ..., 20; K = 10, 20, ..., 100; g<sub>1</sub>, g<sub>2</sub>, g<sub>3</sub>, g<sub>4</sub>; E, WE</td></tr><tr><td>Model-3-S</td><td>L = 2, 3, ..., 20; K = 10, 20, ..., 360; E, WE</td></tr><tr><td>Local R.-S</td><td>L = 2, 3, ..., 20; K = 100, 110, ..., 360; E, WE</td></tr><tr><td>Multi-hourly model</td><td>L = 2, 3, ..., 20, K = 10, 20, ..., 200, R = 0, 1, ..., 6; WE</td></tr></tbody></table>

<table><caption>Table 9: Best results of the models for certain choices of hyperparameters.</caption><thead><tr><th>Model</th><th>best hyperp.</th><th>MAE</th><th>MAPE</th></tr></thead><tbody><tr><td>Naive</td><td>-</td><td>54.74 / 56.73</td><td>5.35 / 5.36</td></tr><tr><td>AR</td><td>p = 9</td><td>50.34 / 52.62</td><td>4.92 / 4.94</td></tr><tr><td>ARIMA</td><td>p = 9, d = 0, q = 3</td><td>49.34 / 51.83</td><td>4.82 / 4.83</td></tr><tr><td>SARIMA</td><td>p = 9, d = 0, q = 3, P = 3, D = 0, Q = 2</td><td>42.09 / 46.08</td><td>4.18 / 4.39</td></tr><tr><td>RF</td><td>See page 12</td><td>43.58 / 46.85</td><td>4.35 / 4.40</td></tr><tr><td>Seasonal RF</td><td>See page 12</td><td>42.91 / 45.99</td><td>4.31 / 4.32</td></tr><tr><td>Model-1</td><td>f<sub>2</sub>, L = 14, K = 40, WE</td><td>45.17 / 47.69</td><td>4.47 / 4.59</td></tr><tr><td>Model-2</td><td>g<sub>4</sub>, L = 13, K = 40, WE</td><td>45.08 / 47.59</td><td>4.39 / 4.47</td></tr><tr><td>Model-3</td><td>L = 15, K = 780, WE</td><td>43.82 / 49.88</td><td>4.33 / 4.70</td></tr><tr><td>Local R.</td><td>L = 19, K = 670, E</td><td>44.19 / 46.74</td><td>4.33 / 4.39</td></tr><tr><td>Model-1-S</td><td>f<sub>2</sub>, L = 4, K = 30, WE</td><td>43.23 / 44.39</td><td>4.43 / 4.23</td></tr><tr><td>Model-2-S</td><td>g<sub>3</sub>, L = 4, K = 40, WE</td><td>43.21 / 44.36</td><td>4.43 / 4.23</td></tr><tr><td>Model-3-S</td><td>L = 4, K = 350, WE</td><td>43.01 / 45.39</td><td>4.40 / 4.32</td></tr><tr><td>Local R.-S</td><td>L = 5, K = 260, E</td><td>41.13 / 43.16</td><td>4.08 / 4.09</td></tr><tr><td>Multi-hourly model</td><td>Multiple Models</td><td>42.88 / 43.28</td><td>4.35 / 4.09</td></tr></tbody></table>

It is seen from Table 9 that the local regression (with similarity) and the multi-hourly model performs better than both the benchmark and other similarity based methods. In particular, let us note that the similarity approach of Model-1-S was used in [17].
<figure><img src="image_14.png" /><figcaption>Figure 11: Diebold-Mariano test results. Dark green means column-model is better than row-model. Black means no significant superiority can be deduced.</figcaption></figure>

Also, one may also like to see whether this conclusion can be shown to hold via a statistical test. When we used the Diablod-Mariano test to compare the benchmarks, we have observed that seasonal random forest model is significantly the best benchmark model. Comparing Model-1-S with benchmarks revealed that Model-1-S is significantly better than the naive, ARIMA and random forest models. However, we have observed no superiority comparing Model-1-S with SARIMA and seasonal random forest. Multi-hourly model significantly outperforms the benchmarks and Model-1-S. Finally, seasonal local regression model exhibited a superior performance against all of the other models shown in Table 9.

# 6 Prediction intervals based on similarity

The purpose of this section is to discuss our approaches on obtaining prediction intervals using similarity of trajectories. The section begins with a general discussion on existing methods in the literature, and then introduces two approaches which we call ST (Similarity Trajectory) and MDST (Model Dependent Similarity Trajectory).

## 6.1 Background

Let $X_t$, $t \ge 1$, be the time series of interest. One naive standard strategy for obtaining a prediction interval is the historical simulation approach. In this case, there is some underlying model (e.g., ARIMA, LSTM,...) that provides point forecasts for the next time instance. Let us write $F_t$ for the model predictions provided by this model. Now suppose that we are at some time $T \in \mathbb{N}$, and that we would like to obtain a prediction interval for time $T+1$ based on observations corresponding to times $1, \dots, T$.<fnref n="5" /> For some predetermined $L \in \mathbb{N}$, the historical simulation approach begins by considering the actual values $X_t$, and the model forecasts $F_t$, for $t \in \{T - L + 1, \dots, T\}$, and looks

<footnote n="5">⁵As in point forecasts, our prediction intervals will be constructed solely based on the earlier observations for the time series of interest, and there will be no exogenous data used. Variations that use similarity of certain other data could improve the results here. In case of the traffic forecasting problem, for example, the data obtained from some other station could be considered as the exogeneous data.</footnote>
at the corresponding errors $\epsilon_t = F_t - X_t$. Then, for given $\alpha \in (0, 1)$, in order to obtain a $(1-\alpha)100\%$ prediction interval for $X_{T+1}$, we choose the sample quantiles $\epsilon(\alpha/2)$, $\epsilon(1 - \alpha/2)$ of the sequence $\epsilon_{T-L+1}, \dots, \epsilon_T$. Then the $(1-\alpha)100\%$ prediction interval is given by

$$ (F_{T+1} + \epsilon(\alpha/2), F_{T+1} + \epsilon(1 - \alpha/2)). $$

We may summarize the discussion on historical simulation prediction intervals, which we call HS, as follows:

HS1. Consider the errors $\epsilon$ made done once a specific model is used for point forecasts.

HS2. Fix some window size $K$, and consider the most recent $K$ errors. Choose the corresponding $\alpha/2$ and $1 - \alpha/2$ sample quantiles.

HS3. Construct the prediction interval by adding these quantiles to the point forecasts provided by the model.

The choice of $L$ in this standard setup can be done by using grid search. Also let us note that a more reasonable way to choose the errors to compute the quantiles could be the use of seasonality. For example, instead of looking at the previous $L$ errors, one could go back further, and look at the $L$ errors corresponding to that instance of previous $L$ days. Such an approach will be called the seasonal historical simulation (HS-S).

Below, we will be using the HS method as a benchmark method in our comparisons. However, there are various other techniques for obtaining prediction intervals that are not discussed below. These include but are not restricted to the distributional prediction interval approach, bootstrap method, and the quantile regression averaging (QRA). Among these, QRA, in which the regressors are the forecasts of certain individual models, is of particular interest to us. In an independent work [1] in preparation, we combine the basic ideas behind QRA and the similar trajectories, and introduce some novel methodologies for providing prediction intervals.

## 6.2 Similarity based prediction intervals

In this subsection we propose two methods for obtaining prediction intervals. The following flowchart provides a general guideline for the tools that are used in these two approaches: ST and MDST.
<figure><img src="image_15.png" /><figcaption>Figure 12: Summary of ST and MDST</figcaption></figure>

The first method we propose is analogous to the case of similarity based point forecasting, and it will be called the similarity of trajectories (ST) approach. It consists of the following steps:

ST1. Fix a distance function, a window size, and choose K most similar trajectories in the past.

ST2. Consider the set of possible candidates {$f_1, \dots, f_K$}, corresponding to these K similar trajectories, as the actual forecast.

ST3. Compute the $f(\alpha/2)$, $f(1 - \alpha/2)$ sample quantiles and form the prediction interval as

$$ (f(\alpha/2), f(1 - \alpha/2)). $$

Note that the described method is different than the historical simulation approach in the sense that it is independent of any model assumptions. No probabilistic assumptions are made, and the given methodology is totally deterministic in nature. When the trajectories are selected among a filtered data according to seasonality (Section 5.3.1), the method will be called ST-S, and when the forecasting is done on an hourly basis as in Section 5.3.2, it will be called multi-hourly model.

Beyond the ST approach, one may wonder how we may combine the HS and ST methods for obtaining prediction intervals. Towards this direction, we introduce and experiment on the following strategy, which we call the model dependent similarity trajectory (MDST):

MDST1. Consider the error sequence $\epsilon$ discussed in HS approach.

MDST2. Fix a distance function, a window size, and choose K most similar trajectories in the past, with respect to the corresponding part of the error sequence.

MDST3. Consider the set of possible candidates {$\epsilon_1, \dots, \epsilon_K$} for the errors.

MDST4. Form the prediction interval by adding the $\alpha/2$ and $1 - \alpha/2$ sample quantiles of the candidate errors to the point forecast provided by the model.

When daily seasonal seasonality is involved, the method will be called MDST-S. This seasonal variation includes the trajectories only starting from the same time of the day as possible candidates, as described in Section 5.3.1. This approach improves the performance significantly as can be seen in the Table 10.
## 6.3 Results

The purpose of this section is to experiment on the performance of using similarity of trajectories in obtaining interval forecasts. Table 10 compares the benchmark and similarity based models using unconditional coverage and the Winkler score. The models are trained to make interval forecasts with 95% coverage. While unconditional coverage shows the accuracy of the models in terms of covering 95% of the observations, the Winkler score is helpful for comparing these models with each other.

<table><caption>Table 10: Prediction interval performance comparison</caption><thead><tr><th>Model</th><th>Best Hyperp.</th><th>Unconditional Coverage</th><th>Winkler Score</th></tr></thead><tbody><tr><td>Naive</td><td>—</td><td>0.9665 / 0.9665</td><td>426.75 / 420.46</td></tr><tr><td>QAR</td><td>p = 9</td><td>0.9655 / 0.9591</td><td>368.91 / 369.95</td></tr><tr><td>SQAR</td><td>p = 8, P = 6</td><td>0.9450 / 0.9386</td><td>337.50 / 366.81</td></tr><tr><td>QRF</td><td>See page 16</td><td>0.9443 / 0.9548</td><td>315.75 / 321.89</td></tr><tr><td>SQRF</td><td>See page 16</td><td>0.9433 / 0.9625</td><td>304.95 / 314.48</td></tr><tr><td>HS</td><td>L = 60</td><td>0.9145 / 0.9168</td><td>388.46 / 374.24</td></tr><tr><td>HS-S</td><td>L = 60</td><td>0.9165 / 0.9208</td><td>335.26 / 327.11</td></tr><tr><td>ST</td><td>L = 9, K = 60</td><td>0.9427 / 0.9358</td><td>333.85 / 329.50</td></tr><tr><td>ST-S</td><td>L = 4, K = 150, R = 5</td><td>0.9499 / 0.9457</td><td>349.43 / 317.32</td></tr><tr><td>Multi-hourly Model</td><td>Multiple Models</td><td>0.9498 / 0.9424</td><td>335.99 / 298.44</td></tr><tr><td>MDST</td><td>L = 8, K = 220</td><td>0.9301 / 0.9308</td><td>340.35 / 326.20</td></tr><tr><td>MDST-S</td><td>L = 8, K = 220, R = 6</td><td>0.9410 / 0.9406</td><td>320.53 / 314.56</td></tr></tbody></table>

The results from Table 10 suggests that the best approach among listed is the multi-hourly model which is the ST-S method in which forecast models are obtained in an hourly basis. None of the other similarity models were better than seasonal random forest but we have seen that the MDST-S and seasonal random forest model exhibit competitive performance with each other.

<figure><img src="image_16.png" /><figcaption>Figure 13: Diebold Mariano Test - Hour based comparison for prediction intervals.</figcaption></figure>
Also, one may also like to see whether these conclusions can be shown to hold via a statistical
test. Following [29], we used Diebold-Mariano test for comparing the performance of interval fore-
casts via their Winkler score. First, we conducted the test on the whole data (as in Section 5.4).
The test resulted in multi-hourly model significantly outperforming all other models and among the
benchmarks, seasonal random forest was the best model. Also, MDST-S turns out to outperform all
models except multi-hourly model and seasonal random forest, and no conclusions could be drawn
between MDST-S and seasonal random forest.

Secondly, using the hour based Diebold-Mariano test approach, again similar to [29], we have
seen that multi-hourly model performs better than all other models in most hours of the day. As
demonstrated in Figure 13, following the multi-hourly model, MDST-S, ST-S and seasonal random
forest were the successful ones. Lastly, let us note that multi-hourly model was superior against
seasonal random forest in only 13 hours of the day, and it was significantly worse in 1 hour of the
day. So the test resulted inconclusive for 10 hours. The competition between the two models in
these hours of the day motivates the search for possible conditions where similarity models work
well and the ones open to improvement.

# 7 Conclusion

The literature on time series forecasting, and in particular on traffic flow forecasting, is very rich, with various approaches and methodologies. Our motivation in this work was to have a general look at the use of similarity of trajectories, and to discuss its possible variations. The experiments consisted of obtaining prediction intervals together with point forecasts.

Our conclusions above can be summarized as follows:

* The methods used in this manuscript seem to yield competitive results against the more sophisticated models such as the random forest. The similarity of trajectories approach has a great flexibility and provides the chance to explain the forecast by a direct reference to historic behavior.

* Using local regression and including seasonality improve the point forecast results significantly. However, each additional step comes up with new hyperparameters and computational issues, and these can be analyzed in an upcoming work.

* Multi-step forecasts using seasonal similarity seem to perform very well, and this can be further analyzed in a separate experimental study.

* Obtaining prediction intervals based on similar trajectories in the same manner seem to be a natural choice, and it yielded results comparable to the ones we used as benchmarks, but we believe that the performance of the proposed strategies are open to improvements.

* Furthermore, in point forecasting, we may customize local regression by replacing linear regressing with other machine learning algorithms (e.g., gradient boosting, support vector machine) for locally obtained data.

In conclusion, the methods we surveyed around are flexible, and are open to further improve-
ments. Our plan in subsequent work is to specialize on certain aspects of the problem and develop
the current strategies in different perspectives via the use of techniques from the statistics and
machine learning literature.

**Acknowledgements** The authors would like to thank Elif Yılmaz for helpful discussions. Second and third authors are supported partially by BAP grant 20B06P.
References

[1] Arslan, İlker, Can Hakan Dağıdır and Ümit Işlak, "On a combination of quantile regression and similar trajectories for predictions intervals", in preparation.

[2] Berndt, Donald J., and James Clifford. "Using dynamic time warping to find patterns in time series." KDD workshop. Vol. 10. No. 16. 1994.

[3] Bolshinsky, Ella, and Roy Friedman. Traffic flow forecast survey. No. CS Technion report CS-2012-06. Computer Science Department, Technion, 2012.

[4] Capocelli, Renato M., ed. Sequences: Combinatorics, Compression, Security, and Transmission. Springer Science & Business Media, 2012.

[5] Chen, Chenyi, et al. "Short-time traffic flow prediction with ARIMA-GARCH model." 2011 IEEE Intelligent Vehicles Symposium (IV). IEEE, 2011.

[6] Diebold, F. X. and Mariano, R. S. (1995). Comparing predictive accuracy. Journal of Business and Economic Statistics, 13: 253-63.

[7] Diebold, Francis X. "Comparing predictive accuracy, twenty years later: A personal perspective on the use and abuse of Diebold–Mariano tests." Journal of Business & Economic Statistics 33.1 (2015): 1-1.

[8] Gilboa, Itzhak, et al. "Economic models as analogies." The Economic Journal 124.578 (2014): F513-F533.

[9] Gilboa, Itzhak, Offer Lieberman, and David Schmeidler. "Empirical similarity." The Review of Economics and Statistics 88.3 (2006): 433-444.

[10] Gilboa, Itzhak, Offer Lieberman, and David Schmeidler. "A similarity-based approach to prediction." Journal of Econometrics 162.1 (2011): 124-131.

[11] Gong, Ruoting, Christian Houdré, and Ümit Işlak. "A central limit theorem for the optimal alignments score in multiple random words." arXiv preprint arXiv:1512.05699 (2015).

[12] Hamilton, James Douglas. Time series analysis. Princeton university press, 2020.

[13] Houdré, Christian, and Ümit Işlak. "A central limit theorem for the length of the longest common subsequences in random words." Electronic Journal of Probability 28 (2023): 1-24.

[14] Hyndman, Rob J., and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2018.

[15] Hyndman, Rob J. "Measuring forecast accuracy." Business forecasting: Practical problems and solutions (2014): 177-183.

[16] Hyndman, Rob J., and Anne B. Koehler. "Another look at measures of forecast accuracy." International journal of forecasting 22.4 (2006): 679-688.

[17] Habtemichael, Filmon G., and Mecit Cetin. "Short-term traffic flow rate forecasting based on identifying similar traffic patterns." Transportation research Part C: emerging technologies 66 (2016): 61-78.

[18] Işlak, Ümit, and Alperen Y. Özdemir. "On an alternative sequence comparison statistic of Steele." Discrete Mathematics & Theoretical Computer Science 22.Combinatorics (2020).
[19] Kumar, Kranti, Manoranjan Parida, and V. K. Katiyar. "Short term traffic flow prediction for a non urban highway using artificial neural network." Procedia-Social and Behavioral Sciences 104 (2013): 755-764.

[20] Kang, Danqing, Yisheng Lv, and Yuan-yuan Chen. "Short-term traffic flow prediction with LSTM recurrent neural network." 2017 IEEE 20th international conference on intelligent transportation systems (ITSC). IEEE, 2017.

[21] Koenker, Roger, and Zhijie Xiao. "Quantile autoregression." Journal of the American statistical association 101.475 (2006): 980-990.

[22] Kumar, S. Vasantha, and Lelitha Vanajakshi. "Short-term traffic flow prediction using seasonal ARIMA model with limited input data." European Transport Research Review 7.3 (2015): 1-9.

[23] Lippi, Marco, Matteo Bertini, and Paolo Frasconi. "Short-term traffic flow forecasting: An experimental comparison of time-series analysis and supervised learning." IEEE Transactions on Intelligent Transportation Systems 14.2 (2013): 871-882.

[24] Liu, Zelei, et al. "A novel trajectory similarity-based approach for location prediction." International Journal of Distributed Sensor Networks 12.11 (2016): 1550147716678426.

[25] Man-Chun, Tan, Feng Luo-Bin, and X. U. Jian-Min. "Traffic flow prediction based on hybrid ARIMA and ANN model." China Journal of Highway and Transport 20.4 (2007): 118.

[26] Meinshausen, Nicolai, and Greg Ridgeway. "Quantile regression forests." Journal of machine learning research 7.6 (2006).

[27] Miglani, Arzoo, and Neeraj Kumar. "Deep learning models for traffic flow prediction in autonomous vehicles: A review, solutions, and challenges." Vehicular Communications 20 (2019): 100184.

[28] Nowotarski, Jakub, and Rafał Weron. "Computing electricity spot price prediction intervals using quantile regression and forecast averaging." Computational Statistics 30.3 (2015): 791-803.

[29] Nowotarski, Jakub, and Rafał Weron. "Recent advances in electricity price forecasting: A review of probabilistic forecasting." Renewable and Sustainable Energy Reviews 81 (2018): 1548-1568.

[30] Pevzner, P. A. "Computational molecular biology: An algorithmic approach A Bradford Book." (2000).

[31] Sankoff, David. "Time warps, string edits, and macromolecules." The Theory and Practice of Sequence Comparison, Reading (1983).

[32] Shirkhorshidi, Ali Seyed, Saeed Aghabozorgi, and Teh Ying Wah. "A comparison study on similarity and dissimilarity measures in clustering continuous data." PloS one 10.12 (2015): e0144059.

[33] Su, Han, et al. "A survey of trajectory distance measures and performance evaluation." The VLDB Journal 29.1 (2020): 3-32.

[34] Tao, Yaguang, et al. "A comparative analysis of trajectory similarity measures." GIScience & Remote Sensing 58.5 (2021): 643-669.
[35] Tian, Yan, et al. "LSTM-based traffic flow prediction with missing data." Neurocomputing 318 (2018): 297-305.

[36] Toohey, Kevin, and Matt Duckham. "Trajectory similarity measures." Sigspatial Special 7.1 (2015): 43-50.

[37] Xia, Ying, and Jungang Chen. "Traffic flow forecasting method based on gradient boosting decision tree." 2017 5th International Conference on Frontiers of Manufacturing Science and Measuring Technology (FMSMT 2017). Atlantis Press, 2017.

[38] Wang, Tianyi. Trajectory similarity based prediction for remaining useful life estimation. University of Cincinnati, 2010.

[39] Wang, Yi, et al. "Combining probabilistic load forecasts." IEEE Transactions on Smart Grid 10.4 (2018): 3664-3674.

[40] Wang, Zhiyuan, Shouwen Ji, and Bowen Yu. "Short-term traffic volume forecasting with asymmetric loss based on enhanced KNN method." Mathematical Problems in Engineering 2019 (2019).

[41] Winkler, Robert L. "A decision-theoretic approach to interval estimation." Journal of the American Statistical Association 67.337 (1972): 187-191.

[42] Wright MN, Ziegler A (2017). "ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R." Journal of Statistical Software, 77(1), 1-17. doi:10.18637/jss.v077.i01.

[43] Yang, Bailin, et al. "Traffic flow prediction using LSTM with feature enhancement." Neurocomputing 332 (2019): 320-327.

[44] Yu, Haiyang, et al. "A special event-based K-nearest neighbor model for short-term traffic state prediction." Ieee Access 7 (2019): 81717-81729.

[45] Zeng, Dehuai, et al. "Short term traffic flow prediction using hybrid ARIMA and ANN models." 2008 Workshop on Power Electronics and Intelligent Transportation System. IEEE, 2008.

[46] Zhang, Lizong, et al. "A hybrid forecasting framework based on support vector regression with a modified genetic algorithm and a random forest for traffic flow prediction." Tsinghua Science and Technology 23.4 (2018): 479-492.

[47] Zhao, Erlong, Pei Du, and Shaolong Sun. "Historical pattern recognition with trajectory similarity for daily tourist arrivals forecasting." Expert Systems with Applications 203 (2022): 117427.