Generalised Network Autoregressive Processes and the GNAR package

This article introduces the GNAR package, which fits, predicts, and simulates from a powerful new class of generalised 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 neighbouring variables, with inclusion controlled by the network structure. The GNAR package is designed to fit 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 summarise the data generating process.
We consider time series observations recorded at different nodes of a network, or graph. Our GNAR package (Leeming, Nason, Nunes, and Knight 2019) and its novel generalised network autoregressive (GNAR) statistical models focus on partnering a network with a multivariate time series and modelling 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.
Networks can provide strong information about the dependencies between variables. Within our generalised network autoregressive (GNAR) model, each node depends on its previous values as in the univariate autoregressive framework, but also may depend on the previous values at its neighbours, neighbours of neighbours, and so on. Our GNAR modelling framework is flexible, allowing for different types of network, networks that change their structure over time (time-varying networks), and also can be powerfully applied in the important practical situation where the time series feature missing observations. Driven in part by the increased popularity and recent research activity in the field of statistical network analysis, there has been a concurrent growth in software for analysing such data. An exhaustive list of these packages is beyond the scope of this article, but we review some relevant ones here.
Existing software in this area predominantly focusses on the various models for networkstructured data. In the static network setting, these include packages dedicated to latent space network models, such as collpcm (Wyse, Ryan, and Friel 2017), HLSM (Adhikari, Junker, Sweet, and Thomas 2018), latentnet (Krvitsky, Handcock, Shortreed, Tantrum et al. 2018b) amongst others; exponential random graph models and their variants, for example ergm (Handcock, Hunter, Butts et al. 2018), GERGM (Denny, Wilson, Cranmer, Desmarais, and Bhamidi 2018) or hergm (Schweinberger, Handcock, and Luna 2018); and block models in e.g., blockmodels (Leger 2015). For dynamic networks, packages for time-varying equivalents of these network models are also available, see e.g., the tergm package (Krvitsky, Handcock, Hunter, Goodreau et al. 2018a) or dynsbm (Matias and Miele 2018). There are also a multitude of more general packages for network analysis, e.g., for network summary computation or implementations of methodology in specific applications of interest.
Despite this, software dedicated to the analysis of time series and other processes on networks is sparse. A number of packages implement epidemic (e.g., SIR) models of disease spread, notably epinet (Groendyke, Welch, and Hunter 2018), EpiLM/EpiLMCT Almutiry, Warriyar, and Deardon 2018) and hybridModels (Marquez, Grisi-Filho, and Amaku 2018); these use transmission rates to model processes as opposed to temporal and network dependence through time series models as in GNAR. Similarly, the NetOrigin software (Manitz and Harbering 2018) is dedicated to source estimation for propagation processes on networks, rather than fitting time series models. Packages such as networkTomography (Blocker, Koullick, and Airoldi 2014) deal with time-varying models for (discrete) count processes or flows on links of a fixed routing network; the tnam package (Leifeld and Cranmer 2017) fits models using spatial (and not network-node) dependence. Both of these are in contrast to the GNAR package, which implements time series models which account for known time-varying network structures.
Other packages can implicitly develop network-like structured time series models through penalised or constrained variable selection, such as autovarCore (Emerencia 2018), nets (Brownlees 2017), sparsevar (Vazoller, Frattarolo, and Billio 2016), as well as the vars package (Pfaff 2008). Packages that take a graphical modelling approach to the dependence structure within time series include gimme (Lane, Gates, Fisher, Molenaar et al. 2019), graphicalVAR (Epskamp 2018), mgm (Haslbeck 2019), mlVAR (Epskamp, Deserno, and Bringmann 2019), and sparseTSCGM (Abegaz and Wit 2016). These approaches also differ fundamentally from the GNAR models since the network is constructed during analysis, as opposed to GNAR, which specifically incorporates information on the network structure into the model a priori. The vars package features in Section 4.2, where we highlight the differences between the GNAR models and this existing class of techniques.
Section 2 introduces our model, and demonstrates how GNAR can be used to fit network models to simulated network time series in Section 2.4. Order selection and prediction are discussed in Section 3, which includes an example of how to use BIC to select model order for a wind speed network time series in Section 3.2. An extended example, concerning constructing a network to aid GDP forecasting, is presented in Section 4. Section 5 discusses different network modelling options that could be chosen, and presents a summary of the article. All results were calculated using version 3.5.1 of the statistical software R (R Core Team (2019)).

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 neighbouring nodes at previous time steps. These neighbouring 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 neighbourhoods
We introduce the notion of neighbours and stage-neighbours in the graph structure as follows; for a subset A ⊂ K the neighbour set of A is given by N (A) = {j ∈ K/A : i j; i ∈ A}. These are the first neighbours, or stage-1 neighbours of A. The rth stage neighbours of a node i ∈ K are given by N (r) , for r = 2, 3, ... and N (1) (i) = N ({i}). Figure 1 shows an example graph, where node E has stage-1 neighbour A, stage-2 neighbour D, and stage-3 neighbours B and C. Neighbour 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 neighbour 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 neighbour set and also encodes any edge-weight information. Formally, the values of the connection weights from a node i ∈ K to its stage-r neighbour j ∈ N (r) (i) will be the reciprocal of the number of stage-r neighbours; ω 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-normalised weight between these nodes, denote this µ i, ∈ R + . To find the length of connection between a node i and its stage-r neighbour, 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 neighbour). 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 neighbour 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 neighbour 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 neighbour 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 generalised 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 generalised 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 neighbour dependence for time lag j, 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 neighbours, 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 parametrisation through the use of time-varying weights and neighbours. 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 behaviour 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 neighbouring 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 neighbours 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> 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.

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") 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.

