\section{GNAR Model and Methods}
\label{sec:gnar methodology}
%\addtolength{\textheight}{.5in}%
A network time series $\mathcal{X} := (\boldsymbol{X}_t, \mathcal{G})$ is a stochastic process composed of a multivariate time series $\boldsymbol{X}_t \in \mathbb{R}^d$,
$t=1, \ldots, T$, and an underlying network $\mathcal{G} = (\mathcal{K}, \mathcal{E})$, where $\mathcal{K} = \{1, \dots, d\}$ is a set of nodes, $\mathcal{E}  \subseteq \mathcal{K} \times \mathcal{K}$ is a set of edges, and $\mathcal{G}$ is an undirected graph, which has $d$ nodes
(vertices). Each univariate time series $X_{i, t} \in \mathbb{R}$ is associated to node $i \in \mathcal{K}$ in $\mathcal{G}$.
\par 
We review the global--$\alpha$ GNAR model for analysing network time series given our particular focus on similar interactions across nodes in the network. We assume that the network is static throughout, 
although, GNAR processes can handle time--varying networks; see \cite{gnar_org, gnar_paper}. We next
introduce two matrices that enable to express the GNAR model more compactly.
\subsection{Weights and \texorpdfstring{$r$}{r}--stage adjacency matrices}
GNAR processes are parsimonious models as they exploit network structure. A key notion is that of $\rstage$ neighbours. We say that nodes $i$ and $j$ are $\rstage$ neighbours if and only if the shortest path on the network $\mathcal{G}$ between them has `distance' equal to $r$, in the sense that the number of edges on the shortest path is equal to $r$, and define $\mathcal{N}_r (i) \subseteq \mathcal{K}$ as the set of $\rstage$ neighbours of node $i$.
We denote this higher--order structure by introducing the $\rstage$ adjacency matrix $\bmat{S}_r$.
%The $\bmat{S}_r$ matrix indicates whether or not nodes $i$ and $j$ are $\rstage$ neighbors, hence, each entry $[\mathbf{S}_r]_{ij}$ either equals one or zero, where $[\bmat{S}_r]_{ij} = 1$ to signal that $j \in \mathcal{N}_r (i)$, since this relationship is symmetric we also have that $[\bmat{S}_r]_{ji} = 1$.
\begin{definition} \label{def: rstage adjacency}
    Let $\delta(i, j)$ be the number of edges between $i$ and $j$ on the shortest path between them, and $(\boldsymbol{X}_{t}, \mathcal{G})$ be a network time series as above.
    Define the \mbox{$r$--stage} adjacency matrix $\mathbf{S}_r$ $\in \mathbb{R}^{d \times d}$ to
    have entries $[\mathbf{S}_r]_{ij}$, where each entry $[\mathbf{S}_r]_{ij} = 1$ if and only if node $i$ is an \mbox{$r$--stage} neighbour of node $j$, otherwise $[\mathbf{S}_r]_{ij} = 0$. That is, $j \in \mathcal{N}_r (i) \iff \delta(i,j) = r \iff [\mathbf{S}_r]_{ij} = 1$.
\end{definition}
The matrix $\bmat{S}_1$ is the regular adjacency matrix and  nodes cannot be $\rstage$ neighbours for different choices of $r$.
Figure~\ref{fig: rstage neighbors} illustrates $\mathbf{S}_r$ by plotting the set of $\rstage$ neighbours for $r = 1, 6$ present in the network $\mathcal{G}$ associated to the COVID--19 (network) time series. There are no neighbours at seven stages or higher for this network time series.
\begin{figure}
        \centering
        \begin{subfigure}{0.55\textwidth}
            \includegraphics[width=\textwidth]{images/1_neighbors.pdf}
            \caption{NHS Trusts network 1-stage neighbours.}
            \label{fig: 1 stage neighbors}
        \end{subfigure}
        
        \begin{subfigure}{0.55\textwidth}
            \includegraphics[width=\textwidth]{images/6_neighbors.pdf}
            \caption{NHS Trusts network 6-stage neighbours.}
            \label{fig: 6 stage neighbors}
        \end{subfigure}
        \caption{(a) network with directly connected nodes. (b) network corresponding to 6-stage neighbours. NHS Foundation Trust Codes: REF=Royal Cornwall Hospitals;
        RVW=North Tees \& Hartlepool; RXP=County Durham \& Darlington;
        RNN=North Cumbria Integrated Care; RR7=Gateshead Health;
        RTR=South Tees Hospitals; RTF=Northumbria Healthcare;
        R0B=South Tyneside \& Sunderland; RTD=Newcastle upon Tyne Hospitals.}
        \label{fig: rstage neighbors}
    \end{figure}
The set of $\rstage$ neighbours of node $i$ can be found by looking at the $i$th row of $\bmat{S}_r$, furthermore, each $\bmat{S}_r$ is a symmetric matrix that can be computed sequentially from previous $\rstage$ adjacency matrices; see the separate supplementary material document section~C for a more thorough exposition. 
\par
We also consider  weights between nodes $i$ and $j$: each weight $w_{ij} \in [0, 1]$ quantifies the relevance node $j$ has on node $i$ with respect to neighbourhood regression. If the network $\mathcal{G}$ does not have weighted edges, a GNAR model assigns equal importance to each node $j$ in the set of $\rstage$ neighbours in the sense that if node $j$ is an $\rstage$ neighbour of node $i$, then $w_{ij} = \{ | \mathcal{N}_{r} (i) | \}^{-1} $. 
\par
If $\mathcal{G}$ has weighted edges $\tilde{w}_{ij}$, then each weight is normalised, so that the weights associated to each $\rstage$ neighbourhood sum to one, as follows $w_{ij} = \tilde{w}_{ij} \{ \sum_{l \in \mathcal{N}_r (i)} \tilde{w}_{il} \}^{-1}$. Hence, $\sum_{j \in \mathcal{N}_r (i)} w_{ij} = 1$ for all $\rstage$ neighbourhood sets. The connection between these weights and inverse distances is explored in Section \ref{sub: subsect oracle selection}.


A GNAR model performs autoregression for each nodal time series as well as neighbourhood regression. By neighbourhood regression we mean that for each nodal time series $X_{i,t}$ we compute its autoregression at lag $k$ with a convex linear combination $Z_{i, t}^{r}$ of its $\rstage$ neighbours, which for a fixed $r \in \{1, \dots, r_{\max}\}$ is given by
\begin{equation} \label{eq: node neigh regression}
    Z_{i, t}^{r} := \sum_{j \in \mathcal{N}_r (i)} w_{ij} X_{j, t},
\end{equation}
where $r_{\max} \in \mathbb{N}$ is the longest shortest path in the network $\mathcal{G}$,
i.e.\ $\delta(i, j) \leq r_{\max}$ for all node pairs $(i, j)$,
$i=1, \ldots, d$ and $t=1, \ldots, T$.
Since node $j$ can only be part of one neighbourhood regression with respect to all other nodes in $\mathcal{G}$ we can identify the unique connection weight between nodes $i$ and $j$ as $w_{ij}$ for all nodes in $\mathcal{G}$. To help express neighbourhood regression in matrix notation we introduce the weights matrix $\bmat{W}$ for GNAR models. 
\begin{definition} \label{def: gnar weight matrix}
     Let $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ be a network, then the weights matrix is the matrix $\mathbf{W}$ $\in \mathbb{R}^{d \times d}$ with entries $[\mathbf{W}]_{ij} := w_{ij}$ for $i,j \in \mathcal{K}$.
\end{definition}
Note that the diagonal entries are equal to zero because there are no self--loops in $\mathcal{G}$, and that $\bmat{W}$ is not necessarily symmetric because nodes can have different degrees of relevance across the network.
\subsection{GNAR Model}
GNAR models were introduced by \cite{gnar_org}. We review them here and
additionally 
introduce the use of a Hadamard operator, which makes the specification more compact and,
importantly, permits us to write GNAR as a vector process.

Let $\odot$ denote the Hadamard (component--wise) product, which paired with $\bmat{S}_r$ and $\bmat{W}$ can be used to select the set of $\rstage$ neighbours for each node $i$ and compute the corresponding neighbourhood regression $Z_{i, t}^{r}$. The {\em node--wise} representation of a global-$\alpha$ GNAR model with maximum lag equal to $p \in \mathbb{N}$ and maximum $\rstage$ depth $s_k \in \{1, \dots, r_{\max} \}$ at each lag is 
\begin{equation} \label{eq: node-wise}
    X_{i, t} = \sum_{k = 1}^{p} \big ( \alpha_k X_{i, t - k} + \sum_{r = 1}^{s_k} \beta_{kr} Z_{i, t - k}^{r} \big ) + u_{i, t},
\end{equation}
where $i=1, \ldots, d$, $t= 1, \ldots, T$, the $\alpha_k \in \mathbb{R}$ are `standard' autoregressive parameters and the $\beta_{kr} \in \mathbb{R}$ are neighbourhood autoregressive parameters for $r = 1, \dots, s_k$ at each lag $k = 1, \dots, p$. The model is denoted
$\GNAR (p, [s_1, \dots, s_p] )$, which has $p$ autoregressive terms and for each one of these there are $s_k$ neighbourhood regression terms
given by (\ref{eq: node neigh regression}).
The following notational convenience is used
when consecutive $s_k$ are identical in the
$[s_1, \dots, s_p]$ vector. For example, we write
 $[3^{(2)} 7^{(5)}]$ for $[3, 3, 7, 7, 7, 7, 7]$ or $[3^{(4)}]$ for
