Generalized Network Autoregressive Processes and the GNAR Package

This article introduces the GNAR package, which ﬁts, predicts, and simulates from a powerful new class of generalized network autoregressive processes. Such processes consist of a multivariate time series along with a real, or inferred, network that provides information about inter-variable relationships. The GNAR model relates values of a time series for a given variable and time to earlier values of the same variable and of neighboring variables, with inclusion controlled by the network structure. The GNAR package is designed to ﬁt this new model, while working with standard ‘ ts ’ objects and the igraph package for ease of use.


Introduction
Increasingly within the sciences, networks and network methodologies are being used to answer research questions. Such networks might be observed, such as connections in communication network or information flows within, or they could be unobserved: inferred networks that can explain a process or effect. Given the increase in the size of data sets, it may also be useful to infer a network from data to efficiently summarize the data generating process.
We consider time series observations recorded at different nodes of a network, or graph. Our GNAR package (Leeming, Nason, and Nunes 2020) and its novel generalized network autoregressive (GNAR) statistical models focus on partnering a network with a multivariate time series and modeling them jointly. One can find an association network, see, e.g., Chapter 7 of Kolaczyk (2009), or Granger causality network, e.g., Dahlhaus and Eichler (2003), between different variables by analysing a multivariate time series and its properties. However, here we assume the existence of an underlying network and use it during the analysis of the time series, although sometimes its complete structure is unknown.

Network time series processes
We assume that our multivariate time series follows an autoregressive-like model at each node, depending both on the previous values of the process at that node, and on neighboring nodes at previous time steps. These neighboring nodes are included as part of the network structure, as defined below.

Network terminology and notation
Throughout we assume the presence of one or more networks, or graphs, associated with the observed time series. Each univariate time series that makes up the multivariate time series occurs, or is observed at, a node, or location on the graph(s). These nodes are connected by a set of edges, which may be directed, and/or weighted.
We denote a graph by G = (K, E), where K = {1, . . . , N } is the set of nodes, and E is the set of edges. A directed edge from node i ∈ K to j ∈ K is denoted i j, and an undirected edge between the nodes is denoted i j. The edge set of a directed graph is E = {(i, j) : i j; i, j ∈ K}, and similarly for the set of un-directed edges.

Stage-r neighborhoods
We introduce the notion of neighbors and stage-neighbors in the graph structure as follows; for a subset A ⊂ K the neighbor set of A is given by N (A) = {j ∈ K/A : i j; i ∈ A}. These are the first neighbors, or stage-1 neighbors of A. The rth stage neighbors of a node i ∈ K are given by N (r) (i) = N {N (r−1) (i)}/[{∪ r−1 q=1 N (q) (i)} ∪ {i}], for r = 2, 3, . . . and N (1) (i) = N ({i}). Figure 1 shows an example graph, where node E has stage-1 neighbor A, stage-2 neighbor D, and stage-3 neighbors B and C. Neighbor sets for this example include N (1) (D) = {A, B, C}, and N (3) (E) = {B, C}. In the time-varying network setting, a subscript t is added to the neighbor set notation.

Connection weights
Each network can have connection weights ω ∈ [0, 1] associated with every pair of nodes. This connection weight can depend on the size of the neighbor set and also encodes any edge-weight information. Formally, the values of the connection weights from a node i ∈ K to its stage-r neighbor j ∈ N (r) (i) will be the reciprocal of the number of stage-r neighbors; ω i,j = |N (r) (i)| −1 , where |·| denotes the cardinality of a set. In Figure 1 the connection weights would be, for example, ω E,A = 1, ω A,E = ω A,D = 0.5. Connection weights are not necessarily symmetric, even for an un-directed graph. We note that this choice of these inverse distance weights is one of many possibilities, and some other means of creating connection weights could be used.
When the edges are weighted, or have a distance associated with them, we use the concept of distance to find the shortest path between two vertices. Let the distance from node i to be denoted d i, ∈ R + , and if there is an un-normalized weight between these nodes, denote this µ i, ∈ R + . To find the length of connection between a node i and its stager neighbor, k, we sum the distances on the paths with r edges from i to k and take the minimum (note that there are no paths with fewer edges than r as k is a stage-r neighbor). If the network includes weights rather than distances, we find the shortest r length path where d i, = µ −1 i, . Then the connection weights between node i and its stage-r neighbor k are either for a network with weights. This definition means that all nodes will have connection weights that sum to one for any non-empty neighbor set, whether they are in a sparse or dense part of the graph.

