{'fraction_non_alphanumeric': 0.035506804793824905, 'fraction_numerical': 0.013243956936827138, 'mean_word_length': 4.579329102447869, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Variational Bayesian Inference for Network Autoregression Models

We develop a variational Bayesian (VB) approach for estimating large-scale dynamic network models in the network autoregression framework. The VB approach allows for the automatic identification of the dynamic structure of such a model and obtains a direct approximation of the posterior density. Compared to Markov Chain Monte Carlo (MCMC) based sampling approaches, the VB approach achieves enhanced computational efficiency without sacrificing estimation accuracy. In the simulation study conducted here, the proposed VB approach detects various types of proper active structures for dynamic network models. Compared to the alternative approach, the proposed method achieves similar or better accuracy, and its computational time is halved. In a real data analysis scenario of day-ahead natural gas flow prediction in the German gas transmission network with 51 nodes between October 2013 and September 2015, the VB approach delivers promising forecasting accuracy along with clearly detected structures in terms of dynamic dependence.


Introduction
Networks have emerged and become available in various fields, such as energy transmission, logistics and transportation, and financial systems. Networks are dynamic in terms of their temporal dependence, and they often have large scales. The understanding and inference of network dynamics have profound implications for operations and decision making in modern industries. Tremendous growth and heterogeneity in both nodes/edges and dependence over time are the key characteristics of such networks. However, conventional statistical methods either assume that networks are static or consider only low-dimensional temporal data. This creates a need for efficient computational approaches that are able to reveal the essential dependence structures in high dimensions and simultaneously deliver accurate inferences with a low computational cost.
Industrial networks contain series of temporal-spatial data collected over time. While the nodes/edges are often fixed or possess trivial changes, the lead-lag temporal dependence issue can no longer be ignored in network inference. Graph theory has been widely used for unraveling structural information in large-scale network analysis. For example, Fan et al. (2009) and Guo et al. (2011) proposed a sparse graphic network. Liu et al. (2012) proposed the semiparametric Gaussian copula graphical model. Despite their efficiency, these graphical models assume static networks, and the evolution of the network dependence is not considered in their estimation processes.
The temporal dependence of a network can be represented in the vector autoregressive (VAR) modeling framework. In the VAR framework, each node is considered as one time series, and the network dependence is measured as the lead-lag cross-correlations among multiple time series. Both the theoretical properties and empirical performance of VAR have been well studied with respect to multivariate data; see Lütkepohl (2007) and Bańbura et al. (2010). However, the application of VAR for large-scale network analysis is still challenging.
Given a network with models, which correspond to time series, and supposing that the dynamics depend on the last lags, to the result is that there are 2 unknown coefficients in the VAR model. When the number of nodes becomes large or the temporal dependence increases, VAR is overparameterized and this leads to low estimation accuracy or even infeasibility with regard to model inference.
A unique feature in industrial networks is that their lead-lag temporal dependences, such as concurrent dependence among their nodes, are less dense than the networks themselves.
Individual networks are also much sparser than other social networks. It is conceivable that large-scale industrial networks, are driven by a few essential and cohesive connections among their nodes to facilitate network evolution. This motivates the modeling of large-scale dynamic networks using sparse VAR. In particular, penalties are imposed on the parameter space of the VAR framework with various possible types of structural assumptions. For the purpose of both estimation and interpretability, structural sparsity can be enabled in elements, groups and lags. Lag sparsity investigates the effect of time-lagged information, while group sparsity highlights the impacts of certain nodes on others. In addition to the universal effect on a group of lag coefficients, a sparse element illustrates a single effect. Basu & Michailidis (2015) investigated the theoretical properties of ℓ 1 -regularized estimates, where multiple time series were assumed to be stable Gaussian processes. Melnyk & Banerjee (2016) established bounds on the non-asymptotic estimation error of the Lasso-type estimator for structured VAR parameters. Nicholson et al. (2014) proposed several structures for VAR and Lasso, group Lasso and sparse group penatly functions to achieve sparsity in the elements and groups of a network; see also Hsu et al. (2008), Song & Bickel (2011), and Chen et al. (2020). Moreover, VAR models can be easily extended to the above three kinds of sparsity by building up a hierarchical lag structure in the autoregression model via the inclusion of high-dimensional exogenous variables.
The estimation of the structured of a VAR framework faces two challenges. First, the sparse structure needs to be specified to avoid the overfitting of the high-dimensional models.
Second, the framework should be able to adapt to different kinds of stochastic behaviors, as empirical data are likely non-Gaussian. Bayesian methods are natural choices because they deliver stable performances without prespecified assumptions about the structure and distribution of the model. Stochastic search variable selection (SSVS), for example, is the most commonly used Bayesian variable selection approach. It introduces latent indicators embedded in the priors and stochastically searches subsets by generating posterior samples with the Markov Chain Monte Carlo (MCMC) algorithm; see George & McCulloch (1993). Geweke (1996) proposed the component-wise Gibbs method for the purpose of improving computational efficiency. By using the spike-and-slab prior, it avoids computing inverse matrices and thus reduces the computational cost. However, the Gibbs sampler requests either random or systematic updating of the coefficients. Chen et al. (2011) proposed the stochastic matching pursuit (SMP) algorithm, which updates the coefficients of each step for obtaining the best fit based on the current residual vector. In terms of structural selection, Farcomeni (2010) introduced Bayesian-constrained variable selection. Chen et al. (2016) proposed the groupwise Gibbs sampler. As a proof of concept, Chu et al. (2019)  To the best of our knowledge, this is the first attempt to derive variational inference for a large-scale dynamic network analysis in a structured NAR/VAR framework.
This paper is organized as follows. Section 2 introduces the dynamic network model in the VAR framework. Several types of structural assumptions are also demonstrated.
Section 3 presents the proposed variational Bayesian algorithms for large-scale dynamic network inference. Section 4 investigates the finite-sample performance of the proposed VB approach. Section 5 reports the network inference for day-ahead gas flow forecasting with 51 high-pressure nodes in the natural gas transmission network in Germany. Section 6 provides a brief conclusion of our work.

Model
Let Y ∈ R 1× denote a vectorized time series of networks with nodes at time over the time period [1, ]. Without loss of generality, Y is demeaned. We consider a dynamic network model in a vector autoregression framework Y which is assumed to depend on the past values of the network at lags 1 to , i.e., Y = Y −1 1 + · · · + Y − + , = + 1, + 2, . . . , .
Here, we let ℓ be an × coefficient matrix for lag-ℓ, ℓ = 1, 2, · · · , , which is used to measure the lead-lag temporal dependence in the network. The number of coefficients in B ℓ grows quadratically with the number of nodes . Given = 51 in the gas transformation network, there are 51 2 = 2601 unknown coefficients for each ℓ . Obviously, the diagonal elements of ℓ indicate the serial dependence of each node on its own lag-ℓ value and the dependences of off-diagonal elements on other nodes' lag-ℓ values. The term { } = +1 is a sequence of serially uncorrelated 1 × random vectors, with a mean vector of zero and a covariance matrix Σ. Then, the dynamic network model can be represented in matrix form as follows: In this paper, we consider a structured NAR/VAR framework, i.e., dynamic dependence is sparsely motivated by a large-scale industrial network, such as the German gas transformation network. Figure  patterns, and we can treat them as a group. Therefore, if one node affects another node, the others are all influenced by this node. As shown in Figure 2, node M4 affects other nodes, and nodes M3, M5 and M6 do not affect other nodes. As such, the coefficient matrix can be separated into a diagonal matrix reflecting the dynamic dependence of each node on its own, and a sparse matrix showing the temporal dependence of each node on others. In network analysis, the columns that have the same row sparsity are grouped together.
• SG Structure: For the segmentation grouping (SG) structure in ℓ , the nodes in the same segment interact with each other but are independent from the other nodes. This is associated with the empirical observation that there are different types of nodes in the network. Figure 3 illustrates the SG structure among municipal and industrial nodes. The network is divided into disjoint segments of the municipal and industrial types. The estimation process can be conducted segment by segment. • NG Structure: Consider the "no grouping" (NG) structure in ℓ . Each node has its own impact on the other nodes over time. This is a scenario where no regular pattern is identified among the nodes. Figure 4 shows an example with different types of nodes.
In fact, there are no similar patterns among them. In this case, we separate them during the inference procedure.
It is easy to see that the first and third types of structures, UG and NG, are two special cases of the SG structure when the numbers of segments are 1 and , respectively. While the three types of structures are defined for each coefficient matrix, they may appear in different time lags. The identification of a proper structure can help not only increase the estimation accuracy but also enhance computational time of the algorithm.

Variational Approximation Algorithm
We introduce the variational Bayesian approach for NAR/VAR structure inference in highdimensional scenarios, and it is expected to automatically choose the proper structure and perform estimation at a low computational cost. Here, we derive only the VB method for the segmentation grouping structure because it nests the other two types of structures as special cases of segmentation grouping. Thus, it is a derivation under a general setup. To simplify the procedure, we assume that all coefficient matrices share the same segmentation grouping structure. Denote = 1 , 2 , · · · , as an overall index set for the columns in each ℓ , where is the index set of the th segment with size 0 < | | < and =1 | | = , and is the total number of segments (groups). In this section, we start by introducing the Bayesian structure selection approach, and then we mention the details of the proposed VB method.

The Bayesian structure selection algorithm
For the Bayesian hierarchical model used to select segmentation structures, Chu et al. (2019) introduced two indicators ℓ, and ℓ, , to identify the active structures. These indicators denote the active nonzero coefficients of the th segment in the th row in ℓ with respect to the lag values of itself and others. For its own lags, ℓ, = 1 denotes that the element of the th row and the th column in the coefficient matrix for lag-ℓ, ℓ, , , is nonzero, and ℓ, = 0 indicates a coefficient of zero, i.e., ℓ, , = 0. In addition, ℓ, , = 1 indicates that the th row in ℓ and the columns of the th segment for the off-diagonal elements in ℓ are all nonzero, i.e. , , ≠ 0, and , , = 0 otherwise. Consider the Bayesian structure selection approach for determining the segmentation structure in the NAR/VAR model.
Based on these prior assumptions regarding the coefficients and the additional inverse-Wishart prior for the covariance matrix of the NAR/VAR model, coefficient inference can be performed by the vector autoregression Gibbs sampler (VAGSA, Chu et al., 2019). Basically, in the VAGSA, we need to iteratively generate the posterior samples of the two indicators and coefficients for the further inference. Similar to other MCMC algorithms, the VAGSA is computationally expensive, especially when the number of nodes increases. Instead of generating the posterior samples as an approximation of the posterior distribution, the variational Bayesian approach is adopted in our study to directly obtain the approximation of the posterior density function. Here, an EM-type algorithm is used to solve the corresponding optimization problem, and this is expected to significantly improve the computational efficiency of our approach.

The Variational inference procedure
Before introducing the variational Bayesian method, we first reparametrize the NAR/VAR model. Recall that in the VAGSA, the coefficient prior is a mixture distribution with a normal distribution and a delta function at point zero. When an indicator is equal to zero, we can simply set the corresponding coefficient to zero. To simplify the model structure, we reparametrize the coefficient as the product of the coefficient and the indicator. That is, ℓ, , = ℓ, ℓ, , and ℓ, , = ℓ, , ℓ, , . Here, the priors of the indicators, , and , , , are still independent Bernoulli distributions with probabilities 1 and 2 , respectively. In addition, the prior of ℓ, , is chosen as a normal distribution with mean zero and variance 2 , and the prior of ℓ, , comes from the multivariate normal distribution with a mean zero vector and a covariance matrix 2 I | | . Among them, the priors of ℓ, , and ℓ, , are independent of ℓ, and ℓ, , . Thus, ℓ, ℓ, , and ℓ, , ℓ, , have the same effects as those shown in Eqs. (2) and (3).
Instead of generating the posterior samples directly, the variational Bayesian approach identifies the best approximate distribution of the true posterior for the further Bayesian inference. According to the ordinary variational Bayesian approach (Bishop, 2006), the Kullback-Leibler divergence (KL divergence) is used to measure the dissimilarity of the true posterior distribution from the approximated posterior distribution. Let = 1 , 2 , Σ, 2 be the set of parameters and η, γ, B be the set of the indicator variables and coefficient matrices, where η is the vector of ℓ, , , γ is the vector of ℓ, , and B = 1 , 2 , . . . , .
Given the prior assumptions, the posterior density function of η, γ, B is proportional to the the following joint density function: Define (η, γ, B) as an approximate posterior density function of ( B, η, γ |Y , X). Since we have According to Eq. (5) where ℓ, ( ℓ, , , ℓ, ) and ℓ, , ( ℓ, , , ℓ, , ) are chosen from the corresponding prior distributions. Then, a variational extension of the EM algorithm is adopted in the estimation process by maximizing the corresponding ( ) with respect to the parameters. In the Estep, one would take the expectation of ( ) with respect to η, γ and , and then in the M-step, one optimizes ( ) with respect to . Iterate the two steps until the lower bound ( ) converges.
In the E-step, (η, γ, ) is updated as follows: Here, X ( ) ℓ and X (− ) ℓ denote the th column of X ℓ and matrix X ℓ excluding the th column, respectively. The same definition structure is used for Y and . In addition, (− , ) In the M-step, we take the derivative of ( ) with respect to and then set it as a zero vector. Thus, the components in can be updated as the solutions of the normal equations, i.e., More details regarding both steps are shown in the supplementary materials.
Since we iterate the E-and M-steps sequentially, a natural choice of a stopping criterion is that the difference between the values of ( ) in two consecutive iterations is less than a certain threshold. Then one can approximate the posterior inclusion probabilities as follows: ( ℓ, , = 1|Y , X,ˆ) ≈ ( ℓ, , = 1|ˆ) = 2,ℓ, , , Thus, based on the median probability criterion (Barbieri & Berger, 2004), ℓ, , or ℓ, , is identified as non zero if 1, , or 2, , , is larger than or equal to 1/2. Finally, the one-step head predictionŶ +1 can be obtained byŶ +1 = Yˆ1 + Y −1ˆ2 + . . . , Y +1−ˆ, where thê ℓ are estimated based on the identified active structures. The EM-type method is a locally optimal approach and may be sensitive to the initial status. Realistic initial values of may provide poor estimation results due to improper choices of the initial prior probabilities. In our study, the initial values of the prior probabilities, 1 and 2 , are set to be small, e.g., 1 = 2 = 0.01. The initial values of the coefficient matrix, , are obtained via the least-squares estimation method, and the initial value of Σ is set as the sample variance of Y divided by 2. The threshold value of the stopping criterion is usually set to a prespecified value, which is less than or equal to 10 −6 .

Simulation
This section investigates the finite-sample performance of the proposed VB approach compared with that of a known data generation process. We first revisit the simulation studies of VAGSA in Chu et al. (2019). Given dynamic networks, we care not only about the structure selection ability but also about the computational efficiency of our approach. For mediumsized examples, i.e., = 10 and 20, we directly compare the performances of the VB method and the VAGSA. When the dimensionality increases to = 50, the VAGSA requires an extremely long computational time, and thus, we illustrate only the results of the proposed VB method.

Medium-sized network examples
We follow the simulation setups in Chu et al. (2019) and generate medium-sized network series with = 10 and 20 for a fair comparison. We assume that there are multiple lags in the VAR framework. Three dependence structures are considered in the simulations. In addition, the error terms are generated from a multinormal distribution with zero mean.
For each simulation, we generate data with = 301, where the first 300 samples are used for model training and the last point 301 is used to demonstrate the prediction ability of the proposed method. Each simulation is replicated = 100 times, and the segmentation structure is assumed to be known. We simply fix the number of lags to 10 and set the threshold value of the stopping criteria to 10 −8 in all simulations.
As mentioned before, we also implement the VAGSA for the simulation cases for com-

Large-scale network examples
To illustrate the feasibility of the proposed VB approach for high-dimensional scenarios, we conduct simulations with = 50 nodes. The setups of the large-scale simulations are similar to those used in Section 4.1. Here, we also consider the three different grouping structures, and the details of three cases are shown as follows.
In this simulation example, we generate data with = 701. The first 700 samples are used for model training, and the last sample 701 is used to illustrate the prediction ability of the model. where is the number of replicates and

Simulation results
Last, the average CPU time for a replication is also reported in seconds. Here, we run the R code for the proposed VB method, and we run the VAGSA based on its MATLAB code. The good performance of the VB method continues for the large-scale networks because it is computationally possible. Figures 6 and 7 show the average estimates of the coefficient matrices for two different covariance matrices. These coefficient matrices are all close to the true matrices. Use the UG case as an example. The structure identification results are perfect because the TPR = 100%, the FPR = 0% and the AMS is still close to the true value. In the NG case with a covariance matrix of, Σ 50 , the corresponding TPR is 92%, which is slight lower than the TPRs in the other cases. This may be because the VB method is used to approximate the true posterior distribution, and thus, active variables with small coefficients might not be detected. Overall, the proposed VB method has good performance in identifying the active structures for large-scale network cases. In addition, the average CPU times show that the case with a higher level structure, namely, UG, requires the least computational time, followed by the SG and NG cases.

Empirical Study
We apply the VB approach to estimate the temporal dependences of the German natural gas flow networks and perform day-ahead forecasting with the NAR/VAR framework. The data cover two years, from 1st October 2013 to 30th September 2015. The gas flows contain both inflows (supply) and outflows (demand) recorded 7 days a week at 51 distribution nodes belonging to 4 categories with different functions. There are 34 municipal nodes, that provide gas to local residential areas and small business districts. The 11 industry nodes are responsible for factory production. Moreover there is 1 border node, which is important in Germany, because these nodes serve as network transfer points for natural gas imported  In addition, note that the data in municipal nodes have been seasonally adjusted via the daily temperature, and 3 nodes, O3, O4, and O5, in the "others" category have also been seasonally adjusted.
We adopt the NAR/VAR model with the segmentation structure to analyze the gas network data. Here, the segmented grouping structures are defined based on the types of nodes that are involved. That is, each type of node is treated as a group. In addition, we set the number of lags as 14 to incorporate social dependence up to two weeks ahead.
We use the data from Oct 1st, 2013, to March 31, 2015, as the in-sample training data.
The threshold value of the stopping criterion for the proposed VB method is set to be 10 −6 , and the other initial settings are the same as those used we set in the simulation studies in Section 4. After learning the model via the VB method, we perform a one-step ahead forecast. That is, at each daily point, we shift one more day, use the expanded sample to retrain the model coefficients based on the active segmentation structures, and perform a one-day ahead forecast. Note that we repeat the iterative forecast for the remaining 0.5 years.
Overall, we obtain 183 forecasts. To illustrate the performance of the proposed VB method, in addition to the mean absolute percentage estimator (MAPE), the normalized MSE (NRMSE) for the 51 nodes is also used, and both measures are defined as follows: Figures 9 displays the estimated coefficient matrices for lag-1 to lag-3, and these three coefficient matrices are quite sparse. For lag-1, it can be observed that the 4th and 14th rows, corresponding to municipal nodes M4 and M14, affect the other municipal nodes. In the 35th row, an industry node, I1, influences the other industry nodes but not the other three types. The 51st row is the status of node O5; the others node only have relations with municipal nodes. On the other hand, the only border node is only related to itself and has no interactions with the others. For lag-2, only the 35th row I1 influences the nodes in the other category. For lag-3, some nodes have minor influences by themselves. As the lag increases, there are fewer and fewer nodes that may influence others. Furthermore, there is no more serial dependence in the network after lag-10. Basically, the VB method delivers good performances for both the in-sample training period from 1st October 2014 to 31st March 2015 and the out-of-sample forecasting period from Figure 9: The coefficient matrix for lag-1 to lag-3. There are several possible future research directions. The first concerns the parameter turning strategy for the proposed VB method. In our experience, the performance may be sensitive to the chosen initial parameters, especially for the prior probabilities 1 and 2 .
According to Ormerod et al. (2017)  A Appendix  In this extra supplementary document, we provide more details about the EM algorithm for the VB approach with respect to the NAR/VAR model.