$[3, 3, 3, 3]$, {\em etc}. The latter concise specification
is useful when describing an initial full model before, e.g., embarking
on backward variable deletion.

Also, we assume that the $u_{i, t}$ is independent and identically distributed (IID) white noise with mean zero and variance $\sigma^2 > 0$ for all nodes. The model given by (\ref{eq: node-wise}) is a more compact representation  of the model in \cite{gnar_org}.
\par
Before expressing the model in (\ref{eq: node-wise}) by a vector-wise representation we introduce the $\rstage$ linear regression vector time series 
\begin{equation} \label{eq: vector projection}
    \boldsymbol{Z}_{t}^{r} := (\bmat{W} \odot \bmat{S}_r) \boldsymbol{X}_{t}.
\end{equation}
\par
Each $i$th entry in $\boldsymbol{Z}_{t}^{r}$ is equal to the $\rstage$ neighbourhood regression for node $i$ given by (\ref{eq: node neigh regression}). Now we can write the vector notation version of the model in (\ref{eq: node-wise}) as
\begin{equation} \label{eq: vector wise}
    \boldsymbol{X}_{t} = \sum_{k = 1}^{p} \big ( \alpha_k \boldsymbol{X}_{t - k} + \sum_{r = 1}^{s_k} \beta_{kr} \boldsymbol{Z}_{t - k}^{r} \big ) 
    + \boldsymbol{u}_{t}, 
