\section{GNAR Properties and  Useful New  Interpretations}
\label{sec: graph model connection}
%\addtolength{\textheight}{.5in}%
GNAR processes generalise graphical models for multivariate time series by introducing higher--order interactions between nodes, and
 resulting in a parsimonious specification that leverages the similarities of the individual node--wise processes.
For our COVID--19 pandemic analysis later, it is reasonable to expect that the SARS--CoV--2 virus
behaves similarly in different locations across England.
This section  introduces
novel  explicit connections between GNAR and graphical models for multivariate time series  and 
further explains
 new  selection and shrinkage properties that GNAR models implicitly possess due to their
incorporation of prior knowledge.
%% subsection
%%%%%%%%
\subsection{Connection to Graphical Models for Time Series}
A goal of GNAR  modelling is to include possible interactions between nodes that are not directly connected in the network. One possible way of generalising the notions in \cite{graph_series} is to extend edge--based interactions by assigning $\rstage$ adjacency based on a cross--spectral hierarchy presented in this subsection. GNAR introduces higher-order interactions into the graphical model by allowing $\rstage$ neighbourhood regression for $r \in \{1, \dots, r^*\}$, where $r^*$ is the largest active $\rstage$ neighbourhood regression. 
\par
These higher--order interactions can be interpreted as weaker dependence between nodes the further separated the nodes are in the graph, in the sense that $\rstage$ neighbourhood regression is less relevant the larger $r$ is and ultimately non-relevant if $r > r^*$. Here, the idea of an edge between nodes is extended to membership in $\rstage$ adjacency sets. Denote by $\mathcal{N}_r$ the set of $\rstage$ neighbours, so if $j \in \mathcal{N}_r (i)$, then, by symmetry and shortest path uniqueness, $(i, j) \in \mathcal{N}_r$, moreover, note that $\mathcal{N}_1$ is the ordinary set of edges. 


An intuition for GNAR is that the cross--correlation between $X_{i, t + h}$ and $X_{j, t}$ at all lags $h$ should be strongest if
$j \in \mathcal{N}_1 (i)$, drop for $j \in \mathcal{N}_r (i)$ where $r \in \{2, \dots, r^*\}$,  and that $X_{i, t+ h}$ and $X_{j, t}$ do not heavily influence each other if $j \notin \mathcal{N} (i)$, where $\mathcal{N} (i) := \cup_{r = 1}^{r^*} \mathcal{N}_r (i)$ is the $r^*$-borough of node $i$ (i.e., collection of neighbourhoods). This motivation relates $\rstage$ neighbourhood regression to the cross--spectrum and the inverse spectral matrix, as the new result shows next.
\begin{theorem} \label{th: pop_sparsity gnar}
    Let $\boldsymbol{X}_t$ be a stationary $\GNAR(p, [s_1, \dots, s_p])$ process with a static network structure $\mathcal{G} = (\mathcal{K}, \mathcal{E})$, full rank spectral matrix $\bmat{f} (\omega)$ and maximum active $\rstage$ depth $r^*$. Then, the inverse spectral matrix $\bmat{S}(\omega)$ and the node-wise distances $\delta(i, j)$ computed on the network $\mathcal{G}$ satisfy
        \par
        $\text{a. }$ There exists a partial correlation graph $\tilde{\mathcal{G}} = (\mathcal{K}, \tilde{\mathcal{E}})$ with the same set of nodes as $\mathcal{G}$ such that
        $$(i, j) \notin \tilde{\mathcal{E}} \textit{ if and only if }  \delta(i, j) \geq 2 r^* + 1.$$
        \par
        $\text{b.}$ There exists a cross-spectral hierarchy $\xi^{(1)} > \cdots > \xi^{(r^*)} >  \xi^{(r^* + 1)} = 0$ and an active $\rstage$ neighbourhood regression which satisfy
         $$\delta(i, j) \in \{2r - 1, 2r\} \textit{ if and only if  } \xi^{(r)} \leq |[\bmat{S}(\omega)]_{i j}| < \xi^{(r - 1)},$$
         for all $\omega \in (-\pi, \pi]$ and for all $r \in \{1, \dots, r^*\}$.
