Cyril Furtlehner

Jean-Marc Lasgouttes

Alessandro Attanasi

Marco Pezzulla

Guido Gentile
email: guido.gentile@uniroma1.it

Short-term Forecasting of Urban Traffic using Spatio-Temporal Markov Field

Keywords: Traffic prediction, Gaussian belief propagation algorithm, statistical learning, Markov random field

come

Introduction

The rapidly evolving sector of traffic management information systems has made the problem of traffic forecasting at the urban network level a central and key issue. Traffic management systems aim to monitor, allow real-time scenario evaluations, and find control actions that keep the traffic flow as fluid as possible. They generally operate on main axes or principal sub-networks of the transportation infrastructure, and typically rely on traffic flow models, which in principle are able to predict the onset and the propagation of congestion. To be accurate, these models must have a high level of detail, which requires calibrating an enormous amount of parameters. But this calibration is so challenging [START_REF] Jha | Development and calibration of a large-scale microscopic traffic simulation model[END_REF][START_REF] Osorio | High-dimensional offline origin-destination (OD) demand calibration for stochastic traffic simulators of large-scale road networks[END_REF] and time-consuming that the effectiveness of flow models can be greatly reduced. In this respect, data driven models replace the calibration effort with a simpler tuning of the hyperparameters.

Moreover, the massive increase in available traffic data over the last decade has triggered a surge of interest in traffic forecasting based on data-driven models. General references on these various approaches can be found in [START_REF] Vlahogianni | Short-term traffic forecasting: Where we are and where we're going[END_REF]. These can be divided into two main categories:

• parametric: autoregressive models, probabilistic models, Bayesian models, models based on Markov random fields (MRF);

• non-parametric: k nearest neighbors (k-NN), random forest, Gaussian process, support vector regression, neural networks.

Parametric models, based on ordinary statistical considerations, are more traditional. They are sometimes preferred to their non-parametric counterparts owing to their interpretability. Machine learning (ML) can potentially offer a very large variety of non-parametric models with a wide range of complexity and potential efficiency. Many methods are targeted toward independent modeling of road segments. Methods attempting to leverage spatial dependencies have been recently growing considerably in number [START_REF] Ermagun | Spatiotemporal traffic forecasting: review and proposed directions[END_REF]. If we focus more specifically on forecasting models which attempt to address the problem at the network scale, the requirements we find for such models to be deployed in online applications are the following:

• accuracy: predictions should be significantly better than a simple persistent predictor combined with an historical time-of-day-dependent average, for instance, used when data are incomplete;

• missing data: we cannot expect to have a complete information of the network state at any time, which means that both the training and the running of the model have to be able to take missing data into account;

• scaling: the model should scale up to high systems size, i.e. networks of the size of an urban area, where the number of road segments involved can be around one or two hundred thousand. If we think in terms of detectors, this requirement might be lower. Currently, the number of effective detectors covering a given urban area is smaller by one or two orders of magnitude than the number of road segments (see e.g. [START_REF] Loder | Understanding traffic capacity of urban networks[END_REF]).

Traditional methods based on autoregressive models [START_REF] Levin | On forecasting freeway occupancies and volumes (abridgment)[END_REF][START_REF] Lippi | Short-term traffic flow forecasting: An experimental comparison of time-series analysis and supervised learning[END_REF][START_REF] Fusco | Short-term speed predictions exploiting big data on large urban road networks[END_REF] or hybrid ML methods (see e.g. [START_REF] Sun | Network-scale traffic modeling and forecasting with graphical lasso and neural networks[END_REF]) have recently been adapted to this context, without being able to scale up well to large network sizes. More recently, deep learning approaches have been proposed [START_REF] Lv | Traffic flow prediction with big data: A deep learning approach[END_REF][START_REF] Yu | Spatiotemporal recurrent convolutional networks for traffic prediction in transportation networks[END_REF][START_REF] Li | Diffusion convolutional recurrent neural network: Data-driven traffic forecasting[END_REF] showing convincing performance compared to traditional methods at the price of heavy computational costs [START_REF] Polson | Deep learning for short-term traffic flow prediction[END_REF][START_REF] Manibardo | Deep learning for road traffic forecasting: Does it make a difference?[END_REF] and requiring complete data or to be combined with imputation methods [START_REF] Duan | An efficient realization of deep learning for traffic data imputation[END_REF][START_REF] Zhang | Multistep speed prediction on traffic networks: A deep learning approach considering spatio-temporal dependencies[END_REF].

Beside that, it is shown in [START_REF] Kataoka | Traffic data reconstruction based on Markov random field modeling[END_REF] that a simple Gaussian multivariate model with very few parameters can already perform better than simple interpolation methods at large scale with a cost that is linear w.r.t. system size, thanks to mean-field techniques of statistical physics [START_REF] Yedidia | Generalized belief propagation[END_REF], on which the present paper is also relying. More refined models, like that proposed in [START_REF] Hara | Network-wide traffic state estimation using a mixture Gaussian graphical model and graphical lasso[END_REF] based on a mixture of sparse Gaussian random variables, cannot cope with inference on large scale graphs.

Our Approach

