{'fraction_non_alphanumeric': 0.04345350103991284, 'fraction_numerical': 0.019857383381202336, 'mean_word_length': 4.593269630245119, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
RNN with Particle Flow for Probabilistic Spatio-temporal Forecasting

Spatio-temporal forecasting has numerous applications in analyzing wireless, traffic, and financial networks. Many classical statistical models often fall short in handling the complexity and high non-linearity present in time-series data. Recent advances in deep learning allow for better modelling of spatial and temporal dependencies. While most of these models focus on obtaining accurate point forecasts, they do not characterize the prediction uncertainty. In this work, we consider the time-series data as a random realization from a nonlinear state-space model and target Bayesian inference of the hidden states for probabilistic forecasting. We use particle flow as the tool for approximating the posterior distribution of the states, as it is shown to be highly effective in complex, high-dimensional settings. Thorough experimentation on several real world time-series datasets demonstrates that our approach provides better characterization of uncertainty while maintaining comparable accuracy to the state-of-the art point forecasting methods.


Introduction
Spatio-temporal forecasting has many applications in intelligent traffic management, computational biology and finance, wireless networks and demand forecasting. Inspired by the surge of novel learning methods for graph structured data, many deep learning based spatio-temporal forecasting techniques have been proposed recently (Li et al., 2018;Bai et al., 2020). In addition to the temporal patterns present in the data, these approaches can effectively learn and exploit spatial relationships among the time-series using various Graph Neural Networks (GNNs) (Defferrard et al., 2016;Kipf & Welling, 2017). Recent works establish that graph-based spatio-temporal models outperform the graphagnostic baselines (Li et al., 2018;Wu et al., 2020). In spite of their accuracy in providing point forecasts, these models have a serious drawback as they cannot gauge the uncertainty in their predictions. When decisions are made based on forecasts, the availability of a confidence or prediction interval can be vital.
There are numerous probabilistic forecasting techniques for multivariate time-series, for example, DeepAR (Salinas et al., 2020), DeepState (Rangapuram et al., 2018), Deep-Factors , and the normalizing flow-based algorithms in Rasul et al., 2021). Although these algorithms can characterize uncertainty via confidence intervals, they are not designed to incorporate side-knowledge provided in the form of a graph.
In this work, we model multivariate time-series as random realizations from a nonlinear state-space model, and target Bayesian inference of the hidden states for probabilistic forecasting. The general framework we propose can be applied to univariate or multivariate forecasting problems, can incorporate additional covariates, can process an observed graph, and can be combined with data-adaptive graph learning procedures. For the concrete example algorithm deployed in experiments, we build the dynamics of the state-space model using graph convolutional recurrent architectures. We develop an inference procedure that employs particle flow, an alternative to particle filters, that can conduct more effective inference for high-dimensional states.
The novel contributions in this paper are as follows: 1) we propose a graph-aware stochastic recurrent network architecture and inference procedure that combine graph convolutional learning, a probabilistic state-space model, and particle flow; 2) we demonstrate via experiments on graph-based traffic datasets that a specific instantiation of the proposed framework can provide point forecasts that are as accurate as the state-of-the-art deep learning based spatio-temporal models. The prediction error is also comparable to the existing deep learning based techniques for benchmark non-graph multivariate time-series datasets; arXiv:2106.06064v1 [stat.ML] 10 Jun 2021 3) we show that the proposed method provides a superior characterization of the prediction uncertainty compared to existing probabilistic multivariate time-series forecasting methods, both for datasets where a graph is available and for settings where no graph is available.