\end{theorem}
See the supplementary material document (Section D) for a proof of Theorem~\ref{th: pop_sparsity gnar}, and \cite{brockwell_davis, ts_forecast_book} for definitions of $\bmat{S}(\omega)$ and $\bmat{f}(\omega)$. Theorem~\ref{th: pop_sparsity gnar} shows that $X_{i,t}$ and $X_{j, t + h}$ are uncorrelated at all lags $h$ given all the other nodes if and only if nodes $i$ and $j$ do not have common active $\rstage$ neighbours. It extends the notion that nodes without edges between them in a graph are uncorrelated given all the other nodes in said graph, by proposing that if the distance between nodes $i$ and $j$ is larger than $2 r^* + 1$, then $X_{i, t}$ and $X_{j, t + h}$ are uncorrelated given all the other nodal time series at all lags, which we express as nodes $i$ and $j$ do not collide in $\rstage$ neighbourhood regression.  
The ideas in this subsection underpin the node relevance statistics in Section~\ref{sec: node influence}. 
\par
Moreover, Theorem~\ref{th: pop_sparsity gnar} exhibits that the higher-order autocorrelation structure a GNAR process induces is equivalent to the process having a hierarchical dependence structure, which can be identified from the inverse spectral matrix (i.e., $\bmat{S}(\omega)$ acts similarly as the concentration matrix does for Gaussian graphical models). This property allows us to interpret the relevance different nodes $j$ have with respect to node $i$ by looking for the $\rstage$ neighbourhood regressions at which $j$ is active (i.e., if $j \in \mathcal{N}_r (i)$, then the smaller $r$ is the more relevant $j$ is). Note that if we restrict GNAR to \mbox{$1$-stage} neighbourhood regression, then higher-order interactions are not considered and the GNAR induced correlation structure is equivalent to a graphical model for multivariate time series. In that case, Theorem~\ref{th: pop_sparsity gnar} recovers results of \cite{graph_series} for partial correlation graphs.
%% subsection
%%%%%%%%
\subsection{Oracle Selection and Shrinkage} \label{sub: subsect oracle selection}
Prior knowledge of the network can be interpreted as having an oracle solution for selecting active nodes in each $\rstage$ neighbourhood regression. A GNAR formulation exploits this by proposing a reparameterization of a constrained VAR process in a parsimonious manner that is related to autoregressive coefficient shrinkage, which reduces estimator variance and improves model interpretability.
 We explain this
novel interpretation next.

\subsubsection*{\textit{For Selection}}
Before performing selection, each node--wise autoregression includes all possible nodes as predictors and the equivalent non-constrained VAR($p$) model is
\begin{equation} \label{eq: var}
     X_{i, t} = \sum_{k = 1}^{p} \big( \phi_{k}^{ii} X_{i , t- k} + \underset{j \neq i}{\sum} \phi_{k}^{ij} X_{j , t- k}
    \big )+ u_{i, t},
\end{equation}
where $\phi_{k}^{ij}$, $\phi_{k}^{ii} \in \mathbb{R}$ are autoregressive coefficients and the $u_{i,t}$ are IID white noise. 


Model~\eqref{eq: var} has $d$ non-zero coefficients for each node-wise autoregression at every lag, hence, the total number of unknown parameters is $p\, d^2$. Following \cite{statLearn}, we say that a node-wise regression is $m$-sparse if only a subset of $m$ predictor nodes have a non-zero autoregressive coefficient. By definition we have that $m \leq d$, however, our focus is on models which highlight a small subset of relevant node predictors (i.e., $m \ll d$). Next, we constrain the VAR model~\eqref{eq: var}
by assuming that only nodes which satisfy $\delta(i, j) \leq r^*$, where the distances are computed on the underlying network, have a non-zero coefficient. Essentially, for all lags $k \in \{1, \dots, p\}$ we impose the constraint
\begin{equation} \label{eq: var constrain}
    \phi_{k}^{ij} \neq 0 \textit{ if and only if } \delta(i, j) \leq r^*,
\end{equation}
by symmetry, we also have that $\phi_{k}^{ij} \neq 0$ if and only if $\phi_{k}^{ji} \neq 0$ at all lags $k$.
\par
Applying  constraint \eqref{eq: var constrain} to the VAR model~\eqref{eq: var} and noting that $j \in \mathcal{N}(i)$ if and only if $\delta(i, j) \leq r^*$ gives
\begin{equation} \label{eq: sparse var}
    X_{i, t} = \sum_{k = 1}^{p} \big( \phi_{k}^{ii} X_{i , t- k} + \underset{j \in \mathcal{N} (i)}{\sum} \phi_{k}^{ij} X_{j , t- k}
    \big )+ u_{i, t},