This paper proposes a different direction that copes with missing data, based on Markov random field (MRF) modeling. A graphical model is used to couple observations on segments at different locations and different time steps, such that missing and future observations are treated alike. One difficulty is to train the model offline from historical data; another one is to be able to run it in real time, i.e. to perform the probabilistic inference of all the missing/future variables from the observed ones. In a previous study, we investigated [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF] the possibility of building an MRF which could encode both spatial and temporal dependencies and come with a linear computational time when running the probabilistic inference task. We considered two types of models involving either latent binary variables [START_REF] Martin | Latent binary MRF for online reconstruction of large scale systems[END_REF] or Gaussian copula models. Both come with an associated learning algorithm to generate models adapted respectively to Generalized Belief Propagation for binary variables [START_REF] Furtlehner | Cycle-based cluster variational method for direct and inverse inference[END_REF] and Gaussian Belief Propagation for realvalued variables [START_REF] Martin | GMRF estimation under topological and spectral constraints[END_REF].

Gathering all these elements leads us to propose an approach having the following features, which, to our knowledge, are not simultaneously provided by any other method:

• provide accurate predictions, by exploiting spatial and temporal dependency;

• deal with up to 80% of missing data without drastic degradation of performances;

• deal with online constraints for large network sizes, owing to the linear scaling of belief propagation;

• deal with inhomogeneous data, e.g. speed, flow, travel time, etc., thanks to the definition of a generic traffic index;

• incorporate in an economical way, through the aforementioned traffic index, contextual data like e.g. time of day and seasonal dependencies;

• provide consistent confidence intervals along with the predicted value.

The paper is organized as follows: Section 3 gives the necessary materials, respectively the Gaussian belief propagation algorithm used for inference, the -IPS algorithm used to construct the model and the copula encoding. The workflow of the method is described in Section 4. The results of experiments performed on two different traffic datasets with various features are then reported and analyzed in Section 5. Comparison will be made in particular with results from single detector times series [START_REF] Attanasi | A hybrid method for real-time short-term predictions of traffic flows in urban areas[END_REF] forecast models to measure the advantage gained from spatial dependencies. Section 6 presents a statistical analysis of the errors of the model in which we try to separate intrinsic uncertainty of the traffic phenomena from systematic errors based on modeling issues. Finally, Section 7 outlines future work.

Material and methods

Gaussian Belief propagation

We consider a set V of nodes associated to discrete random variables x = {x i , i ∈ V}, obeying a joint probability distribution of the form

P(x) = ij∈E ψ ij (x i , x j ) i∈V φ i (x i ), (1)
where E ⊂ V × V is a set of edges and φ i and ψ ij are functions associated respectively to a single variable x i and to an unordered pair of variables ij ∈ E.

The ψ ij are called the "factors" while the φ i are there by convenience and could be reabsorbed in the definition of the factors. E together with V define the graph G, which will be assumed to be connected, and ∂i denotes the set of neighbors of node i in G.

In the case where the graph does not have loops (i.e., it is a tree), the set of marginal distributions, called the belief, b i (x) = P(x i = x) associated to each variable i can be computed efficiently. The Belief Propagation algorithm (BP) [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF] does this for all variables in one single procedure, by taking into account that the computation of each of these marginals involves intermediate quantities called the messages m i→j (x j ) "sent" by node i to node j, and which are also necessary to compute other marginals. The idea of BP is to compute all these messages simultaneously, using the relation between them as a fixed point equation. The following message updates

m i→j (x j ) ← xi φ i (x i )ψ ij (x i , x j ) k∈∂i,k =j m k→i (x i ), (2)
are iterated until a fixed point is reached. Then the beliefs are given by

b i (x i ) = 1 Z i φ i (x i ) j∈∂i m j→i (x i ), (3)
where Z i is a normalization constant.

When the graph does have loops, it is still possible to apply the algorithm above. In this case, when convergence occurs, b i (x) is only an approximation of P(x i = x).

In practice, these equations are only usable in two situations: when each x i takes a finite number of values, and when x is a Gaussian vector. These are the only cases where the message update formulas can be parameterized easily. The case of binary values, which correspond to an Ising model, has been studied in [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF][START_REF] Martin | Latent binary MRF for online reconstruction of large scale systems[END_REF]. Our focus here is on Gaussian variables, for which the algorithm takes on a specific form, referred to as Gaussian Belief Propagation (GaBP) [START_REF] Bickson | Gaussian Belief Propagation: Theory and Application[END_REF]. For such variables, the factors are naturally parameterized as

ψ ij (x i , x j ) = exp -A ij x i x j , (4)
φ i (x i ) = exp - 1 2 A ii x 2 i + h i x i , (5)
where the precision matrix A = (A ij ) is the inverse of the correlation matrix of variable x and h = (h i ) is such that A -1 h is the expected value of the variable x. Update rule (2) makes sense for these factors with messages of the form

m i→j (x j ) = exp - (x j -µ i→j ) 2 2σ i→j . (6)
Information is then propagated via the 2-component real vector (µ i→j , σ i→j ) representing bias and variance with the following update rules:

µ i→j ←- 1 A ij h i + k∈∂i,k =j µ k→i σ k→i , (7)
σ i→j ←-- 1 A 2 ij A ii + k∈∂i,k =j σ -1 k→i . (8)
At convergence, the beliefs take the form:

b i (x) = σ i 2π exp - (x -µ i ) 2 2σ i , (9)
with

µ i = σ i h i + j∈∂i µ j→i σ j→i , (10)
σ -1 i = A ii + j∈∂i σ -1 j→i . (11)
There is at most one fixed point, even in the presence of loops. It might not necessarily be stable but, if convergence occurs, µ i is equal to the expectation E(x i ) at each node i [START_REF] Weiss | Correctness of belief propagation in Gaussian graphical models of arbitrary topology[END_REF]. The computation time is O(N log(N )) with the system size N .

-IPS algorithm for model selection

In order to use GaBP, it is first necessary to build a model in the form (4)- [START_REF] Loder | Understanding traffic capacity of urban networks[END_REF]. Since GaBP may often encounter convergence issues, especially with non-sparse structures, it can be of practical interest to construct off-line a Gaussian MRF (GMRF) that is compatible with GaBP. Some features of the matrix A can be desirable to obtain convergence, albeit without much guarantee: sparsity, specific topological properties like the absence of short loops or spectral properties, e.g. walk-summability [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF]. By combining various methods proposed in the context of sparse inverse covariance matrix estimation [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Scheinberg | Learning sparse Gaussian Markov networks using a greedy coordinate ascent approach[END_REF], one way to do that has been developed in [START_REF] Martin | GMRF estimation under topological and spectral constraints[END_REF] in the form of the -IPS algorithm1 . -IPS incorporates the desired features explicitly, by combining an approach based on the iterative proportional scaling (IPS) procedure [START_REF] Speed | Gaussian Markov distributions over finite graphs[END_REF], with block updating techniques used in [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF].

The rationale is, instead of inverting the covariance matrix, to construct the matrix A link by link, by selecting the link that maximizes the log-likelihood of the corresponding distribution, and ensuring at each step that the constraints are satisfied. Several constraints can be considered:

• -LOOP forbids loops of size up to ;

• -FLOOP is the same, but only for frustrated loops, i.e. loops with a negative product of partial covariances (-A ij );

• WS (walk-summability) and WWS (weak walk-summability) are spectral constraints, which involve definite positiveness of some matrices.

Each modification is accepted only if it satisfies the constraints. The best candidate link can thus be discarded if the constraints are violated by this addition. In this paper, we will use the 5-FLOOP-IPS algorithm, which prohibits frustrated loop with size up to 5.

In practice, -IPS alternates many link additions, corresponding to a significant increase in mean connectivity, with block update procedures, which allow to fine-tune the coefficients. Sparse precision matrices of good likelihood are generated in O(N 3 ) steps, with the advantage of having the complete optimization path available, by means of graphical models of intermediate connectivity. Note finally that we have discarded the standard way of generating sparse precision matrices based on the Lasso penalty [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Hsieh | Sparse inverse covariance matrix estimation using quadratic approximation[END_REF]. There are two reasons for that: firstly the L 1 norm penalty suffers from a modeling bias, due to excessive penalization of truly large magnitudes entries of A; secondly it is not flexible enough for the kind of constraints we are interested in, in order to produce graphical models compatible with GaBP of high likelihood [START_REF] Martin | GMRF estimation under topological and spectral constraints[END_REF].

Gaussian copula model of traffic indices

A key feature of our method is a mapping of raw data, that can correspond to flow or speed for instance, to a standard normal variable, a kind of properly normalized traffic index on which the prediction is performed. The joint probability measure of these traffic indices is approximated through a Gaussian copula. We define this index in the following way. Let t be a discretized time, measured in time steps δ t of fixed length. N t represents the number of such time steps contained in a single day. For a given t, the time of day τ ∈ {0, . . . N t -1} is given by τ = t modulo N t . For instance, we have N t = 96 if we consider δ t = 15 min time steps. At each time step, we assume the system to be represented by N v variables X t i , i ∈ {1, . . . N v } corresponding to traffic detectors. Now, for each variable X i and each time of day τ , we build from the historical data a running average Xτ i and a variance V τ i . If the dataset is clustered into a certain number of weekly or seasonal patterns representing contextual data, then these quantities will be estimated for each cluster labeled by some extra index . Then, for each variable index i, time t (and associated time of day τ ) and cluster label , we map X t, i to the following variable

U t, i def = X t, i -Xτ, i V τ, i , (12)
which represents a centered and normalized variable for given τ and . From this transformed historical data, we build for each index i a single cumulative distribution function

F i (x) = P (U i < x), (13)
which for a given realization U t, i = x represents the traffic index. The purpose of this index is to encapsulate all average time-dependent trends, week-day and seasonal dependencies, while the Gaussian copula will represent the fluctuations around these trends. In order to build a sparse Gaussian copula of all the indices, we first transform each variable U t, i into a normal centered variable via the following mapping:

Y t, i = F -1 N (0,1) • F i (U t, i ), (14)
where F N (0,1) is the cumulative distribution of a standard normal variable. The copula model corresponding to n + h + 1 time layers (n past layers, 1 present and h future layers) is then obtained by considering the vector

Z t = (Y t+k, i , i = 1 . . . N v , k = -n, . . . , 0, . . . h) (15)
and constructing its associated sparse multivariate approximate model with -IPS. This model will be used to generate predictions Ŷi , which in turn can be converted into predictions Xi of the original variables by inverting [START_REF] Manibardo | Deep learning for road traffic forecasting: Does it make a difference?[END_REF][START_REF] Li | Diffusion convolutional recurrent neural network: Data-driven traffic forecasting[END_REF].

Workflow of the method

The complete workflow of our method (see Figure 1) comprises two separate tasks: creating the model, and using it to produce detailed forecasts. We give the details of our implementation below.

Model Building

O F F L I N E Model Architecture Data Encoding (copula) Covariance Matrix

GMRF Model (*-IPS)

Figure 1: Overall workflow of the proposed methodology. The "Model Building" top part of the picture (in blue) represents the offline job for creating our model starting from an historical database. This offline part can be scheduled with a frequency of days/weeks. The "Imputation & Prediction" bottom part of the picture (in orange) is the continuous online part where a forecast is produced, by running the model on the latest real-time data.

t t + h t -np + 1 t + h + nf -1

Model building (MB, offline)

This task is supposed to be run daily, weekly or even monthly, in order to refresh the model with recent data. It typically takes several hours to complete.

(MB0) Model architecture Choose an architecture for the model, that is:

• n p : number of past and present layers of variables for which observations may be available; if t is the present time, this covers times tn p + 1, . . . , t.

• n f : number of future layers we want to predict, n f ≥ 1;

• h: horizon of prediction of the first future layer.

From this, we construct by concatenation of (n p + n f )N v variables the dynamical configurations (see Figure 2):

{X t-np+1 i , . . . , X t i , X t+h i , . . . , X t+h+n f -1 i , i = 1, . . . , N v }. ( 16
)
This subtask is only supposed to be done once, when setting up the original model. There is no need to tweak the values later when updating historical data.

(MB1) Data encoding (copula) Compute the required statistics of the model from the training set: the mean and variance Xτ, i and V τ, i for each variable i, time of day τ and label for the day (when a clustering is in use). Then, for each i, compute the c.d.f. F i of U t, i aggregated over all t and approximate it by an invertible monotonous linear piecewise function.

(MB2) Covariance matrix Generate using [START_REF] Manibardo | Deep learning for road traffic forecasting: Does it make a difference?[END_REF] a training set of vectors

{Y t-np+1 i , . . . , Y t i , Y t+h i , . . . , Y t+h+n f -1 i , i = 1, . . . N v } (17)
and compute the corresponding

(n p + n f )N v × (n p + n f )N v covariance matrix.
Because the data is incomplete, some regularization can be required if some modes are associated to large negative eigenvalues. In such a case the eigenvalue is replaced by its absolute value.

(MB3) GMRF model ( -IPS) The model is built using the -IPS algorithm [START_REF] Martin | GMRF estimation under topological and spectral constraints[END_REF], which takes the regularized covariance matrix as input and generates an almost continuous set of models with increasing mean connectivity d. The variant that we use (5-FLOOP-IPS) avoids frustrated short loops of sizes up to 5 (see Section 3.2).

Imputation and prediction (IP, online)

At this point, the model is ready to use for prediction tasks. It is defined by a multivariate Gaussian distribution P (Y ) in which the variables corresponding to observations are clamped, while the other (future and non-observed past variables) are inferred using GaBP, as described in Section 3.1. We insist also on the fact that the same model is used independently of day-time or any other contextual data (see Section 3.3).

(IP1) Real-Time data ingestion Initialize the GMRF with observations, using strong beliefs: if no observation exists for node i, then h i = 0 by construction (see Section 3.3); otherwise, x i is fixed to its observed value xi and detached from the factor graph and neighboring variables see their local field modified:

h j ← h j -A ij xi , for all j ∈ ∂i (18)
Another possibility not used here is to use soft beliefs: the observable xi comes with an uncertainty either given by some variance σi .

(IP2) Inference (GaBP) Run GaBP according to equations ( 7)-( 8) until convergence is reached. The predictions are then obtained from [START_REF] Lv | Traffic flow prediction with big data: A deep learning approach[END_REF] for each variable to be predicted, along with an estimation of the uncertainty provided by the variance [START_REF] Yu | Spatiotemporal recurrent convolutional networks for traffic prediction in transportation networks[END_REF]. Note that, for observed variables inserted with soft constraints, the value of the variance [START_REF] Yu | Spatiotemporal recurrent convolutional networks for traffic prediction in transportation networks[END_REF] provided by the fixed point may give some indication of whether the input observation is trustworthy or not.

Practical considerations

There are a couple of practical issues to be aware of when using this workflow.

Model tuning There are basically two hyper-parameters of the model to tune, namely the number of past layers n p and the mean connectivity d * . Indeed, h is imposed by the task and n f is an optional choice, defaulting to 1. d * can be loosely selected somewhere in the region where the log likelihood monitored by -IPS as a function of d becomes flat. A more systematic procedure, leading to similar parameters, relies on a grid search, based on accuracy of predictions measured on the training set. At low connectivity, we observe similar pattern in all datasets: a sharp rise in accuracy, followed by an almost flat plateau while increasing both d and n p . Their values are therefore chosen at the beginning of the flat domain in order to get a good trade-off between speed and precision with a loss in accuracy below 1% w.r.t. the would-be optimal setting while computational time is reduced by a significant factor.

Complexity and affordable systems sizes Since -IPS has a cubic complexity, the model building task (MB) has a running time that scales like

(n p + n f )N v 3 .
In practice, the total number of variables

(n p + n f )N v should
not exceed the order of 10 4 so as to be able to train the model in a reasonable time, which means that, if we consider n p = 2 and n f = 1, we can deal with systems of roughly 3000 detectors.

Concerning the imputation and prediction task (IP), GaBP can converge on networks of size 10 5 within a few tenth of a second on an ordinary CPU, while typically, the online constraint for traffic information is to generate a forecast every few minutes, according to the refresh rate of the data.

Experiments

In order to check the robustness of the method, we have tested our system on two different datasets corresponding to urban traffic, both for imputation and prediction. Horizon of prediction range from 15 minutes to 5 hours, with input observation covering up to one past hour. The datasets considered here are the cities of Vienna and Turin. More extensive results, including figures for the Turin dataset and an additional dataset about highways in the San Francisco Bay Area, are contained in the research report [START_REF] Furtlehner | Spatio-temporal probabilistic short-term forecasting on urban networks[END_REF].

Error metrics

We consider a set of metrics commonly used in traffic analysis. General purpose metrics are the root mean square error (RMSE), the mean absolute error (MAE), and the mean absolute percentage error (MAPE). Letting X be the variable to predict and X its predictor, the metrics are defined as follows.

RMSE def = E[( X -X) 2 ] ( 19
) MAE def = E[| X -X|] (20)
MAPE def = 100 × E | X -X| X (21)
Besides, for traffic flows, we use the GEH statistic, introduced in the 1970s by Geoffrey E. Havers, and used today as a standard indicator by Transport for London, one of the leading transport authorities worldwide [START_REF] Smith | Traffic Modeling Guidelines Version 3.0[END_REF]. When X is an hourly traffic flow, GEH is defined as

GEH def = 2(X -X) 2 X + X . (22)
This number has been empirically designed to yield comparable numbers for a broad range of arterial road capacities. Typically, the prediction is considered good if GEH is below 5, and the statistic of interest is thus the percentage of GEH < 5.

The expectations above are taken uniformly over all detectors and all time steps for which there is an observation to compare against the prediction. Note that MAPE is not defined when X = 0, which can happen in particular for flow measurements. To regularize this, we simply threshold the denominator in MAPE to X min = 10 in practice.

Baseline predictors

In order to provide some reference points of comparison for our results we use some simple predictors and when possible, i.e. for the Vienna dataset, other results from the literature. Our simple baseline predictors are the following:

• Mean(t) is the historical average Xτ, i for detector i at time of day τ and label , as defined in Section 3.3.

• t 0 is the persistent predictor, sometimes called the random walk predictor.

It is equal to the last observed value within the past time window of the detector, and to Mean(t) if no such value is found.

• k-NN is a k nearest neighbors predictor. Given the past and present observations, which can comprise many time layers of incomplete data, we look in the training set for the k closest states to the current ones by means of an L 2 discrepancy. The hyper-parameter k needs to be optimized.

For the Vienna dataset, other methods, based respectively on Bayesian networks (BN), neural networks (NN) and clustering of day profiles, have been previously tested [START_REF] Attanasi | A hybrid method for real-time short-term predictions of traffic flows in urban areas[END_REF] on only 36 count locations. A direct comparison (see Table 2) shows that GaBP improves significantly on these techniques, especially for longer forecast horizons.

Vienna dataset

Our first dataset concerns urban traffic in the Vienna (Austria) urban area, as depicted in Figure 3. It consists in flow data collected over 4 years (2011-2014) from 292 count locations (CLOC) measuring the number of vehicles every 15 minutes. Since these detectors are quite often down, we arbitrarily select those which are up during at least 10% of the time. This leads to retaining 266 CLOCs q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q CLOC status q q discarded used out of 292. With this setting, at any given time, we have typically around 65% of the detectors which are up. This means that, in addition to forecasting a certain number of layers of variables ahead in time, our model performs an imputation of at least 35% of missing observations (more when we randomly hide observations). The first 3 years (2011-2013) are used as a training set, while 2014 is used for testing. As a preprocessing step, we have performed a clustering analysis of the daily traffic conditions using the Affinity Propagation algorithm [START_REF] Frey | Clustering by passing messages between data points[END_REF]. This analysis yields 10 meaningful clusters, presented in Table 1. This clustering will be used in two ways: first it will serve to build a baseline predictor X(τ, ) i for each CLOC i, time-step τ = 1, . . . , 96 and cluster index = 1, . . . , 10; secondly our model, as explained before, will be built on the residues U (t, ) i defined in [START_REF] Li | Diffusion convolutional recurrent neural network: Data-driven traffic forecasting[END_REF] from this baseline or some related ones, obtained in the same way after merging some of the smallest clusters into larger ones, in order to have a sufficient number of observations in each cluster to avoid statistical problems when building the model.

For this setting, several architectures lead to similar performances, as long as at least 3 layers in the past are taken into account, corresponding to a look back interval of at least 30 minutes. Specializing the model to one single horizon instead of directly performing the forecast for many features ahead does not seem to help much. The results presented in Figure 4 were obtained with a model with 4 past layers and one single specialized future layer. The results summarized in Table 2 were obtained with a multistep-ahead model having 4 past and 4 future layers. Specialized or multistep-ahead models yield identical performance within error bars. The multistep-ahead model appears to be more advantageous in practice, as it can deliver all horizon predictions in one pass. The only drawback Jan. Feb. Mar. Apr. May Jun. Jul. Aug. Sep. Oct. Nov. Dec. Sun. 1 comes from the cubic increase in computational cost necessary to build it, since it contains more variables.

As our model is able to perform both completion of missing data and forecasting at the same time within a single GaBP fixed point convergence, we first study the dependency of the forecasting performance on the fraction of observed variables. To this end, we randomly ignore some of the available observations in the past and present layers, in order to have a given density ρ of observed variables per layer, as illustrated in Figure 2. As shown in Figure 4 (top), a lot of precision is gained by observing only a small fraction of variables, say 5-10%, and after that the forecasting and imputation errors continue to decrease more slowly but steadily. Conversely, the prediction error of the k-NN model saturates after ρ = 0.2. The forecasting performance of the model as a function of the horizon is shown in Figure 4 (bottom). We see that the GaBP error is 33% smaller than that of Mean(t) on the short-term prediction. In addition, our model provides meaningful information up to 5 hours in advance w.r.t. this baseline, which means that it is able to identify additional traffic patterns not captured by the clustering.

Table 2 provides the numerical results for the 266 CLOCs, as well as comparison with the predictors proposed in [START_REF] Attanasi | A hybrid method for real-time short-term predictions of traffic flows in urban areas[END_REF] tested on 36 of these CLOCs. These methods treat each detector independently, so the comparison gives an indication of the added value of taking spatial dependencies into account. In all cases, we see some gain at all horizons regarding all metrics of evaluation, rather marginal for h = 15 but rapidly sizable when the horizon increases. Visual inspection of individual CLOC errors shows that those having less than 70% of GEH < 5 are either badly behaving sensors clearly sending corrupted data, or sensors for which a systematic bias is present in our Mean(t) baseline, possibly reflecting a long lasting modification in the traffic conditions for the corresponding segment.

Flow (266 CLOCs)

Flow (same 36 CLOCs as [START_REF] Attanasi | A hybrid method for real-time short-term predictions of traffic flows in urban areas[END_REF]) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q qq q q q q q q q q q q q q q q q q qq q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q qq q q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q qq q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q qq q q qq qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q qq qq q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q qq q qq q qq q q qq q q q q q q q qq q q q qq q q q q q q q q q q qq q q q q qq q q q q q q q q q qqq qq q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q qq q qq q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq qq q q q CLOC status q q discarded used Figure 5: Location of the CLOCs in the center of Turin. Only the blue ones have enough data to be usable. Note that there are 148 additional sensors outside of this area that are not displayed here for the sake of clarity.

%GEH < 5 RMSE MAE MAPE %GEH < 5 RMSE MAE MAPE h = 15 BN - - - - - 33

Turin dataset

The second dataset again concerns urban traffic, in the Turin (Italy) urban area. There are 1489 count locations, giving average speed measurements in km/h in addition to the flow measurements as number of vehicles every 15 minutes. The map of these CLOCs in the center of Turin is shown in Figure 5. Two years of data (2015-2016) have been collected. As for the previous dataset, we select detectors based on an activity threshold of 10%. We end up with 685 detectors, of which 566 correspond to flow and 119 to speed measurements. With this setting, the ratio of available data is around 60% in the past layers. Year 2015 is used for training and year 2016 for testing. A straightforward clustering is done here with 4 clusters: one for Monday-Thursday and the other ones respectively for Fridays, Saturdays and Sundays. This leads to a slightly noisier Mean(t) baseline (around 45 in RMSE instead of 40 for Vienna). Table 3 compares different methods to a multistep-ahead GaBP with (n p , n f ) = (4, 4) with 5480 variables. The selected model is sparser than for Vienna, with mean connectivity 4. Yet the behavior and performance are comparable, with e.g. more than 85% of GEH < 5 at a 1 hour horizon. All indicators for the flow are very similar to those of the Vienna dataset, with slightly higher RMSE and lower MAPE and MAE, because of an higher overall capacity bias of the segments in the case of Turin. As far as speed is concerned, we observe for instance at a 30-minute horizon an RMSE improvement of 10% w.r. the aggregated error for the speed is dominated by nighttime prediction errors, where the small number of speed measurements leads to a very noisy signal.

Finally, we examine a subset of the results from a transportation engineering point of view. We consider the set of 22 CLOCs from the dataset for which both flow and speed are available. For each of these CLOCs, the object of interest is the empirical fundamental diagram, i.e. the flow versus the density (obtained as a simple division between the flow and the speed) of all the data of the year 2016. Then, each data point is labeled as describing either a fluid (hypocritical) or a congested (hypercritical) traffic condition. The manual labeling of the traffic conditions has been performed by three different experts, all of them trained to recognize visually the traffic conditions. These three classifications are not 100% identical, since some boundary situations are difficult to classify. In those cases the majority label is used, which is always possible with an odd number of outcomes. The main use for this manual labeling is to assess the expectation that congested traffic conditions are harder to predict than fluid ones. To this end, we measure the performance of the flow forecast against actual data through the GEH indicator, evaluated either on the congested states, the fluid ones or all of them, as shown in

Error analysis and limits of the model

In this section, in order to estimate the margin and possible directions for improvement, we first discuss the main limit of accuracy of our model and then perform a data analysis of the errors made by the model. In particular, we would like to be able to separate components in the errors due to modeling bias from the intrinsic noise of the traffic phenomenon, which cannot be predicted.

Up to now, we have left aside the fact that our prediction µ i based on GaBP in [START_REF] Lv | Traffic flow prediction with big data: A deep learning approach[END_REF] comes with an uncertainty estimate through the variance σ i [START_REF] Yu | Spatiotemporal recurrent convolutional networks for traffic prediction in transportation networks[END_REF]. While µ i is guaranteed to be the exact marginal of the model once GaBP has converged, σ i on the other hand is only an approximate value of the true variance. These values can be exploited to deliver levels of confidence on our predictions. Applying the inverse map of ( 14) to the corresponding confidence interval, defined by say ±1 standard deviation in copula space, we get straightforwardly a confidence interval of our prediction in the prediction space, e.g. flow, speed or occupancy. This confidence interval proves to be very reliable: about 67% of the predictions lay inside it, indifferently for speed or flow in Vienna or Turin, which is very close to the theoretical value of 68.27%, corresponding to one standard deviation in copula space. As shown, for example, on Figure 6, these confidence intervals are quite consistent and meaningful and may help to identify detector errors. Indeed, we see that the detector wakes up most probably in the middle of a time interval after being silent during 6 time steps. Then it delivers a clearly underestimated observation of traffic flow, 2 standard deviations away from our confidence interval. These confidence intervals reflect an estimation of some intrinsic noise of the traffic signal, which might not be possibly predicted. In any case, it is comparable to the mean error that we finally end up measuring on our predictions. That said, there are systematic errors that our model makes which we would like to identify and cure. The main source of error comes from the Gaussian copula hypothesis. Recall first that, after the transformation ( 14) is performed, each Y t, i taken individually is by construction and up to numerical precision, a standard normal variable. However, the joint distribution has no specific reason in general to be multivariate Gaussian. In order to estimate how far from a multivariate Gaussian our model is, we consider the main directions of fluctuations of Y by extracting the dominant eigenmodes of the covariance matrix. Then for each of these modes e .,α , we compute the corresponding cumulative distribution P (z t, α /σ α < x) of the normalized components defined as

z t, α = N i=1 e i,α Y t, i , σ α def = E (z t, α ) 2 . ( 23
)
Figure 7 shows for the Vienna dataset how these distributions compare to the expected cumulative distribution of a standard normal variable. As we can see the alignment is pretty good for most of the dominant modes over more or less 2 standard deviations. Beyond that, the model shows inadequacy in the distributions tails which are clearly non-Gaussian. Tests done without the clustering preprocessing step show as expected degraded performance, with a much more pronounced deviation from the normal distribution and a multimodal behavior. This underlines the importance of the clustering of the data beforehand.

Another way to look at the results is to consider the Fourier spectrum of the error signal, and in particular how the error is distributed along the different frequencies. Figure 8 shows this for the different datasets and types of variables, compared to a white noise predictor, that is a with random decorrelated errors. First, we notice that the departure from the white noise predictor is more pronounced for flow than for speed. What is also noticeable is the finite contribution of very low frequencies. When looking more precisely at the data, it turns out that this corresponds to detectors which, for a long period, send values completely at odds with the ones observed during training. These badly behaving detectors may either correspond to corrupted ones, or to drastic changes of the traffic conditions on the corresponding segment, because of road work for instance. Then there remains the domain corresponding to time scales ranging between roughly 30 minutes and one day, represented by the right part of the plot. Overall, as seen in Figure 8, this corresponds to more than 80% of the total error. It is in this region that we can mostly pinpoint modeling biases.

Conclusions

The method presented in this paper is intended both to provide valuable predictors of traffic variables on large traffic networks, in order to feed traffic management systems for instance, and to provide some understanding of the statistical properties of traffic. With the advent of big data in the domain of traffic management, it is becoming possible to analyze these properties in detail at the macroscopic level. The structure of the generated network could potentially help us to visualize the information flow among the variables. Various sources of improvement can also be expected, like e.g. working with variables corresponding to positions in the fundamental diagrams at the level of each segment. The question of corrupted data was merely touched upon in this paper. Our model offers the possibility, via the message passing structure, of estimating the reliability of an input variable, depending on the information coming from other variables. Integrating this information in the belief propagation schema with the help of the soft belief setting mentioned in Section 4.2 might help us to improve the accuracy of the forecast by discarding suspicious observations.

Figure 2 :

 2 Figure 2: Time layers setting of one basic space time configuration. Each circle represents one single variable at a given time, filled circles correspond to observations.

Figure 3 :

 3 Figure 3: Location of the CLOCs in Vienna. Only the blue ones have enough data to be usable.

Figure 4 :

 4 Figure 4: Vienna dataset. Top: average flow forecasting error for a fixed horizon of 30 minutes, as a function of the fraction ρ of observed variables in the past time layers, averaged over 5000 test samples. The imputation error for missing data is also shown (in blue). A point of comparison is given by the k-NN predictor, optimized with k = 50. Bottom: Average flow forecasting error as a function of forecasting horizon at maximum possible observation rate ρ ∼ 0.65 (average over the full test set).

Figure 6 :

 6 Figure6: Vienna dataset: flow for CLOC "MA_33_23039_R3", together with uncertainty estimate delivered by GaBP. The confidence interval is 1 standard deviation in copula space, corresponding to 68.27% of confidence that the true value is within this interval.

Figure 7 :

 7 Figure 7: Vienna dataset: empirical cumulative distribution of the normalized projection of the (copula transformed) data along the 50 first principal modes against the cumulative distribution of a standard normal variable. Modes are ordered by decreasing eigenvalues λ α .

Figure 8 :

 8 Figure 8: Cumulative power spectrum of the mean errors signal of the GaBP predictor. The vertical line indicates the 1-day time scale. The black dotted line represents the white noise (WN) reference.

Table 1 :

 1 Public holidays are classified as Sundays; the ranges for Saturday and Sunday are by season, so that Cluster 5, for example, starts on March 21 st . Description of the 10 clusters for Vienna dataset according to month and day of the week.

		9	9	9	5	5	5	5	5	5	9	9	9
	Mon. 1	1	1	4	4	4	7	7	8	8	10 10
	Tue.	1	1	1	4	4	4	7	7	8	8	10 10
	Wed. 1	1	1	4	4	4	7	7	8	8	10 10
	Thu. 1	1	1	4	4	4	7	7	8	8	10 10
	Fri.	3	3	3	3	3	3	3	3	3	3	3	3
	Sat.	2	2	2	2	2	2	6	6	6	2	2	2
	1

Table 2 :

 2 Vienna dataset: results for different horizons of forecasting and comparison on a subset of CLOCs with single time-series forecasting based methods (BN and NN from[START_REF] Attanasi | A hybrid method for real-time short-term predictions of traffic flows in urban areas[END_REF]).

Table 3 :

 3 t. the Mean(t) baseline and 20% w.r.t. the t 0 predictor. The results for flow are clearly more impressive than those for speed predictions. This hides important disparities among various days. Moreover, Turin dataset: results for different horizons of forecasting, for both flow and speed measurements.

			Flow (566 CLOCs)		Speed (119 CLOCs)
		%GEH < 5 RMSE MAE MAPE RMSE MAE MAPE
	h = 15
	Mean(t)	75.46	45.09	21.65	31.02	8.40	5.21	13.73
	t 0	81.41	27.58	15.43	20.32	9.24	5.00	11.74
	k-NN	79.78	37.71	18.47	27.19	8.20	4.99	13.28
	GaBP	88.10	27.56 13.32 19.03	7.65	4.24	10.44
	h = 30
	Mean(t)	75.46	45.09	21.65	31.02	8.40	5.21	13.73
	t 0	73.67	35.02	19.22	24.50	9.68	5.35	12.70
	k-NN	77.65	39.70	19.55	29.42	8.34	5.11	13.78
	GaBP	87.10	30.73 14.20 20.68	7.56	4.25	10.70
	h = 60
	Mean(t)	75.46	45.09	21.65	31.02	8.40	5.21	13.73
	t 0	62.00	47.92	26.67	32.22	10.36	5.88	14.12
	k-NN	73.99	42.28	21.32	31.51	8.43	5.15	13.32
	GaBP	85.24	33.90 15.45 22.58	7.67	4.34	11.09

Table 4 .

 4 As expected, the congested states are more difficult to predict for all methods. Plotting the cumulative distributions of GEH for different forecast horizons shows that GaBP achieves consistently a GEH score that is smaller in distribution than the other methods in all traffic state %GEH < 5 All Data Fluid Congested

	h = 15
	Mean(t)	70.93	71.11	63.12
	t 0	79.18	79.37	72.03
	GaBP	88.28	88.36	85.32
	h = 30
	Mean(t)	70.80	70.96	62.93
	t 0	66.71	66.83	61.08
	GaBP	85.71	85.79	81.43
	h = 45
	Mean(t)	71.09	71.26	62.99
	t 0	57.47	57.50	55.73
	GaBP	84.81	84.95	79.90
	h = 60
	Mean(t)	70.94	71.13	62.61
	t 0	50.34	50.28	53.08
	GaBP	83.14	83.28	77.30

Table 4 :

 4 Turin dataset: results with manually labeled traffic state conditions.

conditions.

Available at https://gitlab.inria.fr/bptraffic/star-ips.

Acknowledgments

The authors would like to thank the municipality of Vienna (Austria) and the Italian company 5T, which manages the mobility in Piemonte region (including Turin), for sharing their data with them.