Related Work
Our work is related to (i) multivariate and spatio-temporal forecasting using deep learning and graph neural networks; (ii) stochastic/probabilistic modelling, prediction and forecasting for multivariate time-series; and (iii) neural (ordinary) differential equations. Recently, neural network-based techniques have started to offer the best predictive performance for multivariate time-series prediction (Bao et al., 2017;Qin et al., 2017;Lai et al., 2018;Guo & Lin, 2018;Chang et al., 2018;Li et al., 2019;Sen et al., 2019;Oreshkin et al., 2020;Smyl, 2020). In some settings, a graph is available that specifies spatial or causal relationships between the time-series. Numerous algorithms have been proposed that combine GNNs with temporal neural network architectures (Li et al., 2018;Yu et al., 2018;Huang et al., 2019;Bai et al., 2019;Chen et al., 2019a;Guo et al., 2019;Yu et al., 2019;Zhao et al., 2019;Bai et al., 2020;Huang et al., 2020;Park et al., 2020;Shi et al., 2020;Song et al., 2020;Xu et al., 2020;Wu et al., 2020;Zheng et al., 2020;Oreshkin et al., 2021). Algorithms that take into account the graph provide superior forecasts, if the graph is accurate and the indicated relationships have predictive power. However, none of these algorithms is capable of characterizing the uncertainty of the provided predictions; all are constructed as deterministic algorithms.
Our proposed method is different from this body of work in two important ways. First, we design a probabilistic statespace modeling framework that can incorporate information about predictive relationships that is provided in the form of a graph. Second, our inference procedure employs particle flow, which avoids the need for some of the approximations required by a variational inference framework and is much better suited to high-dimensional states than particle filtering. Our particle flow method has connections to normalizing flows (Kobyzev et al., 2020) and to neural ordinary differential equations (Chen et al., 2018;2019b). In particular, Chen et al. (2019b) address a Bayesian inference task by solving a differential equation to transport particles from the prior to the posterior distribution. However, such flow-based methods were first introduced by Daum & Huang (2007) in the sequential inference research literature. The filtering task is to approximate the posterior distribution of the state trajectory p Θ (x t |y 1:t ) recursively, where x t denotes the state at time t and y 1:t are observations from times 1 to t. A particle filter (Gordon et al., 1993;Doucet & Johansen, 2009) maintains a population of N p samples (particles) and associated weights {x j t , w j t } Np j=1 that it uses to approximate the marginal posterior distribution of x t : Here, δ(·) denotes the Dirac-delta function. Particles are propagated by the application of importance sampling using a proposal distribution; the weights are updated accordingly. When the disparity in the weights becomes too great, resampling is applied, with particles being sampled proportionally to their weights and the weights being reset to 1. Constructing well-matched proposal distributions to the posterior distribution in high-dimensional state-spaces is extremely challenging. A mismatch between the proposal and the posterior leads to weight degeneracy after resampling, which results in poor performance of particle filters in high-dimensional problems (Bengtsson et al., 2008;Snyder et al., 2008;Beskos et al., 2014). Instead of sampling, particle flow filters offer a significantly better solution by transporting particles continuously from the prior to the posterior (Daum & Huang, 2007;Ding & Coates, 2012).
For a given time step t, particle flow algorithms solve differential equations to gradually migrate particles from the predictive distribution so that they represent the posterior distribution after the flow. A particle flow can be modelled by a background stochastic process η λ in a pseudo-time interval λ ∈ [0, 1], such that the distribution of η 0 is the predictive distribution p Θ (x t |y 1:t−1 ) and the distribution of η 1 is the posterior distribution p Θ (x t |y 1:t ).
One approach (Daum et al., 2010), is to use an ordinary differential equation (ODE) with zero diffusion to govern the flow of η λ : For linear Gaussian state-space models, the flow can be expressed in the form: and we can derive analytical expressions for A(λ) and b(λ) (see supplementary material for details). For non-linear and non-Gaussian models, we employ Gaussian approximations and repeated local linearizations.

State-space model
We postulate that y t ∈ R N ×1 is the observation from a Markovian state space model with hidden state X t ∈ R N ×dx . We denote by x t and z t the vectorizations of X t and Z t , respectively. The state space model is: Here v t ∼ p v (·|x t−1 , σ) and w t ∼ p w (·|x t , γ) are the noises in the dynamic and measurement models respectively. ρ, σ and γ are the parameters of the distribution of the initial state x 1 , process noise v t and measurement noise w t respectively. g and h denote the state transition and measurement functions, possibly linear or nonlinear, with parameters ψ and φ respectively. The subscript G in g and h indicates that the functions are potentially dependent on the graph topology. We assume that h G,φ (x t , z t , 0) is a C 1 function in x t , i.e., h G,φ (x t , z t , 0) is a differentiable function whose first derivative w.r.t. x t is continuous. The complete set of the unknown parameters is formed as: Θ = {ρ, ψ, σ, φ, γ}. Figure 2 depicts the graphical model relating the observed variables (y t and z t ) to the latent variables (v t , w t ) and the graph (G). With the proposed formulation, we can modify recurrent graph convolutional architectures when designing the function g. When a meaningful graph is available, such architectures significantly outperform models that ignore the graph. For example, we conduct experiments by incorporating into our general model the Adaptive Graph Convolutional Gated Recurrent Units (AGCGRU) presented in (Bai et al., 2020). The AGCGRU combines (i) a module that adapts the provided graph based on observed data, (ii) graph convolution to capture spatial relations, and (iii) a GRU to capture evolution in time. The example model used for experiments thus employs an L-layer AGCRU with additive Gaussian noise to model the system dynamics g: In this model, we have p v (v t ) = N (0, σ 2 I), i.e., the latent variables for the dynamics are independent. The initial state distribution is also chosen to be isotropic Gaussian, i.e., p 1 (x 1 , z 1 , ρ) = N (0, ρ 2 I). The parameters ρ and σ are learnable variance parameters. The observation model g incorporates a linear projection matrix W φ . The latent variable w t for the emission model is modelled as Gaussian with variance dependent on x t via a learnable softplus function:

Inference
We assume that a dataset D trn is available for training. Although this data may be derived from a single time-series, because our task is to predict y t0+P +1:t0+P +Q using a limited historical window y t0+1:t0+P , we splice the time-series and thus construct multiple training examples, denoted by (y (m) 1:P , y (m) P +1:P +Q ). In the training set, all of these observations are available; in the test set y P +1:P +Q are not. In addition, the associated covariates z 1:P +Q are known for both training and test sets.
Inference involves an iterative process. We randomly initialize the parameters of the model to obtain Θ 0 . Subsequently, at the k-th iteration of the algorithm (processing the k-th training batch), we first draw samples {y i P +1:P +Q } Np i=1 from the distribution p Θ k−1 (y P +1:P +Q |y 1:P , z 1:P +Q ). With this set of samples, we can subsequently apply a gradient descent procedure to obtain the updated model parameters Θ k . We discuss each of these steps in turn as follows.
5.2.1. SAMPLING In a Bayesian setting with known model parameters Θ = Θ k−1 , we would aim to form a prediction by approximating the posterior distribution of the forecasts, i.e., p Θ (y P +1:P +Q |y 1:P ). (For conciseness we drop the timeoffset t 0 ). p Θ (y P +1:P +Q |y 1:P , z 1: Since the integral in eq. (10) is analytically intractable for a general nonlinear state-space model, we take a Monte Carlo approach as follows: Step 1: For 1 t P , we apply a particle flow algorithm (details in Sec. 4) with N p particles for the state-space model specified by eqs. (4), (5) and (6) to recursively approximate the posterior distribution of the states: Here {x j t } Np j=1 are approximately distributed according to the posterior distribution of x t . The generation of each sample x j t involves an associated sampling of the latent variables v 1:t and implies a sampling of w 1:t , but these samples are not required since the proposed model only needs x P to construct the forecast, Step 2: For P + 1 t P + Q, we iterate between the following two steps: a. We sample x j t from p ψ,σ (x t |x j t−1 , y j t−1 , z t ) (for t > P + 1) or from p ψ,σ (x t |x j t−1 , y t−1 , z t ) (for t = P + 1) for 1 j N p . This amounts to a state transition at time t to obtain the current state x j t from the previous state x j t−1 , using eq. (5). b. We sample y j t from p φ,γ (y t |x j t , z t ) for 1 j N p , i.e., we use x j t in the measurement model, specified by eq. (6), to sample y j t .
A Monte Carlo (MC) approximation of the integral in eq. (10) is then formed as: p Θ (y P +1:P +Q |y 1:P , z 1:P +Q ) ≈ Each y j P +1:P +Q is approximately distributed according to the joint posterior distribution of y P +1:P +Q . The resulting algorithm is summarized in Algorithm 1. A block diagram of the probabilistic forecasting procedure is shown in Figure 1.
In each training mini-batch, for each training example, we perform a forward pass through the model using Algorithm 1 to obtain approximate forecast posteriors and then update all the model parameters using SGD via backpropagation.

Experiments
We perform experiments on four graph-based and four nongraph based public datasets to evaluate proposed methods.

Datasets
We evaluate our proposed algorithm on four publicly available traffic datasets, namely PeMSD3, PeMSD4, PeMSD7 and PeMSD8. These are obtained from the Caltrans Performance Measurement System (PeMS) (Chen et al., 2000) and have been used in multiple previous works Guo et al., 2019;Song et al., 2020;Bai et al., 2020;Huang et al., 2020). Each of these datasets consists of the traffic speed records, collected from loop detectors, and aggregated over 5 minute intervals, resulting in 288 data points per detector per day. In non-graph setting, we use Electricity (Dua

Preprocessing
For the PeMS datasets, missing values are filled by the last known value in the same series. The training, validation and test split is set at 70/10/20% chronologically and standard normalization of the data is used as in (Li et al., 2018). We use one hour of historical data (P = 12) to predict the traffic for the next hour (Q = 12). Graphs associated with the datasets are constructed using the procedure in (Huang et al., 2020).

Baselines
To demonstrate the effectiveness of our model, we compare to the following forecasting methods. A detailed description of each baseline is provided in the supplementary material.
The detailed comparison of our approach with all of these models is provided in the supplementary material for space constraints. Here, we show results of a subset, focusing on those with the most competitive performance ( Figure 3).

Hyperparameters and training setup
For our model, we use an L = 2 layer AGCGRU (Bai et al., 2020) as the state-transition function. The dimension of the learnable node embedding is d e = 10, and the number of RNN units is d x = 64. We treat ρ and σ as fixed hyperparameters and set ρ = 1 and σ = 0 (no process noise). We train for 100 epochs using the Adam optimizer, with a batch size of 64. The initial learning rate is set to 0.01 and we follow a decaying schedule as in (Li et al., 2018). Hyperparameters associated with scheduled sampling , gradient clipping, and early stoppng are borrowed from (Li et al., 2018). We set the number of particles N p = 1 during training and N p = 10 for validation and testing. The number of exponentially spaced discrete steps (Li & Coates, 2017) for integrating the flow is N λ = 29. For each dataset, we conduct two separate experiments minimizing the training MAE (results are used to report MAE, MAPE, RMSE, and P50QL) and the training negative log posterior probability (results are used to report CRPS, P10QL, and P90QL). We also experiment with alternative state transition functions, including the DCGRU (Li et al., 2018) and GRU (Chung et al., 2014). For these, the hyperparameters are fixed to the same values as presented above.

Results and Discussion
Comparison with baselines : Results for the point forecasting task are summarized in Table 1. We observe that most of the spatio-temporal models perform better than graph agnostic baselines in most cases. Moreover, the pro-Some of the recent spatio-temporal models such as (Chen et al., 2020;Zhang et al., 2020;Park et al., 2020) do not have publicly available code. Although the codes for (Wu et al., 2020;Song et al., 2020; are available, these works use different datasets for evaluation. We could not obtain sensible results from these models for our datasets, even with considerable hyperparameter tuning. The code for de Bézenac et al., 2020) is not publicly available.
posed AGCGRU+flow algorithm achieves on par or better performance with the best-performing spatio-temporal models, such as GWN, GMAN and AGCRN. We present a comparison of the average rankings across datasets in Figure 3. Our proposed method achieves the best average ranking and significantly outperforms the baseline methods. Table 3 summarizes the results for probabilistic forecasting. We observe that in most cases, the proposed flow based algorithms outperform the competitors. MQRNN also shows impressive performance in predicting the forecast quantiles, as it is explicitly trained to minimise the quantile losses. In particular, comparison of GRU+flow with the DeepAR model reveals that even without a sophisticated RNN architecture, the particle flow based approach shows better characterization of prediction uncertainty in most cases. Figure 4 provides a qualitative comparison of the uncertainty characterization, showing example confidence intervals for 15-minute ahead prediction for the PeMSD7 dataset. We see that the proposed algorithm provides considerably tighter intervals, while still achieving coverage of the observed values.
Generalization of particle flow inference across architectures : Table 2 shows that in comparison to deterministic encoder-decoder based sequence to sequence prediction models, the proposed flow based approaches perform better in almost all cases for three different architectures of the RNN. In each case, both of the encoder-decoder model and our approach use a 2-layer architecture with 64 RNN units.
Comparison to the particle filter : Table 4 demonstrates the effectiveness of particle flow (Daum & Huang, 2007), comparing it to a Bootstrap Particle Filter (BPF) (Gordon et al., 1993) with the same number of particles. The use of the bootstrap particle filter leads to a computationally faster algorithm (requiring approximately 60% of the training time of the particle flow-based method).
Comparison to ensembles : We compare the proposed approach with an ensemble of competitive deterministic forecasting techniques. We choose the size of the ensemble so that the algorithms have an approximately equal execution time. We use AGCRN and GMAN to form the ensembles, as they are the best point-forecast baseline algorithms. From Table 5, we observe that the proposed AGCGRU+flow achieves lower average CRPS compared to the ensembles in all cases.
Point forecasting results on non-graph datasets : We evaluate our proposed flow-based RNN on the Electricity and Traffic datasets, following the setting described in Appendix C.4 in (Oreshkin et al., 2020). We augment the results table in (Oreshkin et al., 2020) with the results from an FC-GRU (a fully connected GRU encoder-decoder) and GRU+flow. We use a 2 layer GRU with 64 RNN units in both cases. We follow the preprocessing steps in (Oreshkin et al., 2020). In the literature, four different data splits have been used for the Electricity dataset, and three different splits have been used for the Traffic dataset. The evaluation metric is P50QL (Normalized Deviation).
In Table 6, we observe that the flow based approach performs comparably or better than the state-of-the-art N-BEATS algorithm for the Electricity dataset, even with a simple GRU as the state transition function. The better performance of the univariate N-BEATS compared to the multivariate models suggests that most time-series in these datasets do not provide valuable additional information for predicting other datasets. This is in contrast to the graphbased datasets, where the performance of N-BEATS was considerably worse than the multivariate algorithms. The proposed flow-based algorithm achieves prediction performance on the Traffic dataset that is comparable to N-BEATS except for one split with limited training data. Across all datasets and split settings, our flow-based approach signifi-cantly outperforms the FC-GRU. The proposed algorithm outperforms TRMF, DeepAR, DeepState and DeepGLO. It outperforms DeepFactors for the Electricity dataset, but is worse for the Traffic dataset (for the same split with limited available training data).

Probabilistic forecasting results on non-graph datasets :
For comparison with state-of-the-art deep learning based probabilistic forecasting methods on standard non-graph time-series datasets, we evaluate the proposed GRU+flow algorithm following the setting in (Rasul et al., 2021). The results reported in Table 1 of (Rasul et al., 2021) are augmented with the results of the GRU+flow algorithm. We use a 2 layer GRU with 64 RNN units in each case. We follow the preprocessing steps as in (Salinas et al., 2019;Rasul et al., 2021). The evaluation metric is (normalized) CRPS sum (defined in the supplementary material), which is obtained by first summing across the different time-series, .083 n/a n/a n/a 0.167 n/a n/a DeepFactors n/a 0.112 n/a n/a n/a 0.225 n/a DeepGLO n/a n/a 0.082 n/a n/a n/a 0. both for the ground-truth test data, and samples of forecasts, and then computing the (normalized) CRPS on the summed data. The results are summarized in Table 7. We observe that the proposed GRU+flow achieves the lowest CRPS sum for all datasets.
Computational complexity : For simplicity, we consider a GRU instead of a graph convolution based RNN and we only focus on one sequence instead of a batch. Our model has to perform both GRU computation and particle flow for the first P time steps and then apply the GRU and the linear projection for the next Q steps to generate the predictions. For an L-layer GRU with d x RNN units and N -dimensional input, the complexity of the GRU operation for N p particles is O((P + Q)N p LN d 2 x ) (Chung et al., 2014). The total complexity of the EDH particle flow (Choi et al., 2011) is O(P N λ N 3 ) for computing the flow parameters and O(P N p N λ N d 2 x ) for applying the particle flow (more details in the supplementary material). The total complexity of the measurement model for N p particles is O (QN p N d 2 x ). Since in most cases N d x and N N p , the complexity of our algorithm for forecasting of one sequence is O (P N λ N 3 ). Many of the other algorithms exhibit a similar O(N 3 ) complexity, e.g. TRMF, GMAN. We specify the execution time and memory usage in the supplementary material. Scaling the proposed methodology to extremely high dimensional settings is of significant importance and can be addressed in several ways.
For spatio-temporal predictions using the graph-based recurrent architectures, this can be done if the graph can be partitioned meaningfully. For non-graph datasets, we can use the cross-correlation among different time-series to group them into several lower-dimensional problems. Alternatively, we can train a univariate model based on all the time-series as in (Rangapuram et al., 2018).

Conclusion
In this paper, we propose a state-space probabilistic modeling framework for multivariate time-series prediction that can process information provided in the form of a graph that specifies (probable) predictive or causal relationships. We develop a probabilistic forecasting algorithm based on the Bayesian inference of hidden states via particle flow. For spatio-temporal forecasting, we use GNN based architectures to instantiate the framework. Our method demonstrates comparable or better performance in point forecasting and considerably better performance in uncertainty characterization compared to existing techniques.

Particle flow for Bayesian inference in a state-space model 8.1. Background
Particle flow is an alternative to particle filters for Bayesian filtering in a state-space model. Recall the first order Markov model specified in eqs. (4), (5), and (6) in Section 5.1 of the main paper.
Here y t is the observation from the state-space model at time t. x t and z t denote the unobserved state variable and observed covariates at time t respectively. The filtering task is to compute the posterior distribution of the state trajectory p Θ (x t |y 1:t , z 1:t ) recursively. Suppose we have a set of N p samples (particles) {x t−1 } Np j=1 which approximates the posterior distribution of x t−1 .
In the 'predict' step, we approximate the predictive posterior distribution at time t as follows: p Θ (x t |y 1:t−1 , z 1:t ) = p ψ,σ (x t |x t−1 , y t−1 , z t ) p Θ (x t−1 |y 1:t−1 , z 1:t−1 )dx t−1 , where, the particles {x j t } Np j=1 from the predictive posterior distribution p Θ (x t |y 1:t−1 , z 1:t ) are obtained by propagating {x j t−1 } Np j=1 through the state-transition model specified by eq. (16). Subsequently, the 'update' step applies Bayes' theorem to compute the posterior distribution at time t as follows: For non-linear state space models, particle filters (Gordon et al., 1993;Doucet & Johansen, 2009) employ importance sampling to approximate the 'update' step in eq. (20). However, constructing well-matched proposal distributions to the posterior distribution in high-dimensional state-spaces is extremely challenging. A mismatch between the proposal and the posterior leads to weight degeneracy after resampling, which results in poor performance of particle filters in high-dimensional problems (Bengtsson et al., 2008;Snyder et al., 2008;Beskos et al., 2014). Instead of sampling, particle flow filters offer a significantly better solution in complex problems by transporting particles continuously from the prior to the posterior (Daum & Huang, 2007;Ding & Coates, 2012;Daum & Huang, 2014;Daum et al., 2017).

Particle flow
In a given time step t, particle flow algorithms (Daum & Huang, 2007;Daum et al., 2010) solve differential equations to gradually migrate particles from the predictive distribution such that they represent the posterior distribution for the same time step after the flow. A particle flow can be modelled by a background stochastic process η λ in a pseudotime interval λ ∈ [0, 1], such that the distribution of η 0 is the predictive distribution p Θ (x t |y 1:t−1 , z 1:t ) and the distribution of η 1 is the posterior distribution p Θ (x t |y 1:t , z 1:t ).
Since particle flow only considers migration of particles within a single time step, we omit the time index t in η λ , y, and z to simplify notation.
In (Daum et al., 2010), an ordinary differential equation (ODE) with zero diffusion governs the flow of η λ : If the predictive distribution and the additive measurement noise is Gaussian and the measurement function h is linear, the Exact Daum-Huang (EDH) flow is given as: where, b(λ) = (I + 2λA(λ))[(I + λA(λ))P H T R −1 y + A(λ)η 0 ] , Hereη 0 andP are the mean vector and the covariance matrix of the predictive distribution respectively. For a general nonlinear state-space model, we usually a Gaussian approximation of the predictive distribution based on sample estimates ofη 0 andP or via the Extended Kalman Filter (EKF).

11:
Apply particle flow to all particles: η j m = η j m−1 + m A(λ m−1 )η j m−1 + b(λ m−1 ) 12: end for 13: Set x j t = η j 1 for 1 j N p Numerical integration is usually used to solve the ODE in Equation (22). The integral between λ m−1 and λ m for 1 m N λ , where λ 0 = 0 and λ N λ = 1, is approximated via the Euler update rule. For the j-th particle, the EDH flow in the m-th pseudo-time interval becomes: where the step size m = λ m − λ m−1 and N λ m=1 m = 1.
We start the particle flow from η j 0 =x j t and after the flow is complete, we set x j t = η j 1 to approximate the posterior distribution of x t as: The overall EDH particle flow algorithm is summarized in Algorithm 2. Figure 5 demonstrates the migration of the particles from the prior to the posterior distribution for a Gaussian predictive distribution and a linear-Gaussian measurement model.

Model training
Algorithm 3 summarizes the learning of the model parameters Θ, described in Section 5.2.2 of the main paper.  Update the system parameters using SGD algorithm: The statistics of the PeMS datasets and the non-graph datasets used in our experiments are summarized in Tables 8 and 9 respectively. The description of the PeMS datasets are provided in Section 6.1 of the main paper. The Electricity dataset contains electricity consumption for 370 clients. The Traffic dataset is composed of 963 time-series of lane occupancy rates. The Taxi dataset contains counts of taxis on different roads and the Wikipedia dataset specifies clicks to web links.

Definitions of evaluation metrics
The point forecasts are evaluated by computing mean absolute error (MAE), mean absolute percentage error (MAPE),  ∈ R N denote the ground truth and the prediction at horizon t for m-th test example respectively. The average MAE, MAPE, and RMSE at horizon t are defined as follows: For comparison among the probabilistic forecasting models, we compute the Continuous Ranked Probability Score (CRPS) (Gneiting & Raftery, 2007), and the P10 and P90 Quantile Losses (QL) (Salinas et al., 2020;Wang et al., 2019). Let F (·) be the Cumulative Distribution Function (CDF) of the forecast of the true value x ∈ R. We denote by 1{x z} the indicator function that attains the value 1 if x z and the value 0 otherwise. The continuous ranked probability score (CRPS) is defined as: CRPS is a proper scoring function, i.e., it attains its minimum value of zero when the forecast CDF F is a step function at the ground truth x. The average CRPS at horizon t is defined as the average marginal CRPS across different time-series.
where F (m) t,i (·) is the marginal CDF of the forecast at horizon t for i-th time-series in m-th test example.
Let F (m) t,sum (·) is the CDF of the sum of the forecasts of all time-series at horizon t in m-th test example. The (normalized) CRPS sum is defined as: For a given quantile α ∈ (0, 1), a true value x, and an αquantile predictionx(α) = F −1 (α), the α-quantile loss is defined as: The average (normalized) quantile loss (QL) is defined as follows: The P10QL metric is obtained by setting α = 0.1 in eq. (34); the P90QL metric corresponds to α = 0.9 and the ND (P50QL) metric is obtained using α = 0.5.

Baseline algorithms
For the experiments on the PeMS road traffic datasets, we compare the proposed AGCGRU+flow algorithm with four different classes of forecasting techniques, listed as follows: Graph agnostic statistical and machine learning based point forecasting models: • HA (Historical Average): uses the seasonality of the historical data.
Spatio-temporal point forecast models: • DCRNN (Li et al., 2018): Diffusion Convolutional Recurrent Neural Network, combines diffusion convolution with GRU to form an encoder-decoder architecture for sequence to sequence prediction.
• ASTGCN (Guo et al., 2019): Attention Spatial-Temporal Graph Convolutional Network, spatial and temporal attentions to learn recent and seasonal patterns.
• GMAN (Zheng et al., 2020): Graph Multi-Attention Network, multiple spatio-temporal attention blocks to form an encoder-decoder architecture, transform attention between encoder and decoder.
• LSGCN (Huang et al., 2020): Long Shortterm Graph Convolutional Network, a novel attention mechanism and graph convolution, integrated into a spatial gated block.
• N-BEATS (Oreshkin et al., 2020): Neural Basis Expansion Analysis for Interpretable Time-Series, an univariate model, built using backward and forward residual connections and deep stack of fully-connected layers.
Deep-learning based probabilistic forecasting methods: • DeepAR (Salinas et al., 2020): RNN based probabilistic method using parametric likelihood for forecasts.   Tables 10, 11, 12, and 13. We observe that statistical models such as HA, ARIMA, and VAR and basic machine learning models such as SVR, FNN, and FC-LSTM show poor predictive performance as they cannot model the complex spatio-temporal patterns present in the real world traffic data well. Graph agnostic deep learning models such as DeepGLO and N-BEATS perform better than the statistical models, but they cannot incorporate the graph structure when learning. FC-GAGA has lower forecasting errors as it is equipped with a graph learning module. The spatio-temporal graph-based models (especially AGCRN, GMAN, GWN, and LSGCN) display better performance. These models either use the observed graph or learn the graph structure from the data. In general, the deep learning based probabilistic forecasting algorithms such as DeepAR, DeepFactors, and MQRNN do not account for the spatial relationships in the data as well as the graph-based models, although MQRNN is among the best performing algorithms. DeepAR and DeepFactors aim to model the forecasting distributions and thus do not perform as well in the point forecasting task. The training loss function (negative log likelihood of the forecasts) does not match the evaluation metric. However, MQRNN shows better performance, possibly because it does target learning the median of the forecasting distribution along with other quantiles. The proposed AGCGRU+flow algorithm demonstrates comparable prediction accuracy to the best-performing spatio-temporal models and achieves the best average ranking across the four datasets. Figure 6 demonstrates that the proposed AGCGRU+flow has lower average MAE in most of the nodes compared to the second best performing AGCRN algorithm, for all four PeMS datasets. Some qualitative visualization of the confidence intervals for 15-minute ahead predictions for the PeMSD3, PeMSD4, PeMSD7, and PeMSD8 datasets are shown in Figures 7, 8, 9, and 10 respectively. We observe that the confidence intervals from the proposed algorithm are considerably tighter compared to its competitors in most cases, whereas the coverage of the ground truth is still ensured.

Detailed results for comparison with particle filter
In Table 4 of the main paper, we compare the average MAE and average CRPS of the proposed AGCGRU+flow with a Bootstrap Particle Filter (BPF) (Gordon et al., 1993) based approach. Tables 14 and 15 provide the detailed comparison both in terms of point forecasting and probabilistic forecasting metrics. We observe that the proposed AGCGRU+flow algorithm outperforms the particle filter based approach in most cases.

Effect of number of particles
For this experiment, we consider three different settings with varying number of particles N p = 1/10/50 for testing. The model is trained using 1 particle in each case. From Table 16, we observe that increasing the number of particles cannot improve the point forecasting accuracy significantly, whereas the results in Table 17 show that characterization of the prediction uncertainty is improved as more particles are used to form the approximate posterior distribution of the forecasts.

Effect of different learnable noise variance at each node
In this experiment, we compare the proposed state-space model with different learnable noise variance at each node (parameterized by the softplus function in eq. (9) in the main paper with fixed and uniform noise standard deviation γ = 0.01/0.05/0.10 at all nodes. Other hyper-parameters and the training setup remain unchanged. The results in Table 18 demonstrate that the learnable noise variance approach is not particularly beneficial in comparison to a uniform, fixed variance approach in most cases. However, we note that the probabilistic metrics reported in Table 19 are the lowest for the learnable noise variance model in all cases. This suggests that different time-series in these road traffic datasets have different degrees of uncertainty which cannot be effectively modelled by the uniform, fixed noise variance approach.

Detailed comparison with deterministic encoder-decoder models
In Table 2 of the main paper, we compare the average MAE of the proposed flow based approaches with those of the deterministic encoder-decoder based sequence to sequence prediction models for three different RNN architectures. In Table 20, we report the MAPE and RMSE, in addition to the MAE. We see that the particle flow based RNN models outperform the corresponding deterministic encoder-decoder models in most cases.   In Table 5 of the main paper, we compare the average CRPS of the proposed AGCGRU+flow algorithm with ensembles of AGCRN and GAMN. From Table 21, we observe that our approach is comparable or slightly worse compared to the ensembles in terms of the MAE, MAPE and RMSE of the point forecasts. However, the proposed AGCGRU+flow shows better characterization of the prediction uncertainty compared to the ensemble methods in almost all cases, as shown in Table 22.

Comparison with a Variational Inference (VI) based approach
Although there is no directly applicable baseline forecasting method in the literature that incorporates VI, RNNs, and GNNs, we can derive a variational approach using equivalent GNN-RNN architectures and compare it to the particle flow approach. We wish to approximate p Θ (y P +1:P +Q |y 1:P , z 1:P +Q ). So, the ELBO is defined as follows: L(Θ, Ω) = E qΩ log p Θ (y P +1:P +Q , x 1:P |y 1:P , z 1:P +Q ) − log q Ω (x 1:P |y 1:P +Q , z 1:P +Q ) .
Experimental details : We treat ρ, σ and σ as hyperparameters and set ρ = 1 and σ = σ = 0. This implies that q ψ ,σ is a Dirac-delta function and the maximization of ELBO (in eq. (35)) using SGD (SGVI) amounts to mimization of the same cost function as defined in eq. (14) in the main paper. The only difference is that now a) we have two separate AGCGRUs for encoder and decoder and b) there is no particle flow in the forward pass. We call this model AGCGRU+VI and compare it to AGCGRU+flow. The other hyperparameters are set to the same values as for the AGC-GRU+flow algorithm. From Table 23, we observe that for comparable RNN architectures, the flow based algorithm significantly outperforms the variational inference based approach in the point forecasting task. The results in Table 24 indicate that in the probabilistic forecasting task, both particle flow and VI approaches show comparable performance despite AGCGRU+flow having approximately half of the learnable parameters of the AGCGRU+VI model. Table 25 summarizes the run time, GPU usage during training, and the size of the learned model for AGCRN-ensemble, GMAN-ensemble, and the proposed AGCGRU+flow for the four PeMS datasets. We observe that if we choose the ensemble size so that the algorithms have an approximately equal execution time, then the model-size of the ensemble algorithms are comparable to our approach as well. However, our method requires more GPU memory compared to the ensembles during training because of the particle flow in the forward pass.