\end{equation}
where $t=1, \ldots, T$, the $\alpha_k \in \mathbb{R}$ and $\beta_{kr} \in \mathbb{R}$ are the autoregressive coefficients in (\ref{eq: node-wise}) and $\boldsymbol{u}_t$ are (IID) multivariate white noise with mean zero and covariance matrix $\sigma^2 \bmat{I_d}$. The representation  (\ref{eq: vector wise}) highlights the parsimonious structure of a global-$\alpha$ GNAR model. Note that the number of parameters in the model increases not with the dimension of $\boldsymbol{X}_t$ and lag, but rather with the depth of $\rstage$ regression and maximum lag. 
\par
A close look at (\ref{eq: vector wise}) reveals that the \mbox{global-$\alpha$} GNAR model can be written as a constrained 
Vector Autoregressive (VAR) model; see \cite{brockwell_davis}%cite VAR model
, for which the autoregressive matrices 
$$\bmat{\Phi}_k := \{\alpha_k \bmat{I_d} + \sum_{r = 1}^{s_k} \beta_{kr} (\bmat{W} \odot \bmat{S}_r ) \}, $$
are restricted by the network structure. \cite{gnar_paper} exploit this connection to show a more general result with respect to stationarity conditions for GNAR processes with a static network than the one we give below.

\begin{theorem}{[\textbf{\cite{gnar_paper}}]} \label{th: stationary cond}
Let $\boldsymbol{X}_t$ be a \mbox{global-$\alpha$} $\GNAR (p, [s_1, \dots, s_p] )$ process with associated static network $\mathcal{G} = (\mathcal{K}, \mathcal{E})$. If the autoregressive coefficients in (\ref{eq: vector wise}) satisfy 
$$
    \sum_{k = 1}^{p} ( |\alpha_k| + \sum_{r = 1}^{s_k} |\beta_{kr} | ) < 1, 
$$
then $\boldsymbol{X}_t$ is stationary.
\end{theorem}
%ADD stationarity conditions as theorem proved in previous GNAR paper
We end by highlighting that GNAR models can be viewed as
highly parsimonious VAR models with parameter
constraints
informed by the underlying network structure. If a GNAR formulation is appropriate for modelling a particular data set, then it can be an extremely powerful forecasting tool,
as previously noted \citep{gnar_paper,NasonWei22}, which enables interpretation of a large number of interactions. We describe these properties in Section \ref{sec: graph model connection} and illustrate its advantages by analysing the COVID--19 (network) time series in Section \ref{sec: application}.

\subsection{GNAR Model Estimation}
The original GNAR work \citep{gnar_org} estimated parameters using least squares. Expression as a formal linear model was established by
\cite{leeming} and soon implemented by~\cite{gnar_package}. To introduce notation
and to enable us to properly describe our new connection between GNAR modelling
and variable shrinkage and selection in Section~\ref{sec: graph model connection},
we present a slightly modified and expanded
description here.