\end{equation}
where $\mathcal{N}(i) = \cup_{r = 1}^{r^*} \mathcal{N}_{r} (i)$ and the $u_{i,t}$ are IID white noise. Above, each $i$th nodal time series has at most $m^{(i)} := |\mathcal{N}(i)| + 1 \leq d$ non-zero autoregressive $\phi_{k}^{ij}$, $\phi_{k}^{ii}$ coefficients at each lag $k$. Thus, under  constraint~\eqref{eq: var constrain}, each node--wise regression is $m^{(i)}$-sparse at each lag and has at most $p\, m^{(i)}$ unknown parameters. 
\par
We set $m := \max \{ m^{(i)} \}$ and see that all node--wise autoregressions are $m$--sparse and that the model has at most $p\, m^2$ unknown parameters. Furthermore, constraint~\eqref{eq: var constrain}
performs variable selection for all node-wise autoregressions, and reflects
our assumption of closer nodes being more relevant to each other. 


Variable selection can be performed by multiplying the vector time series by the sum of $\rstage$ adjacency matrices. Let $\bmat{S} = \sum_{r = 1}^{r^*} \bmat{S}_r$, and note that $[\bmat{S}]_{ij} \neq 0$ if and only if $\delta(i, j) \leq r^*$, and that $[\bmat{S}]_{ii} = 0$ for all nodes $i$, hence, model~\eqref{eq: sparse var} is equivalent to
\begin{equation} \label{eq: selection equivalence}
    X_{i, t} = \sum_{k = 1}^{p} \big( \phi_{k}^{ii} X_{i , t- k} + \underset{j = 1}{\sum^d} \phi_{k}^{ij} [\bmat{S}]_{ij} X_{j , t- k}
    \big )+ u_{i, t},
\end{equation}
which has the same active nodes (i.e., $\phi_k^{ij} [\bmat{S}]_{ij} \neq 0 \iff \delta(i, j) \leq r^*$) as the GNAR vector-wise representation~\eqref{eq: vector wise}. This establishes variable selection equivalence between a GNAR formulation and a constrained VAR, where the number of selected nodes for each node-wise autoregression satisfies $m^{(i)} = \sum_{j = 1}^{d} [\bmat{S}]_{ij}$.
\par
Moreover, fixing $r^*$ in a GNAR formulation is equivalent to imposing a $\ell_0$-ball constraint on the autoregressive coefficients, which can be efficiently approximated by an $\ell_1$-norm constraint; 
see Chapter~7.2.1 in~\cite{wainwright_2019}. Further, decreasing $r^*$ reduces the number of nodes included in node-wise autoregressions, which results in a sparser model by noting that at all lags $k$
\begin{equation*}
    0 < \normx{ \bmat{\Phi}_k \odot \bmat{S}_1 }_1 \leq \nu^{(1)} \leq \cdots \leq  \normx{ \bmat{\Phi}_k \odot \bmat{S} }_1 \leq \nu <  \normx{ \bmat{\Phi}_k }_1,