Edge or node covariates
A further important innovation permits a covariate that can be used to encode edges effects (or nodes) into certain types. Our covariate will take C ∈ N discrete values and be indexed by c. A more general covariate could be considered, but we wish to keep our notation simple in the definition that follows. For example, in an epidemiological network we might have two edge types: one that carries information about windborne spread of infection and the other carries information about identified direct infections. The covariates do not change our neighbor sets or connection weight definitions, so we have the property for all i ∈ K and r ∈ N such that N (r) (i) is non-empty.

The generalized network autoregressive model
Consider an N × 1 vector of nodal time series, X t = (X 1,t , . . . , X N,t ) , where N is considered fixed. Our aim is to model the dependence structure within and between the nodal series using the network structure provided by (potentially time-varying) connection weights, ω. For each node i ∈ {1, . . . , N } and time t ∈ {1, . . . , T }, our generalized autoregressive model of order (p, [s] where p ∈ N is the maximum time lag, [s] = (s 1 , . . . , s p ) and s j ∈ N 0 is the maximum stage of neighbor dependence for time lag j, with N 0 = N ∪ {0}, N (r) t (i) is the rth stage neighbor set of node i at time t, ω (t) i,q,c ∈ [0, 1] is the connection weight between node i and node q at time t if the path corresponds to covariate c. Here, we consider a sum from one to zero to be zero, i.e., 0 r=1 (·) := 0. The α i,j ∈ R are "standard" autoregressive parameters at lag j for node i. The β j,r,c ∈ R correspond to the effect of the rth stage neighbors, at lag j, according to covariate c = 1, . . . , C. Later, we derive conditions on the model parameters to achieve process stationarity over the network. Here the noise, {u i,t }, is assumed to be independent and identically distributed at each node i, with mean zero and variance σ 2 i . Our model meaningfully enhances that of the arXiv publication Knight, Nunes, and Nason (2016) by now additionally including different autoregressive parameters, connection weights at each node and, particularly, parameters β that depend on covariates. Note that the IID assumption on the noise {u i,t } could of course be relaxed to include correlated innovations.
We note that crucially, the time-dependent network topology is integral to the model parametrization through the use of time-varying weights and neighbors. These features yield a model that is sensitive to the network structures and captures contemporaneous as well as autoregressive relationships, as defined by Equation 1. The network should therefore be viewed not as an estimable quantity, but as a time-dependent known structure.
In the GNAR model, the network may change over time, but the covariates stay fixed. This means that the underlying network can be altered over time, for example, to allow for nodes to drop in and out of the series but model fitting can still be carried out. Practically, this is extremely useful, as shown by the example in Section 4. Our model allows for the α parameters may be different at each node, however the interpretation of the network regression parameters, β j,r,c , is the same throughout the network.
A more restrictive version of the above model is the global-α GNAR(p, [s]) model, which has the same autoregressive covariate at each node, where the α i,j are replaced by α j . This defines a process with the same behavior at every node, with differences being present only due to the graph structure.

GNAR network example
Networks in the GNAR package are stored in a list with two components edges and dist. The edges component is itself a list with N slots each containing a vector whose entries are indices to their neighboring nodes. For example, if 3 4 denotes an undirected edge between nodes 3 and 4 then the vector edges[[3]] will contain a 4 and edges [[4]] will contain a 3. If the network is undirected this will mean that each edge is "double counted" in summary information. A directed edge 3 4 would be listed in edges [[3]] as a 4, but not edges [[4]] if there is no edge in the opposite direction. The dist component is of the same format as edges, and contains the distances corresponding to the edge links, if they exist. For example, in an un-weighted setting, the connection weights are such that all neighbors of a node have equal effect on the node. This is achieved by setting all entries of the dist component to one, and the software calculates the connection weight from these. A GNAR network is stored in a 'GNARnet' object, and an object can be checked using the is.GNARnet function. The S3 methods plot, print, and summary are available for 'GNARnet' objects. Figure 1 shows an example that is stored as a 'GNARnet' object called fiveNet and can be reproduced using R> library("GNAR") R> library("igraph") R> set.seed(1) R> plot(fiveNet, vertex.label = c("A", "B", "C", "D", "E")) The basic structure of the 'GNARnet' object is, as usual, displayed with

Converting a network to GNARnet form
Our GNARnet format integrates with other methods of specifying a network via a set of functions that generate a 'GNARnet' from others, such as an 'igraph' object.
We can also use the GNARtoigraph function to extract graphs involving higher-order neighbor structures, for example, creating a network of third-order neighbors.
In addition to interfacing with igraph, we can convert between 'GNARnet' objects and adjacency matrices using functions as.matrix and matrixtoGNAR. We can produce an adjacency matrix for the 'fiveNet' object with R> as.matrix(fiveNet) and an example converting a weighted adjacency matrix to a 'GNARnet' object is R> set.seed(920) R> adj <-matrix(runif(9), ncol = 3, nrow = 3) R> adj[adj < 0.3] <-0 R> print(matrixtoGNAR(adj)) WARNING: diagonal entries present in original matrix, these will be removed GNARnet with 3 nodes edges:1--3 3--1 3--2 edges of unequal lengths Note that as the output suggests, the matrixtoGNAR function removes self-loops (represented by ones in the diagonal of the associated adjacency matrix). This is because our model in Section 2 accounts for nodal self-dependence by allowing node-specific autoregressive parameters α i,j .

Example: GNAR model fitting
The fiveNet network has a simulated multivariate time series associated with it of class 'ts' called fiveVTS. The pair together are a network time series. The object can be loaded in the usual way using the data function. GNAR contains functions for fitting and predicting from GNAR models: GNARfit and the predict method, respectively. These make use of the familiar R command lm, since the GNAR model can be essentially re-formulated as a linear model, as we shall see in Section 3 and Appendix B. As such, least squares variance/standard error computations are also readily obtained, although other, e.g., HAC-type variance estimators could also be considered for GNAR models.
Suppose we wish to fit the global-α network time series model GNAR(2, [1, 1]), a model with four parameters in total. We can fit this model with the following code.
R> data("fiveNode", package = "GNAR") R> answer <-GNARfit(vts = fiveVTS, net = fiveNet, alphaOrder = 2, + betaOrder = c(1, 1)) R> answer In this fit, the global autoregressive parameters areα 1 ≈ 0.206 andα 2 ≈ 0.021 and the β network parameters areβ 1,1,1 ≈ 0.503 andβ 2,1,1 ≈ −0.095. Also, the network edges only have one type of covariate so C = c = 1. We can just look at one node. For example, the model at node A is The model coefficients can be extracted from a 'GNARfit' object using the coef function as is customary. The 'GNARfit' object returned by GNARfit function also has methods to extract fitted values and the residuals. For example, Figure 2 shows the first node time series and the residuals from fitting the model. By altering the input parameters in the GNARfit function, we can fit a range of different GNAR models and the reader can consult Appendix C for further examples.

Example: GNAR data simulation on a given network
The following example demonstrates network time series simulation using the network in Figure 1.

Missing data and changing connection weights with GNAR models
Standard multivariate time series models, including vector autoregressions (VAR), can have significant problems in coping with certain types of missingness and imputation is often used, see Guerrero and Gaspar (2010), Honaker and King (2010), Bashir and Wei (2016). While in VAR modeling successful solutions have been found to cope with specific missingness scenarios, such as implemented in the gimme R package (Lane et al. 2020), however, if a variable has e.g., block missing data, the coefficients corresponding that variable can be difficult to calculate, and impossible if their partner variable is missing at cognate times. In addition, due to computational burden gimme is limited to modeling a single time lag. A key advantage of our parsimonious GNAR model is that it models via neighborhoods across the entire data set. If a node is missing for a given time, then it does not contribute to the estimation of neighborhood parameters that the network structure suggests it should, and there are plenty of other nodes that do contribute, generally resulting in a high number of observations to estimate each coefficient. In GNAR models, missing data of this kind is not a problem.
The flexibility of GNAR modeling means that we can also model missing data as a changing network, or alternatively, as changing connection weights. In the situation where the overall network is considered fixed, but when observations are missing at particular nodes, the connections and weightings need altering accordingly. Again, using the graph in Figure 1, consider the situation where node A does not have any data recorded. Yet, we want to preserve the stage-2 connection between D and E, and the stage-3 connection between E and both B and C. To do this, we do not redraw the graph and remove node A and its connections, instead we reweight the connections that depend on node A. As node A does not feature in the stage-2 or stage-3 neighbors of E, the connection weights ω E,D , ω E,B , ω E,C do not change, but the connection weight ω E,A drops to zero in the absence of observation from node A. Similarly, the stage-1 neighbors of D are changed without A, so ω D,A drops to zero and the other two connection weights from node D increase accordingly; ω D,B = ω D,C = 0.5. Missing data of this kind is handled automatically by the GNAR functions using customary NA missing data values present in the vts (vector time series) component of the overall network time series. For example, inducing some (artificial) missingness in the fiveVTS series, we can still obtain estimates of model parameters: R> fiveVTS0 <-fiveVTS R> fiveVTS0[50:150, 3] <-NA R> nafit <-GNARfit(vts = fiveVTS0, net = fiveNet, alphaOrder = 2, + betaOrder = c(1, 1)) R> layout(matrix(c(1, 2), 2, 1)) R> par(mar = c(4.1, 4.1, 0.75, 2.1), cex.axis = 0.75) R> plot(ts(fitted(nafit)[, 3]), ylab = "Node C fitted values") R> par(mar = c(4.1, 4.1, 0.75, 2.1), cex.axis = 0.75) R> plot(ts(fitted(nafit)[, 4]), ylab = "Node D fitted values") As shown in Figure 4, after removing observations from the time series at node C, its neighbor, node D, still has a complete set of fitted values.

Stationarity conditions for a GNAR process with fixed network
Theorem 1 Given an unchanging network, G, a sufficient condition for the GNAR model (1) to be stationary is (2) The proof of Theorem 1 can be found in Appendix A.
For the global-α model this condition reduces to We can explore these conditions using the GNARsim function. The following example uses parameters whose absolute value sums to greater than one and then we calculate the mean over successive time periods. The mean increases rapidly indicating nonstationarity.

Benefits of our model and comparisons to others
Conditioned on a given network fixed in time and with a known (time-dependent) weightand neighborhood structure, the GNAR model can be mathematically formulated as a specific restricted VAR model, where the restrictions are imposed by the network and thus impact model parametrization, as mathematically encoded by Equation 1. This is explored in more depth in Appendix B and contrasts with a VAR model where any restrictions can only be imposed on the parameters themselves.
An unrestricted VAR model with dimension n has O(n 2 ) parameters, whereas a GNAR model with known network (usually) has O(n) parameters, and a global-α GNAR model can have O(1) parameters. The large, and rapidly increasing, number of parameters in VAR often make it a challenging model to fit and non-problem-specific mathematical constraints are often used to mitigate those challenges. Further, the large number of VAR parameters usually mean that it fits multivariate time series well, but then performs poorly in out-of-sample prediction. An example of this is shown in Section 4.
Our model has similarities with the network autoregression introduced by Zhu, Pan, Li, Liu, and Wang (2017), motivated by social networks.
In our notation, the Zhu et al. (2017) model can be written as a special case as where β 0 is a global intercept term, Z i is a vector of node-specific covariates with corresponding parameters γ, ω i is the reciprocal of the out-degree of node i, and the innovations are independent and identically distributed, with zero mean, such that var(u i,t ) = σ 2 . Hence, the Zhu et al. (2017) model without intercept and node-specific covariates is a special case of our GNAR model, with max j∈{1,...,p} s j = 1, i.e., dependencies limited to stage-1 immediate neighbors, and un-weighted edges.
Our model is designed to deal with a time-varying network, and our β j,r,c parameters can include general edge-based covariate information. A further important advantage is that our GNAR model in Section 2.2 can express dependence on stage-r neighbor sets for any r.
An earlier model with similarities to the generic network autoregression is the dynamic Bayesian network (DBN) model considered in Spencer, Hill, and Mukherjee (2015). Their model can be written as where β 0,i is a node-specific intercept term, the other β parameters describe the network autoregression, and u i,t ∼ N (0, σ 2 i ). The DBN model is also a constrained VAR model, but with no univariate autoregression terms, and the network autoregression only includes the stage-1 neighbors. Unlike our model and the Zhu et al. (2017) model, there are no restrictions on the parameters other than parameters only being present when there is an edge between two nodes. The Spencer et al. (2015) framework does not allow for a range of networks, as their underlying network is assumed to be a Directed Acyclic Graph. With these assumptions, the network and parameters are inferred by considering potential predictors for each node in turn. A key difference between our model and the Spencer et al. (2015) model is that we assume that the behavior of connected nodes is the same throughout the network, whereas the DBN model allows for different β parameters for different connections, including allowing a change of sign.
The benefits of the GNAR model compared to these, and other models, include the ability to deal with a time-changing network, missing observations, and using network information to reduce the number of parameters. As detailed in Section 2.6, we can incorporate missing data information with the GNAR model by allowing the connection weights to change. Allowing for a changing network structure enables us to model new nodes being added to the system, or connections between nodes changing over time. Adding autoregressive parameters to neighbors with stage greater than one results in our model being able to capture more network relationships than just those of immediate neighbors.

Estimation
In modeling terms, our GNAR model is a linear model and we employ standard techniques such as least squares estimation to fit them and to provide statistically consistent estimators, as verified in Appendix B. An important practical consideration for fitting GNAR models is the choice of model order. Specifically, how do we select p and s?

Order selection
We use the Bayesian information criterion (BIC) proposed by Schwarz (1978) to select the GNAR model order. Under the assumption of a constant network, and that the innovations are independent and identically distributed white noise with bounded fourth moments, this criterion is consistent, as shown in Lütkepohl (2005). The BIC allows us to select both the lag and neighborhood orders simultaneously by selecting the model with smallest BIC from a set of candidates.
For a general candidate GNAR(p, [s]) model with N nodes, the BIC is given by whereΣ p,s = T −1Û Û ,Û is the residual matrix from the NAR(p, [s]) fit, and M is the number of parameters. In the general case M = N p + C p j=1 s j , and in the global-α model M = p + C p j=1 s j . The covariance matrix estimate,Σ p,s , is also the maximum likelihood estimator of the innovation covariance matrix under the assumption of Gaussian innovations.
GNAR enables us to easily compute the BIC for any model by using the BIC method for 'GNARfit' objects. For example, on the default model fitted by GNARfit, and an alternative model that additionally includes second-order neighbors at the first lag into the model, we can compare their BICs by

R> BIC(GNARfit())
[1] -0.003953124 R> BIC(GNARfit(betaOrder = c(2, 1))) [1] 0.02251406 Whilst we focus on the BIC for model selection for the remainder of this article, the GNAR package also include functionality for the Akaike information criterion (AIC) proposed by Akaike (1973) as AIC(p, s) = ln|Σ p,s | + 2T −1 M, whereΣ p,s is as defined in Equation 3 and M is again the number of model parameters. Similar to above, the AIC can be obtained by using the code R> AIC(GNARfit()) [1] -0.06991947 R> AIC(GNARfit(betaOrder = c(2, 1))) [1] -0.05994387 Similar to the BIC, the model with the lowest AIC is preferred. Note that the likelihood of the data associated to the model fit can also be obtained using e.g., logLik(GNARfit()).
Various models can be tried to obtain a good fit whilst, naturally, attending to the usual aspects of good model fitting, such as residual checks. A thorough simulation study that displays the numerical performance of our proposed method appears in Section 4.5 of Leeming (2019).

Model selection on a wind network time series
GNAR incorporates the data suite vswind that contains a number of R objects pertaining to 721 wind speeds taken at each of 102 weather stations in England and Wales. The suite contains the vector time series vswindts, the associated network vswindnet, a character vector of the weather station location names in vswindnames and coordinates of the stations in the two column matrix vswindcoords. The data originate from the UK Met Office site http://wow.metoffice.gov.uk/ and full details can be found in the vswind help file in the GNAR package. Figure 5 shows a picture of the meteorological station network with distances created by [1] -233.3848 R> BIC(GNARfit(vts = vswindts, net = vswindnet, alphaOrder = 1, + betaOrder = 0, globalalpha = FALSE)) [1] -233.1697 Interestingly, the model with the single α gives the better fit, as judged by BIC. The single α model with alphaOrder = 2 and betaOrder = c(0, 0) gives a lower BIC of −243, so we investigate this next. Note that this model also gives the lowest AIC score. In particular, we explore a set of GNAR(2, [b1, b2]) models with b1, b2 ranging from zero to 14 using the following code: R> BIC.Alpha2.Beta <-matrix(0, ncol = 15, nrow = 15) R> for(b1 in 0:14) + for(b2 in 0:14) + BIC.Alpha2.Beta[b1 + 1, b2 + 1] <-BIC(GNARfit(vts = vswindts, + net = vswindnet, alphaOrder = 2, betaOrder = c(b1, b2))) R> contour(0:14, 0:14, log(251 + BIC.Alpha2.Beta), + xlab = "Lag 1 Neighbor Order", ylab = "Lag 2 Neighbor Order") The results of the BIC evaluation for incorporating different and deeper neighbor sets, at lags one and two, are shown in the contour plot in Figure 6. The minimum value of the BIC occurs in the bottom-right part of the plot, where it seems incorporating five or sixth-stage neighbors for the first time lag is sufficient to achieve the minimum BIC, and incorporating further lag one stages does not reduce the BIC. Moreover, increasing the lag two neighbor sets beyond first stage neighbors would appear to increase the BIC for those lag one neighbor stages greater than five (the horizontal contour at 0 in the bottom right hand corner of the plot). A fit of a possible model is R> goodmod <-GNARfit(vts = vswindts, net = vswindnet, alphaOrder = 2, + betaOrder = c(5, 1)) R> goodmod We investigated models with alphaOrder equal to two, three, four and five, but with no neighbors. As judged by BIC, alphaOrder = 3 gives the best model. We could extend the example above to investigate differing stages of neighbors at time lags one, two and three. However, a more comprehensive BIC investigation would examine all combinations of neighbor sets over a large number of time lags. This would be feasible, but computationally intensive for a single CPU machine, but could be coarse-grain parallelized. Further analysis would proceed with model diagnostic checking and further modeling as necessary.

Constructing a network to aid prediction
Whilst some multivariate time series have actual, and sometimes obvious, networks associated with them, our methodology can be useful for series without a clear or supplied network. We propose a network construction method that uses prediction error, but note here that our scope is not to estimate an underlying network, but merely to find a structure that is useful in the task of prediction. Here, we use a prediction error measure, understood as the sum of squared differences between the observations and the estimates: N i=1 (X i,t −X i,t ) 2 . The predict S3 method for GNAR models takes an input GNARfit model object and from this predicts the nodal time series at the next timepoint, similar to the S3 method for the 'Arima' class. This allows for a "ex-sample" prediction evaluation. The predict function outputs the prediction as a vector. For example, to predict the series at the last timepoint R> prediction <-predict(GNARfit(vts = fiveVTS[1:199,], net = fiveNet, + alphaOrder = 2, betaOrder = c(1, 1)))

R> prediction
Time Series: Start = 1 End = 1 Frequency = 1 Series 1 Series 2 Series 3 Series 4 Series 5 1 -0.6427718 0.2060671 0.2525534 0.1228404 -0.8231921 For a small-dimensional multivariate series, any and all potential un-weighted networks can be constructed and the corresponding prediction errors compared using the predict method. Next, we consider the larger data setting where it is computationally infeasible to investigate all possible networks. Erdős-Rényi random graphs can be generated with N nodes, and a fixed probability of including each edge between these nodes, see Chapter 11 of Grimmett (2010) for further details. The probability parameter controls the overall sparsity of the graph. Many random graphs of this type can be created, and then our GNAR model can be used for within-sample prediction. The prediction error can then be used to identify networks that aid prediction. We give an example of this process in the next section.

OECD GDP: Network structure aids prediction
We obtained the annual gross domestic product (GDP) growth rate time series for 35 countries from the OECD website (OECD 2018). The series covers the years 1961-2013, but not all countries are included from the start. The values are annual growth rates expressed as a percentage change compared to the previous year. We differenced the time series for each country to remove the gross trend.

Finding a network to aid prediction
This section considers the case where we observe data up to t = 51, and then wish to predict the values for each node at t = 52. We begin by exploring "within-sample" prediction at t = 51, and identify a good network for prediction. We use randomly generated Erdős-Rényi graphs using the GNAR function seedToNet. To demonstrate this, the GNAR package contains the gdp data and a set of seed values, seed.nos so that the random graphs can be reproduced for use with the time series object 'gdpVTS' here.
As well as investigating which network works best for prediction, we also need to identify the number of parameters in the GNAR model. Initial analysis of the autocorrelation function For the GDP example, we simulate 10,000 random un-directed networks, each with connection probability 0.15, and predict using the GNAR model with the orders above. Hence, this example requires significant computation time (about 90 minutes on a desktop PC), so only a segment of the analysis is included in the code below. For computational reasons, we first divide through by the standard deviation at each node so that we can model the residuals as having equal variances at each node. The function seedSim outputs the sum of squared differences between the prediction and original values, and we use this as our measure of prediction accuracy.
R> net921 <-seedToNet(seed.no = seed.nos[921], nnodes = 35, + graph.prob = 0.15) R> layout(matrix(c(1), 1, 1)) R> plot(net921, vertex.label = colnames(gdpVTS), vertex.size = 0) The network generated from seed.nos [921] is plotted in Figure 10, where all countries have at least two neighbors, with 97 edges in total. This "921" network was constructed with GDP prediction in mind, so we would not necessarily expect any interpretable structure in our found network (and presumably, there were other networks with not too dissimilar predictive power). However, the USA, Mexico and Canada are extremely well-connected with eight, eight and six edges, respectively. Sweden and Chile are also well-connected, with eight and seven edges, respectively. This might seem surprising, but, e.g., the McKinsey Global Institute MGI Connectedness Index, see Manyika, Lund, Bughin, Woetzel, Stamenov, and Dhingra (2016), ranks Sweden and Chile 18th and 45th respectively out of 139 countries, and each country is most connected within their regional bloc (Nordic and South America, respectively). Each of these edges, or subgraphs of the "921" network could be tested to find a sparser network with a similar predictive performance, but we continue with the full chosen network here.
Using this network, we can select the best GNAR order using the BIC. The model that minimized BIC in this case was the sixth model, GNAR(2, [2, 0]), a model with two autoregressive parameters and network regression parameters on the first two neighbor sets at time lag one.

Results and comparisons
We use the previous section's model to predict the values at t = 52 and compare its prediction errors to those found using standard AR and VAR models. The GNAR predictions are found by fitting a GNAR(2, [2, 0]) model with the chosen network (corresponding to seed.nos [921]) to data up to t = 51, and then predicting values at t = 52. We first normalize the series, and then compute the total squared error from the model fit. The fitted parameters of this GNAR model wereα 1 −0.42,β 1,1 −0.13,β 1,2 0.28, and α 2 −0.33.
We compared our methods with results from fitting an AR model individually to each node using the forecast.ar() and auto.arima() functions from version 8.13 of the CRAN forecast package (Hyndman, Athanasopoulos, Bergmeir, Caceres, Chhay, O'Hara-Wild, Petropoulos, Razbash, Wang, and Yasmeen 2020), for further details see Hyndman and Khandakar (2008). Due to our autocorrelation analysis from Section 4.1 we set the maximum AR order for each of the 35 individual models to be p = 2. Conditional on this, the actual order selected was chosen using the BIC. [1] 8.065491 Our VAR comparison was calculated using version 1.5-3 of the CRAN package vars (Pfaff 2008). The missing values at the beginning of the series cannot be handled with current software, so are set to zero. The number of parameters in a zero-mean VAR(p) model is of order pN 2 . In this particular example, the dimension of the observation data matrix is T × N , with T < 2N , so only a first-order VAR can be fitted. We fit the model using the VAR function and then use the restrict function to reduce dimensionality further, by setting to zero any coefficient whose associated absolute t statistic value is less than two. R> library("vars") R> gdpVTSn2.0 <-gdpVTSn2 R> gdpVTSn2.0[is.na(gdpVTSn2.0)] <-0 R> varforecast ], p = 1, + type = "none")), n.ahead = 1) This results in forecast vectors for each node, so we extract the point forecast (the first element of the forecast vectors) and compute the prediction error as follows Our GNAR model gives a lower prediction error than both the AR and VAR results, reducing the error by 29% compared to AR and by 78% compared to VAR. Table 1 summarizes these results and also shows the number of parameters fitted. It is clear that GNAR is particularly parsimonious.

Model
Prediction error at t = 51 Prediction error at t = 52 GNAR(2, [0, 0]) 11.8 8.1 Individual AR (2) 18.6 11.3 VAR (1) 115.0 120.4 [1] 120.4467 Table 2 shows that the GNAR model is again the best performing, although in the two-step ahead prediction the fitted model is a special case of GNAR model with no neighborhood parameters. Tables 1 and 2 indicate that the VAR model works particularly poorly here, despite using thresholding to reduce the number of parameters. This example highlights that, for a multivariate series with many observations per time point, the VAR framework is restricted by the number of parameters that have to be fitted per time lag, thus reducing the AR-order, p, it can capture. In addition, we were unable to find software to fit VAR models with for missing data at the start of a series.

Results in
We end this section by noting that using Erdős-Rényi graphs are not the only type of network that could be used to aid prediction. As suggested by a referee, models such the Chung-Lu model (Aiello, Chung, and Lu 2001;Chung and Lu 2002) could also be used to simulate random networks for this task; these graphs would allow for more flexible network generation, for example using node-specific connection probabilities proportional to a country's size.

Discussion and summary
The GNAR package can be used to model network time series using a network autoregressive structure. Estimation under the proposed model is informed by the, potentially time-varying, structure of the network, assumed known. Network time series models are in an early stage of development, but there is enormous potential, especially as network data are increasingly being collected and analyzed in many fields. As far as possible, we attempt to integrate our methods with existing valuable R functionality, such as its linear modeling capability and the fit/summary/predict methods that are familiar with R users.
Within our model a network is formed using edges of all covariates simultaneously, and the connection weights of this single network can be calculated e.g., as described in Section 2.1.2. Another approach is to consider a separate network for each covariate, and then calculate connection weights for each of these networks. This would result in different (known) weightings, ω, and consequently different fitted coefficients, β. The single network approach is more appropriate for sparse networks and when different types of edge are closely related. In comparison, when covariates relate completely separate link information between the nodes, use of different networks would be appropriate.
When covariates are present, the neighbor set structure is more complex, as different edge types can be included in a path between nodes. For example, in a network with event and proximal edges, network paths between stage-2 neighbors could include edges event-event, event-proximal/proximal-event, or proximal-proximal. These different types of path could be represented separately in the model using additional β parameters. We note that the number of such parameters would increase greatly for large covariate cardinality C or high neighbor set stage s j , so, in these cases, the large number of additional parameters may not enhance the model. Our model permits regression on any non-empty stage neighbor set, so models with high s j can be fitted. For large s j , the neighbor sets may not be scientifically interpretable so small s j is recommended, to favor parsimony and interpretability.
Trend is another factor that can seriously affect modeling and estimation, just as in the regular time series situation. However, trend can be successfully modelled and estimated by using second-generation wavelet (lifting) techniques before stochastic modeling, as in Nunes, Knight, and Nason (2015).
With the option of having different covariates and high order neighborhood structures included, our GNAR model as presented in Section 2 is incredibly flexible. In this article a sufficient condition for stationarity and consistency of the fitted parameters have been shown for the fixed network scenario. In addition, practical suggestions for order selection, and connection weights in the case of missing data have been discussed.

A. Proof of stationarity conditions for the GNAR model
A sufficient condition for stationarity of the GNAR model (1) with a static network is |β j,r,c | < 1 ∀i ∈ 1, . . . , N.
Proof: First Gerschgorin's theorem and a corollary are presented without proof, both taken from Varga (1962).
Theorem Let A = (a i,j ) be an arbitrary n × n complex matrix, and let Λ i ≡ n j=1,j =i |a i,j |, Then, all of the eigenvalues λ of A lie in the union of the disks |z − a i,i | ≤ Λ i , 1 ≤ i ≤ n.
Since the disk |z − a i,i | ≤ Λ i is a subset of the disk |z| ≤ |a i,i | + Λ i , we have the immediate result of Corollary 1 If A = (a i,j ) is an arbitrary n×n complex matrix with eigenvalues λ i , 1 ≤ i ≤ n, We can write the static-network GNAR process X t = (X 1,t , . . . , X N,t ) as a VAR process, by writing X t = φ 1 X t−1 + · · · + φ p X t−p + u t , where φ k are n × n matrices such that r,c) , where matrices W (r,c) have entries [W (r,c) ] ,m = ω ,m,c I{m ∈ N (r) ( )} and u t is the vector of errors at time t. We use the notation [·] ,m to denote the , m entry of a matrix.
From Brockwell and Davis (2006), for example, we have that if det(I N − φ 1 z − · · · − φ p z p ) = 0, for all z ∈ C such that |z| ≤ 1, then the VAR model has exactly one stationary solution. Using Lemma 2.1 from Tsay (2014) we have det(I N − φ 1 z − · · · − φ p z p ) = det(I N p − Φz), where Φ is the N p × N p companion matrix defined as where I N and 0 N are the N × N identity and zero matrices, respectively. 1 Thus we require that the roots of det(I N p − Φz) are outside of the unit circle for stationarity, or equivalently, that the eigenvalues of Φ lie inside the unit circle.
We investigate the eigenvalues of Φ using Corollary 1.
For |Φ ,m | ≤ 1, and, using Corollary 1, we have that the spectral radius of Φ is at most one.
We next check whether an eigenvalue with modulus 1 is possible.
Therefore v k = λ p−k v p ∀k ∈ {1, . . . , p} and, replacing this on both sides of (i), we have that p k=1 φ k λ p−k v p = λ p v p . This results in the equation p k=1 φ k λ −k v p = v p , which can be written in matrix form as Ψ v p = v p , where Ψ is the N × N matrix Ψ = p k=1 φ k λ −k . Hence, if Φ has an eigenvalue of modulus 1, then Ψ must have 1 as an eigenvalue.
We again use Corollary 1 for the eigenvalues of Ψ , under the assumption |λ| = 1. Under condition (2) this is smaller than 1, so, by Corollary 1, no eigenvalues of Ψ have modulus 1 or greater. This contradicts the assumption that an eigenvalue of Φ, λ, exists such that |λ| = 1. Hence, the eigenvalues of Φ are inside the unit circle under condition (2) and the GNAR model is stationary.

B. Parameter estimate consistency
We employ least squares estimation for the GNAR model parameters and establish their consistency using results from Lütkepohl (2005). The column form of the static-network GNAR(p, [s]) model can be written in a VAR framework as where the matrices φ i contain the network information. In matrix form the GNAR model is X = BZ + U , where X = [X p+1 , . . . , X T ], B = [φ 1 , . . . , φ p ], Z = [Z p , . . . , Z T −1 ], with Z t = [X t , . . . , X t−p+1 ], and U = [u p+1 , . . . , u T ]. The constraints imposed to form a GNAR model can be written linearly as vec (B) = Rγ, where R is the constraint matrix embedding the network structure of dimension pN 2 × M , γ is an unrestricted parameter vector of length M , where M is defined as in Section 3.1, and vec is the operator that stacks the columns of a matrix into a vector. Using the estimated generalized least squares estimator, we apply results from Section 5.2 of Lütkepohl (2005) to obtain consistency for the GNAR parameters. Let ⊗ denote the Kronecker product and plim denote limit in probability.
Proposition 1 Suppose {X t } is an N -dimensional, stationary GNAR(p) process with a static network, whose innovations {u t } are independent white noise with finite fourth moment, and covariance matrix Σ u . Then, given an estimator of the innovation covariance matrixΣ u , such that plimΣ u = Σ u , the estimated generalized least squares estimator of the unrestricted parameters, γ = {R (ZZ ⊗Σ −1 u )R} −1 R(Z ⊗Σ −1 u ) vec(X ),