A global--$\alpha$ $\GNAR (p, [s_1, \dots, s_p] )$ model has $p$ autoregressive $\alpha_k$ coefficients and $\sum_{k = 1}^{p} s_k$ neighbourhood regression $\beta_{kr} $ coefficients. For comparison, a VAR($p$) model has $p d^2$ parameters, so as long as $q := p + \sum_{k = 1}^{p} s_k < p d^2$ the GNAR model (\ref{eq: vector wise}) needs to estimate far fewer parameters. This is particularly important for settings in which $p d^2$ is larger than the number of time step realisations observed, such as the COVID--19 data we analyse below.
\par
Assume that $T \in \mathbb{N}$ time steps of data $\bmat{X} := [\boldsymbol{X}_{1}, \dots, \boldsymbol{X}_{T} ]$ arising from a GNAR process with known order become available. Our objective is to estimate the unknown autoregressive coefficients $\alpha_k$ and $\beta_{kr}$ after fixing the lag--depth pair $(p, [s_1, \dots, s_k] )$. 
\par
To do this we assume that $\boldsymbol{X}_{t}$ is a stationary $\GNAR (p, [s_1, \dots, s_p] )$ process, and that $n = T - p$ is the number of observations for which there are $p$ previous observed lags available for estimation, and, for $t = p + 1, \dots, T$, define the following.
\begin{align} \label{eq: estimation build-up}
    \boldsymbol{y}_t &:= \boldsymbol{X}_{p+1}, \nonumber \\
    \bmat{Z}^{1:s_1}_{t -1} &:= [\boldsymbol{Z}_{t - 1}^{1}, \dots, \boldsymbol{Z}_{t - 1}^{s_1}], \nonumber \\
    \bmat{R}_{t} &:= [\boldsymbol{X}_{t - 1} , \bmat{Z}^{1:s_1}_{t -1}, \dots, \boldsymbol{X}_{t - p} , \bmat{Z}^{1:s_p}_{t -p}],
\end{align}
where $\boldsymbol{y}_t \in \mathbb{R}^{d}$ is the data vector of `responses', $\bmat{R}_t \in \mathbb{R}^{d \times q}$ is the design matrix at time-step $t$, and $\boldsymbol{Z}_{t - 1}^{s_1}$ are given by (\ref{eq: vector projection}). Furthermore, define the vector of parameters $\boldsymbol{\theta} \in \mathbb{R}^q$ as
$$\boldsymbol{\theta} := (\alpha_1, \beta_{11}, \dots, \beta_{1 s_1}, \alpha_2, \dots, \beta_{p s_p}), $$
which are the unknown linear model coefficients.
\par
Thus, if $\boldsymbol{\theta}$ satisfies the assumptions in Theorem \ref{th: stationary cond}, the realisations
$\boldsymbol{y}_t = \bmat{R}_{t}\, \boldsymbol{\theta} + \boldsymbol{u}_{t}$ are statistically uncorrelated observations of linear models given by \eqref{eq: vector wise}, where $\boldsymbol{u}_t$ are the same as in (\ref{eq: vector wise}) for $t = p + 1, \dots, T$. Next, by concatenating both the column vectors $\boldsymbol{y}_{t}$ into one column vector $\boldsymbol{y} \in \mathbb{R}^{nd}$ and the design matrices $\bmat{R}_{t}$ into one design matrix $\bmat{R} \in \mathbb{R}^{nd \times q}$, we can write (\ref{eq: vector wise}) as the linear model
\begin{equation} \label{eq: T steps linear model}
    \boldsymbol{y} = \bmat{R}\, \boldsymbol{\theta} + \boldsymbol{u}.
\end{equation}
Now, we can couple the $n$ linear regression problems as if we had $n$ independent samples of size $d$ and fit the linear model by ordinary least squares for data with sample size equal to $nd$ and $q$ unknown parameters, thus, the least squares estimator for a $\GNAR (p, [s_1, \dots, s_p] )$ model with a \mbox{global-$\alpha$} specification is given by
\begin{equation} \label{eq: ls estimator}
    \hat{\boldsymbol{\theta}} = \big ( \myprod{\bmat{R}}\,{\bmat{R}} \big )^{-1} \myprod{\bmat{R}}{\boldsymbol{y}}, 
\end{equation}
where $\bmat{R}$ and $\boldsymbol{y}$ are given by (\ref{eq: estimation build-up}).
\par
If we further assume that the $\boldsymbol{u}_t \sim N(0, \sigma^2 \bmat{I_d})$, then $\Hat{\boldsymbol{\theta}}$ given by (\ref{eq: ls estimator}) is also the  conditional  maximum likelihood estimator. We also point out that assuming that all the nodal white noise processes have the same variance might not be sensible, in that case it is possible to adapt (\ref{eq: T steps linear model}) into a generalised least squares problem, and other relaxations
are possible. Nevertheless, throughout this work we maintain the assumptions in Theorem \ref{th: stationary cond} and (\ref{eq: vector wise}), and estimate the unknown coefficients with $\Hat{\boldsymbol{\theta}}$ given by~\eqref{eq: ls estimator}.