\end{equation*}
hence, as $r^*$ decreases the number of active nodes decreases too, which results in a smaller $\nu > 0$ (i.e., a tighter constraint) and a sparser model.
\par
Based on the above, we interpret $\rstage$ neighbours as the relevant predictors each node has, and we can think of maximum $\rstage$ depth (i.e., $r^*$) as a hyperparameter that controls variable selection.
\subsubsection*{\textit{For Shrinkage}}
GNAR accounts for the possibility that nodes in the same $\mathcal{N}_r$ are likely to be correlated, which will affect estimator and predictive performance (under a valid GNAR model); 
see Chapter~3 in~\cite{statLearn}. Motivated by this, a GNAR formulation performs shrinkage on the selected variables at the population level by reparametrizing the constrained VAR model as follows. 
\par
Using the notation and concepts from a GNAR formulation the node--wise autoregressions in a sparse VAR($p$) are given by (\ref{eq: sparse var}). We assume that relevance weights $\sigma_{ij} > 0$ between nodes are available, which quantify the similarity between nodes $i$ and $j$, and can be computed as $\sigma_{ij} := \{h(i,j)\}^{-1}$, where $h(i,j) > 0$ measures a notion of distance between nodes, which does not have to be equal to the distance in the graph (but could be),
for instance, the distance can be a function of some exogenous variable that measures `closeness' between nodes in a non--linear manner. Essentially, the $\sigma_{ij}$ are a form of hard prior knowledge.
\par
The notion of closeness and `similarity' between nodes that the $\sigma_{ij}$ provide can be incorporated into the GNAR framework by reparameterizing the autoregressive coefficients in (\ref{eq: sparse var}) as $\beta_{kr}\, \sigma_{ij} :=  \phi_{k}^{ij}\, [\bmat{S}]_{ij}$, where prior knowledge of $\sigma_{ij}$ makes the $\beta_{kr}$ identifiable. Also, it means that if $\delta(i, j) > r^*$, then $\beta_{kr}\, \sigma_{ij} = 0$. This parameterization assumes that there is a global effect 
$\beta_{kr}$ common across $\rstage$ neighbours, which results in a parsimonious representation, i.e.,\ a GNAR
model is valid.
\par
To satisfy the assumptions of Theorem~\ref{th: stationary cond} (i.e., weights less than or equal to one), we normalise the weights by computing updates as $w_{ij} := \sigma_{ij} ( \sum_{l \in \mathcal{N}_r (i)} \sigma_{il} )^{-1}$, hence, 
$\sum_{j \in \mathcal{N}_r (i)} w_{ij} = 1$ for all nodes at every active $\rstage$ neighbourhood regression. Finally, to transform (\ref{eq: sparse var}) into a model equivalent to a GNAR model shrink the active $\phi_{k}^{ij}$ coefficients as follows
\begin{equation} \label {eq: shrink gnar}
\gamma_{k}^{i j} :=  \phi_{k}^{i j}\, \{ 1 + \sigma_{i j}^{-1}\, v_{r}(i, j) \}^{-1},
\end{equation}
where $ v_r (i, j) := \sum_{l \in \mathcal{N}_r (i) - \{j \} } \sigma_{il}$ - i.e., weight normalisation can be thought of as a constraint on the $\ell_2$-norm of the active coefficients. We further note an interesting connection between the above and parameter updates in linear ridge regression.

\begin{proposition} \label{prop: shrinkgae relation}
Let $\boldsymbol{y} \in \mathbb{R}^n$ be a vector of responses, $\bmat{A} \in \mathbb{R}^{n \times m}$ the design matrix and $\boldsymbol{\phi} \in \mathbb{R}^m$ the vector of unknown linear coefficients. Next define $\hat{\boldsymbol{\phi}} := \argmin \{ || \boldsymbol{y} - \bmat{A}\, \boldsymbol{\phi} ||_2^2 \},$ and $\hat{\boldsymbol{\gamma}} := \argmin \{ || \boldsymbol{y} - \bmat{A}\, \boldsymbol{\gamma} ||_2^2 \}$ such that $|| \boldsymbol{\gamma}||_2 \leq \lambda$, where $\boldsymbol{\gamma} \in \mathbb{R}^m$ and $\lambda > 0$. If the design matrix admits the singular value decomposition $\bmat{A} = \bmat{U}\, \bmat{\Sigma} \,\bmat{V}^T$, then we have that
\begin{equation} \label{eq: shrink sparse ridge}
    \hat{\gamma}_j = \hat{\phi}_j\, (1 + v\, \sigma_j^{-2} )^{-1},
\end{equation}
where $\sigma_{j}^{2}$ are the diagonal entries in $\bmat{\Sigma}^2$ and $v > 0$ is the Lagrange multiplier linked to the constraint $|| \boldsymbol{\gamma}||_2 \leq \lambda$.
\end{proposition}

