\section{New Graphical Aids for Model Selection}
\label{sec:corbit}
%\addtolength{\textheight}{.5in}%
This section introduces a network analogue of the
autocorrelation and partial autocorrelation function
from classical time series.
Our network autocorrelation function (NACF) reveals 
the autocorrelation structure for different choices of maximum $\hlag$ and $\rstage$ depth, which we denote by $(h, r)$.
We introduce the correlation-orbit (Corbit) plot that displays NACF values and
also serves as an effective graphical aid to perform model selection, namely
selecting good values of $(h, r)$. The NACF depends not only on the
stochastic values of the multivariate time series in a network time
series, but also 
takes account of the network structure and weights. 
\par
%The NACF name we choose coincides with the one in the \texttt{SNA} package aimed at social network analysis; 
%see \cite{sna}. 
%%%%subsection
\subsection{GNAR Network Autocorrelation Function}
Recall the familiar
autocorrelation, $\rho(h)$, for a univariate time series $X_t \in \mathbb{R}$
for lags $h$:
\begin{equation}
    \label{eq:regularacf}
        \rho(h) = \frac{\sum_{t =1}^{T - h} (X_{t + h} - \overline{X}) (X_t - \overline{X}) }{\sum_{t =1}^{T} (X_{t} - \overline{X})^2}.
    \end{equation}
This autocorrelation, $\rho\,(h)$, computes the covariance
between (a mean-corrected) $X_{t}$ and its $h$--lagged
value, $X_{t+h}$ and then normalises by the sample variance in the denominator of~\eqref{eq:regularacf}.
The usual $T^{-1}$ covariance and variance scale normalisations in the numerator and denominator of~\eqref{eq:regularacf} cancel,
see \citet[Section~2.7]{chatfield}

Our network autocorrelation, $\rm{nacf}$, has a similar basic structure to $\rho$ with two main differences.
The first difference is that we examine network autocorrelations not only at time lag $h$ but
also for a given \mbox{$r$-stage} neighbour, so $\rm{nacf} = \rm{nacf}(h, r)$. The second
difference occurs because
the GNAR model~\eqref{eq: vector projection} involves not just past values of $\boldsymbol{X}_t$
but {\em weighted} past values of neighbours and hence those weights need to be
incorporated as the
following definition shows.
\begin{definition}[Network Autocorrelation Function (NACF)]\label{def: nacf}
Suppose that $\boldsymbol{X}_t$ is a GNAR process as defined by~\eqref{eq: vector projection}
that satisfies the
conditions in Theorem~\ref{th: stationary cond}.
Further, let $\bmat{S}_r$ and $\bmat{W}$ be the \mbox{$r$-stage} adjacency and weight
matrices as defined in Definitions~\ref{def: rstage adjacency} and \ref{def: gnar weight matrix},
respectively.
The network autocorrelation function of a $\GNAR$ process $\boldsymbol{X}_t$, with autocovariance bound
\newline
$\lambda := \bigg [ \underset{j = 1, \dots, d}{\max} \big \{\sum_{i = 1}^d [(\bmat{W} \odot \bmat{W})]_{i j} \big \} \bigg ]^{1/2},$ is given by
\begin{equation} \label{eq:nacf definition}
    \nacf(h, r) := \frac{\sum_{t=1}^{T - h} ( \boldsymbol{X}_{t + h} - \boldsymbol{\overline{X}})^{T} \big ( \mathbf{W} \odot \bmat{S}_r + \bmat{I_d} \big )
    ( \boldsymbol{X}_{t} - \boldsymbol{\overline{X}})}
    {\sum_{t=1}^{T}  ( \boldsymbol{X}_{t} - \boldsymbol{\overline{X}})^{T} \big \{ \big (1 + \lambda \big) \bmat{I_d} \big \} ( \boldsymbol{X}_{t} - \boldsymbol{\overline{X}})},