Example: GNAR data simulation on a given network
The following example demonstrates network time series simulation using the network in Figure 1.
Both simulations are created using standard normal noise whose standard deviation is controlled using the sigma argument.

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 modelling successful solutions have been found to cope with specific missingness scenarios, such as implemented in the gimme R package (Lane et al. 2019), 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 modelling a single time lag. A key advantage of our parsimonious GNAR model is that it models via neighbourhoods across the entire data set. If a node is missing for a given time, then it does not contribute to the estimation of neighbourhood 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 modelling 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 neighbours 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 neighbours 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.

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 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) weight-and neighbourhood 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 parametrisation, 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 neighbours, 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 neighbour 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 neighbours. 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 behaviour 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 neighbours with stage greater than one results in our model being able to capture more network relationships than just those of immediate neighbours.

Estimation
In modelling 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 neighbourhood 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 neighbours 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 (6) 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). 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:

Model selection on a wind network time series
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 Neighbour Order", ylab = "Lag 2 Neighbour Order") The results of the BIC evaluation for incorporating different and deeper neighbour sets, at lags one and two, are shown in the contour plot in Figure 6. occurs in the bottom-right part of the plot, where it seems incorporating five or sixth-stage neighbours 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 neighbour sets beyond first stage neighbours would appear to increase the BIC for those lag one neighbour 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 neighbours. As judged by BIC, alphaOrder = 3 gives the best model. We could extend the example above to investigate differing stages of neighbours at time lags one, two and three. However, a more comprehensive BIC investigation would examine all combinations of neighbour 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 modelling 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 ], net = fiveNet, + alphaOrder = 2, betaOrder = c(1, 1))) R> prediction 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 1 . 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.
We use the first T = 52 time points and designate each of the 35 countries as nodes to investigate the potential of modelling this time series using a network. In this data set 20.8% (379 out of 1820) of the observations were missing due to some nodes not being included from the start. We model this by changing the network connection weights as described in Section 2.6. In this example, we do not use covariate information, so C = 1. The pattern of missing data along with the time series values is shown graphically in Figure 7, produced by the following code.

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. R> layout(matrix(c(2, 1), 1, 2)) R> par(mar=c(0,1,0,1)) R> plot(net1, vertex.label = colnames(gdpVTS), vertex.size = 0) R> plot(net2, vertex.label = colnames(gdpVTS), vertex.size = 0) Figure 8 shows two of these random graphs.
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 at each node indicated that a second-order autoregressive component should be sufficient, so GNAR models with orders up to p = 2 were tested, and we included at most two neighbour sets at each time lag. The GNAR models are: 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 Erdős-Rényi random graphs constructed from the first two elements of the seed.nos variable with 35 nodes and connection probability 0.15. 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 neighbours, 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 minimised 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 neighbour 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 normalise 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.0 of the CRAN forecast package (Hyndman, Athanasopoulos, Christoph Bergmeir, Caceres et al. 2017), 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-2 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 summarises these results and also shows the number of parameters fitted. It is clear that GNAR is particularly parsimonious.
We repeat the procedure above to perform analysis based upon two-step ahead forecasting. In this case, a different network minimises the prediction error for model GNAR(2,[2,2] [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 neighbourhood 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,

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 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.
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 analysed in many fields. As far as possible, we attempt to integrate our methods with existing valuable R functionality, such as its linear modelling 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 neighbour 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 neighbours 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 neighbour 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 neighbour set, so models with high s j can be fitted. For large s j , the neighbour sets may not be scientifically interpretable so small s j is recommended, to favour parsimony and interpretability.
Trend is another factor that can seriously affect modelling 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 modelling, as in Nunes, Knight, and Nason (2015).
With the option of having different covariates and high order neighbourhood 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. 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 W (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( where I N and 0 N are the N × N identity and zero matrices, respectively. 2 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. |Φ ,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 generalised 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 generalised least squares estimator of the unrestricted parameters, Again, adapting Lütkepohl (2005), we have the following result for a consistent estimator of the innovation covariance matrix in the GNAR setting.

Proposition 2 A consistent estimator of Σ u is given bỹ
whereBZ are the fitted values from estimating the parameters using the least squares estimator γ = {R (ZZ ⊗ I N )R} −1 R (Z ⊗ I N ) vec(X ).
Estimating the parameters withγ involves using the linear constraints, but assumes independent and identically distributed innovations across nodes.