In view of Proposition \ref{prop: shrinkgae relation}, and by comparing (\ref{eq: shrink gnar}) with (\ref{eq: shrink sparse ridge}), we interpret each $i$th node-wise $\rstage$ neighbourhood regression as a sparse linear regression with shrunken coefficients, which satisfy the weight normalisation constraint $\gamma_{k}^{ij} =  \phi_{k}^{i j}\, \{ 1 + \sigma_{i j}^{-1} \,v_{r}(i, j) \}^{-1} \,[\bmat{S}]_{ij}$ for all nodal pairs at all lags. Moreover, by (\ref{eq: shrink gnar}) and the above, 
the usual estimated GNAR coefficients can thus be re-expressed by a particular transformed and shrunk version of VAR coefficients
that we have just explicitly represented by $\gamma_{k}^{ij} = \beta_{kr}\, w_{ij}$.
\par
As above, the $\beta_{kr}\, w_{ij}$ are the node--wise coefficients in the GNAR parameterization given by~\eqref{eq: node-wise}, which is a constrained reparametrization of the sparse VAR($p$) process~\eqref{eq: sparse var} with shrunken coefficients $\gamma_{k}^{ij}$. This reparametrization results in a parsimonious model with $p + \sum_{k = 1}^{p} s_k \ll p\, m^2 \ll p\, d^2$ unknown parameters. Hence, GNAR processes can be interpreted as constrained VAR models which exploit knowledge of the underlying network structure to perform variable selection and shrinkage at the population level. Our GNAR framework produces parsimonious models that are highly interpretable.
%%New subsection
%%%%%%%%%%%

\subsection{NACF Interpretation}
\label{sub: nacf interpretation}
The connection between GNAR processes, graphical models for multivariate time series and sparse VAR models allows us to interpret the NACF as a measure of the constrained (i.e., 
network--induced) correlation structure in the process. This effect is accounted for by the autocovariance bound $\lambda$ given by Definition \ref{def: nacf}.
\par
The sparseness and shrinkage properties GNAR models have permit us to bound the dot product
$\langle \boldsymbol{x}, ( \mathbf{W} \odot \bmat{S}_r + \bmat{I_d}) \tildesym{x} \rangle $ for all $\rstage$ neighbourhood regressions as follows
\begin{equation} \label{eq: nacf inequality}
    \big | \big \langle \boldsymbol{x}, ( \mathbf{W} \odot \bmat{S}_r + \bmat{I_d}) \tildesym{x}  \big \rangle \big | 
    \leq ( 1 + \lambda) \normx{\boldsymbol{x}} \normx{\tildesym{x}},
\end{equation}
which illustrates the inner workings of the NACF; see the supplementary material document (Proposition~3, in Section~B) for a proof of (\ref{eq: nacf inequality}). Essentially, $\lambda$ is a constant depending on the network structure that globally takes into account and corrects for the selection and shrinkage properties the weights $w_{ij}$ and $\rstage$ adjacency have on the network time series.
\par
The underlying GNAR structure points vectors in the direction of $(\mathbf{W} \odot \bmat{S}_r + \bmat{I_d})$ for which the magnitudes of the cross-covariances (i.e., the centered dot products) between $\boldsymbol{X}_{t + h}$ and $( \mathbf{W} \odot \bmat{S}_r + \bmat{I_d})\boldsymbol{X}_{t}$ are bounded by the magnitude of the autocovariance between independent components being projected onto the hypersphere $\{ \boldsymbol{x}: \normx{\boldsymbol{x} - \overline{\boldsymbol{x}}}_2 \leq (1 + \lambda) \}$. 
\par
By Theorem~\ref{th: pop_sparsity gnar} and Section \ref{sub: subsect oracle selection}, decreasing $r^*$ makes the model sparser and $\lambda$ smaller, and
the Corbit plot exploits this to highlight the lag and stage pairs at which the NACF cuts-off. Also, since the operator norm of $(\mathbf{W} \odot \bmat{S}_r + \bmat{I_d}) $ is 
upper--bounded by $(1 + \lambda)$, the NACF highlights the $(h,r)$ pairs at which the projected vectors point in the direction of the eigenvector associated to the largest eigenvalue of $(\mathbf{W} \odot \bmat{S}_r + \bmat{I_d})$, hence, the NACF is close to $\,p\,m 1$ when the network weights and $\rstage$ adjacency are such that the magnitude of projections compared in (\ref{eq: nacf inequality}) are approximately the same. In that extreme case, future observations can be seen as eigenvectors of the network structure matrix $(\mathbf{W} \odot \bmat{S}_r + \bmat{I_d})$.