\end{equation}
where $\bmat{I_d}$ is the \mbox{$d$-dimensional} identity matrix.
\end{definition}
The $r$th stage neighbour selection and weighting of those in the nacf is achieved by
the term $\mathbf{W} \odot \bmat{S}_r$ in the numerator of~\eqref{eq:nacf definition}.
The direct comparison of $\boldsymbol{X}_{t+h}$ to its predecessors,
analogous to the numerator in~\eqref{eq:regularacf} is enabled by the $\bmat{I_d}$
term in the numerator of~\eqref{eq:nacf definition}. The introduction of the 
$\mathbf{W} \odot \bmat{S}_r$ term in the numerator of~\eqref{eq:nacf definition} can be thought 
of as a potential `inflation' over the standard numerator and so, to ensure that
$-1 \leq \nacf(h,r) \leq 1$ we need to modify the denominator. We achieve this by
computing the all-encompassing autocovariance bound $\lambda$, which includes all
weights and hence bounds any $\bmat{S}_r$-selected neighbour weights that might appear in
the numerator. Incidentially, the nacf for a trivial GNAR model
with a single node network reduces to the ordinary autocorrelation, i.e. $\rm{nacf}(h, r) = \rho\,(h)$.

For model selection we can employ NACF on a similar basis to that for the regular autocorrelation
for univariate time series.
By comparing the NACF values for different choices of $(h, r)$ we can analyse the lag and $\rstage$ depth at which autocorrelation starts to decay. 
We explore a possible NACF interpretation in Section \ref{sub: nacf interpretation}; see the supplementary material document Section~A for the NACF derivation and some of its properties. 
%%subsection
%%%%%%%%
\subsection{Corbit Plot}
We introduce the Corbit plot by studying realisations coming from a stationary \mbox{global-$\alpha$} $\GNAR(2, [1, 1])$ process  with
parameters $\alpha_1 = 0.23, \alpha_2=0.11, \beta_{1, 1} = 0.21$ and $\beta_{2, 1} = 0.12$.
Assume that we do not know the model order:
we can study the network autocorrelation decay by plotting the observed NACF values via the Corbit plot.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.60]{images/gnar_2_11_sim_nacf.pdf}
    \caption{Corbit NACF plot for $T=200$ samples from a stationary \mbox{global--$\alpha$} $\GNAR(2, [1, 1])$, where the underlying network is the fiveNet network included in the \texttt{GNAR} package; see \cite{gnar_package}. See text for description.}
    \label{fig: corbit sim}
\end{figure}
Figure \ref{fig: corbit sim} introduces the Corbit plot, each point corresponds to 
a specific  \mbox{$h$-lag} and \mbox{$r$--stage} pair, and the colour is set by a colour scale based on the overall NACF values. The first ring depicts $1$--stage neighbours (i.e., nodes with one edge between them), the second ring considers \mbox{2--stage} neighbours (i.e., nodes with a shortest path length equal to two). In short, the ring number counting from the inside corresponds to $\rstage$ depth. 
\par
The numbers on the outside ring indicate the time lag used for computing the NACF. Therefore, the value of each point is $\nacf(h, r)$ where $h$ is the lag denoted on the the last ring and $r$ is the ring that corresponds to \mbox{$r$--stage} adjacency. 
\par
The point on the first inner ring with a number one to the right in Figure \ref{fig: corbit sim} is the value of $\nacf(1, 1)$, the second point on the first inner ring is $\nacf(2, 1)$, and the first point on the second inner ring is $\nacf(1, 2)$; this pattern repeats for subsequent rings and lags.
\par
Currently, in our software, there are two possible ways for assigning size to each point. The default choice is to use the NACF absolute value, i.e., $|\nacf(h, r)|$. This choice highlights the different $\nacf(h, r)$ magnitudes. 
An alternative choice is based on the residual sum of squares constraining $\GNAR$ to a specific lag and stage pair produces; see the supplementary material. 
The value of the largest correlation in Figure~\ref{fig: corbit sim} is
about $0.14$, which might seem a bit
small. This is because the correlations
are weighted
by the (usually inverse distance) weights
matrix $\bmat{W}$, which sum to one for
a given node. So, the pairs that make
up the correlation are weighted by numbers
between zero and one, and hence 
overall
they are generally smaller.
The autocovariance bound in
Definition~\ref{def: nacf}
also reduces the correlations a little.
Examples with larger correlations can
be found in the associated RMarkdown document Section~2.
\par
Finally, the point at the centre has NACF value equal to zero and the smallest size, which highlights the larger NACF values and facilitates comparing the NACF values to $\nacf(h, r) = 0$. We note that these choices are not exclusive and other measures of model fit and/or correlation could be used for assigning point size and colour. 
\par
Corbit plots are produced using the \texttt{viridis} R library, which provides a colour scale that is easily perceived by viewers with common forms of colour blindness (\cite{viridis_lib}) and the \texttt{ggplot} package functionality (\cite{ggplot_lib}).
The {\tt corbit\_plot} function in the {\tt GNAR} package
can also produce other versions of the plot --- one with a line between the highest lag and lag~1
to indicate that these lags are not in actuality close, and also a rectangular version, see
Figures~\ref{fig: corbit + line gap} and~\ref{fig:nacfrect} in
Appendix~\ref{app: furtherplots} below.
\par
The Corbit plot in Figure \ref{fig: corbit sim} shows that the NACF decays on/after lags equal to or larger than four across all stages. Also, the NACF drops after the second stage at the first lag and after the first stage for the second lag. Observe that the further separated nodes are in the network the closer the NACF is to zero. Furthermore, across all $\rstage$ depths we see that the NACF decays to zero as the lag increases. This Corbit plot suggests that autocorrelation for the simulated network time series has larger values for $h \in \{1, 2\}$ and $r \in \{1, 2, 3\}$ and decays sharply as the $\hlag$ and $\rstage$ depth increase. 
\par
This is in accordance with univariate ACF plots which decay as the lag increases, and does reflect the known underlying $\GNAR(2, [1, 1])$ structure.
\subsection{GNAR Partial Autocorrelation Function}
The NACF computes the autocorrelation between $\boldsymbol{X}_t$ and its lagged observations for a specific choice of $\rstage$ neighbourhood regression, however, it does not account for the effects previous and intervening lags and/or $\rstage$ neighbours have when computing $\nacf(h, r)$. This makes diagnosing model order by looking at the NACF values on a Corbit plot challenging since we do not know if autocorrelation has not reduced because of the effects previous lags and/or $\rstage$ depths might have. This difficulty could prevent us from observing the underlying GNAR autocorrelation structure if there is one.
\par
Motivated by techniques from univariate time series analysis, we propose the partial network autocorrelation function (PNACF) as a tool for diagnosing GNAR model selection. The PNACF computes the autocorrelation between $\boldsymbol{X}_t$ and $h$-lagged observations of itself for a specific $\rstage$ neighbourhood regression after the linear effects of previous lags and $\rstage$ neighbours have been removed. The PNACF acts as the partial autocorrelation function does for univariate time series by identifying model order from examining the $(h, r)$ pair after which there is a sharp decline in autocorrelation.
\par
Assume that $\boldsymbol{X}_t$ is a GNAR model given by (\ref{eq: vector wise}) and satisfies the conditions in Theorem \ref{th: stationary cond}, then it is possible to remove the effects from previous lags and $\rstage$ neighbours by focusing on the empirical residuals arising from a $\GNAR(h - 1, [(r - 1)^{(h-1)}])$ fit. These residuals correspond to the best linear prediction restricted to lag $(h - 1)$ and maximum $\rstage$ depth $(r - 1)$, which we denote by $\boldsymbol{X}_{t}^{h - 1, r - 1}:= \sum_{k = 1}^{h - 1} ( \hat{\alpha}_k \boldsymbol{X}_{t - k} + \sum_{s = 1}^{r - 1} \hat{\beta}_{kr} \boldsymbol{Z}_{t - k}^{r - 1}  )$; see the Definition~2 in
the separate supplementary material document. 
\par
One possible extension of the PACF from a univariate setting to the GNAR framework is defining the PNACF as the NACF between the residuals arising from the best linear predictions using $(h - 1)$ lags and $(r - 1)$ $\rstage$ neighbourhood regressions for $\hlag$ and $\rstage$ pairs. Unfortunately, computing the PNACF as mentioned above requires us to have prior knowledge of the autoregressive coefficients in (\ref{eq: vector wise}). We circumvent this by inputting the least-squares estimators as if they were the true parameter values, which is valid given the consistency of the estimator given by (\ref{eq: ls estimator}); see \cite{gnar_paper}. The values we get using the imputed parameters should reflect the underlying GNAR autocorrelation structure if there is one. The new forecast value is $\hat{\boldsymbol{X}}_{t}^{h - 1, r - 1} := \sum_{k = 1}^{h - 1} ( \hat{\alpha}_k \boldsymbol{X}_{t - k} + \sum_{s = 1}^{r - 1} \hat{\beta}_{ks} \boldsymbol{Z}_{t - k}^{s}  )$. Our PNACF is defined below.
\begin{definition} \label{def: pnacf}
For a stationary $\GNAR$ process $\boldsymbol{X}_t$ compute the residuals coming from a $\GNAR(h - 1, [(r - 1)^{(h-1)}])$ fit, the corresponding residual mean $\overline{\boldsymbol{u}}$, and the residuals with $\hlag$s between them by
    $\boldsymbol{\hat{u}}_{t + h} = \boldsymbol{X}_{t + h} - \hat{\boldsymbol{X}}_{t + h}^{h - 1, r - 1}$ and
    $\boldsymbol{\hat{u}}_{t} = \boldsymbol{X}_{t} - \hat{\boldsymbol{X}}_{t}^{h - 1, r - 1}.$
Then, the sample partial network autocorrelation function is
\begin{equation} \label{eq: pnacf}
    \pnacf(h, r) := \frac{\sum_{t=1}^{T - h} ( \boldsymbol{\hat{u}}_{t + h} - \boldsymbol{\overline{u}})^{T} \big ( \bmat{W} \odot \bmat{S}_r + \bmat{I_d} \big )
    ( \boldsymbol{\hat{u}}_{t} - \boldsymbol{\overline{u}})}
    {\sum_{t=1}^{T}  ( \boldsymbol{\hat{u}}_{t} - \boldsymbol{\overline{u}})^{T} \big \{ \big (1 + \lambda \big) \bmat{I_d} \big \} ( \boldsymbol{\hat{u}}_{t} - \boldsymbol{\overline{u}})}, 
\end{equation}
where $\lambda$ is the same as in Definition \ref{def: nacf}.
\end{definition}
The PNACF computes the remaining network autocorrelation between residuals after removing the linear effects of previous lags and stages. Intuitively, it will cut--off to zero at every $(h, r)$ pair where $h > p$ and $r > r^*$, where $r^* := \max \{ s_1, \dots, s_p \}$ is the largest active $\rstage$ depth, since these pairs correspond to the sum of cross-correlations between the IID white noise processes $u_{i, t + h}$ and $u_{j, t}$ given by (\ref{eq: node-wise}).
By computing said network autocorrelations, the PNACF highlights the stage and lag at which the network autocorrelation `cuts off'.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.60]{images/gnar_2_21_sim_pnacf.pdf}
    \caption{Corbit PNACF plot for $T=200$ samples from a stationary \mbox{global-$\alpha$} $\GNAR(2, [2, 1])$ which has the fiveNet network as underlying structure; see \cite{gnar_package}. The maximum lag is equal to 20 and maximum $\rstage$ is equal to 3.}
    \label{fig: pnacf sim}
\end{figure}
Figure \ref{fig: pnacf sim} shows that the PNACF cuts--off at lag three across all stages and
indeed it cuts off for all $\rstage$ depths for $h \geq 3$), and that it `cuts off' at stage one for lag two and at stage two for lag one. This Corbit plot suggests fitting a $\GNAR(2, [2, 1])$ which recovers the known data--generating process in this case. Note that the PNACF cut--off mimics the way in which the PACF decays when looking at univariate time series; see \cite{brockwell_davis}. We will introduce the R--Corbit plot below,
which is a development of the Corbit plot that highlights non--stationarities or indicates the effect of covariates. 