\section{GNAR simulations and the Corbit plot}
\label{sec: gnar sims}
This section illustrates the usefulness of the network autocorrrelation function (NACF) and the partial network autocorrelation function (PNACF). All of the models we simulate in this section are simulated using the {\bf fiveNet} network included in the \texttt{GNAR} package; see \cite{gnar_package}. Figure \ref{fig:graph_image} shows the graph associated with the {\bf fiveNet} network. 
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.25]{images/fiveNet.pdf}
    \caption{fiveNet graph, this network will be used for all the simulations in Section \ref{sec: gnar sims}. We use this network to simulate three different GNAR processes with the same graph structure.}
    \label{fig:graph_image}
\end{figure}
We assume that the network only has one covariate, thus, we drop the covariate from weight and coefficient specification. The weights are given by normalising to one the sum of inverse distances between \mbox{$r$-stage} neighbour nodes. Also, in what follows we will not study time-varying networks, thus, we drop the $t$ subscript and write $\bmat{W} = \bmat{W}_{t}$.
\par
Thus, for this GNAR specification, we have that $\delta (i, j) = \tilde{w}_{ij}^{-1}$, so for $j \in \mathcal{N}_r(i)$ the weight is given by $w_{ij} = \tilde{w}_{ij} \big (
\sum_{q \in \mathcal{N}_r(i)} \tilde{w}_{ij} \big )^{-1}$. Therefore, the corresponding matrix for node-distance based weights with only one covariate is
\begin{equation*}
    \mathbf{W} = \begin{pmatrix} 0 & 0.5 & 0.5 & 0.5 & 0.5 \\ 1 & 0 & 0.5 & 0.5 & 1 \\ 1 & 0.5 & 0 & 0.5 & 1 \\ 0.33 & 0.33 & 0.33 & 0 & 1 \\ 1 & 0.5 & 0.5 & 1 & 0 
    \end{pmatrix} 
    ,
\end{equation*}
and the $r$-stage adjacency matrices are:
\begin{align*}
    \mathbf{S_1} &= \begin{pmatrix} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ 0 & 1 & 0 & 1 & 0 \\ 1 & 1 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \end{pmatrix} 
    , \\
    \mathbf{S_2} &= \begin{pmatrix} 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix},  \\
    \mathbf{S_3} &= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 \end{pmatrix} .
\end{align*}
Note that the row sums of $\big ( \bmat{W} \odot \bmat{S}_r \big )$ are equal to one for all $r \in \{1, 2, 3 \}$, as required. We compute these matrices with the following code:
\begin{verbatim}
    W <- weights_matrix(fiveNet, 3)

    stages_tensor <- get_k_stages_adjacency_tensor(
                    as.matrix(GNARtoigraph(fiveNet)), 3).
\end{verbatim}
In what follows, we present three different GNAR simulations. These simulations are chosen to highlight the differences we observe in their respective Corbit plots. We simulate three five nodal network time series as proxies for studying different GNAR autocorrelation structures. The three simulations are shown in Figure \ref{fig:sim_one}. 
\begin{figure}[H]
    \centering
    \includegraphics[height=12cm]{images/five_sims_joint.pdf}
    \caption{Simulation with 150 time-steps combining three stationary GNAR data-generating process. The first 50 time-steps come from a GNAR(2, [2, 2]) with parameters given by \eqref{eq: simulation one}, the subsequent 50 time-steps come from a GNAR(1, [1]) with parameters given by \eqref{eq: simulation two} and the last 50 time-steps come from a GNAR(1, [2]) with parameters given by \eqref{eq: simulaiton three}. All three simulations use the {\bf fiveNet} network as underlying stucture, which is included in the \texttt{GNAR} package.}
    \label{fig:sim_one} 
\end{figure}
\subsection{Simulation one}
The first simulation comes from a local-$\alpha$ GNAR(2, [2, 2]), the code for the simulation is

\begin{verbatim}
sim1 <-  GNARsim(n = 50, net=fiveNet,
            alphaParams = list(c(0.1, 0.12, 0.16, 0.075, 0.21),
            c(0.12, 0.14, 0.15, 0.6, 0.22)), 
            betaParams = list(c(0.1, 0.16), c(0.11, 0.14)))
\end{verbatim}

The data-generating process has the following parameter values
\begin{align} \label{eq: simulation one}
    \boldsymbol{\alpha}_{\cdot 1} &= (0.1, 0.12, 0.16, 0.075, 0.21) , \nonumber \\
    \boldsymbol{\alpha}_{\cdot 2} &= (0.12, 0.14, 0.15, 0.6, 0.22) , \nonumber \\
    \mathbf{A}_1 &= \diag (\boldsymbol{\alpha}_{\cdot 1}) , \\
    \mathbf{A}_2 &= \diag (\boldsymbol{\alpha}_{\cdot 2}) , \nonumber \\
    \mathbf{B}_1 &= 0.1 \big ( \mathbf{W} \odot \mathbf{S}_1 \big ) + 0.16 \big ( \mathbf{W} \odot \mathbf{S}_2 \big ) , \nonumber \\
    \mathbf{B}_2 &= 0.11 \big ( \mathbf{W} \odot \mathbf{S}_1 \big ) + 0.14 \big ( \mathbf{W} \odot \mathbf{S}_2 \big ) \nonumber . 
\end{align}
%\newpage
Thus, the simulations are fifty realisations of the stationary VAR process
\begin{equation*}
    \boldsymbol{X}_t = \left ( \bmat{A}_1 + \bmat{B}_1 \right ) \boldsymbol{X}_{t - 1} +  \left ( \bmat{A}_2 + \bmat{B}_2 \right ) \boldsymbol{X}_{t - 2} + \boldsymbol{u}_t,
\end{equation*}
where $\boldsymbol{u}_t \sim \mathcal{N} (0, \sigma^2 \bmat{I_5})$ are independent for different time steps $t$ (i.e., a homoscedastic white noise process). We study the network autocorrelation of the simulations by analysing the Corbit plots below. 
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{images/nacf_corbit_sim1.pdf}
    \caption{NACF Corbit plot for 50 realisations simulated from the model given by \eqref{eq: simulation one}. The first inner ring corresponds to 1-stage neighbours, the second to 2-stage neighbours and so on. The numbers on the outermost ring denote lags. Each point is $\nacf(h, r)$ and size is assigned as $|\nacf(h, r)|$. The point at the centre has the smallest size and value $\nacf(\cdot, \cdot) = 0$ for aiding comparison.
    }
    \label{fig: nacf sim1}
\end{figure}
An alternative option for assigning size is to use the coefficient of determination computed from the linear model
\begin{equation*}
    E ( X_{i, t + h} |  \boldsymbol{X}_{t} ) =  \alpha_h X_{i, t} + \beta_{hr} \left \langle [(\bmat{W} \odot \bmat{S}_r)]_{i \cdot}, \boldsymbol{X}_{t} \right \rangle . 
\end{equation*}
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{images/pnacf_corbit_sim1.pdf}
    \caption{PNACF Corbit plot for 50 realisations simulated from the model given by \eqref{eq: simulation one}. The first inner ring corresponds to 1-stage neighbours, the second to 2-stage neighbours and so on. The numbers on the outermost ring denote lags. Each point is $\pnacf(h, r)$ and size is assigned as $|\pnacf(h, r)|$. The point at the centre has the smallest size and value $\pnacf(\cdot, \cdot) = 0$ for aiding comparison.}
    \label{fig: pnacf sim1}
\end{figure}
We remark that the Corbit plot in Figure \ref{fig: pnacf sim1} reflects the known model order $\GNAR(2, [2, 2])$, i.e., partial autocorrelation cuts-off after the second lag, and after the second stage at both the first and second lag. The code for producing these Corbit plots is
\begin{verbatim}
    corbit_plot(vts = sim1, net = fiveNet, max_lag = 10, 
        max_stage = 3)
    corbit_plot(vts = sim1, net = fiveNet, max_lag = 10, 
        max_stage = 3, partial = "yes")
\end{verbatim}

\subsection{Simulations two and three}
Now, we simulate two contrasting GNAR processes to highlight the different correlation structure between the simulations using Corbit plots. We will contrast two models, one is predominantly autoregressive (i.e., the autoregressive parameter $\alpha$ is larger than the neighbourhood regression parameter $\beta$) and the other one has the opposite property (i.e., its neighbourhood regression parameter $\beta$ is larger than the autoregressive parameter $\alpha$). The first model is a GNAR(1, [1]) with the following specification
\begin{align} \label{eq: simulation two}
    \alpha_{1} &= 0.83, \nonumber \\
    \mathbf{A}_1 &= \alpha_1 \bmat{I}, \nonumber \\
    \mathbf{B}_1 &= 0.12 \big ( \mathbf{W} \odot \mathbf{S}_1 \big ) \nonumber, \\
    \boldsymbol{X}_t &= \left ( \bmat{A}_1 + \bmat{B}_1 \right ) \boldsymbol{X}_{t - 1} + \boldsymbol{u}_t.
\end{align}
The second model is a GNAR(1, [2]) with the following specification
\begin{align} \label{eq: simulaiton three}
    \alpha_{1} &= 0.12 , \nonumber \\
    \mathbf{A}_1 &=  \alpha_1 \bmat{I} , \nonumber \\
    \mathbf{B}_1 &= 0.35 \big ( \mathbf{W} \odot \mathbf{S}_1 \big ) + 0.47 \big ( \mathbf{W} \odot \mathbf{S}_2 \big ) \nonumber, \\
    \boldsymbol{X}_t &= \left ( \bmat{A}_1 + \bmat{B}_1 \right ) \boldsymbol{X}_{t - 1} + \boldsymbol{u}_t.
\end{align}
In the first model \eqref{eq: simulation two} we would expect network autocorrelation to decay mostly with respect to lag, whereas, in  the second model \eqref{eq: simulaiton three} we would expect decay to take longer and highlight stage depth dependency. To analyse these propositions we produce the corresponding Corbit plots for the NACF and the PNACF using realisations from these two simulations. The code for producing these simulations is:
\begin{verbatim}
    sim2 <- GNARsim(n = 50, alphaParams = list(rep(0.83, 5)), 
        betaParams = list(c(0.12)))

    corbit_plot(vts = sim2, net = fiveNet, max_lag = 10, 
        max_stage = 3)
    corbit_plot(vts = sim2, net = fiveNet, max_lag = 10,
        max_stage = 3, partial = "yes")
  ###########################################################

    sim3 <- GNARsim(n = 50, alphaParams = list(rep(0.12, 5)),
        betaParams = list(c(0.35, 0.47)))

    corbit_plot(vts = sim3, net = fiveNet, max_lag = 10, 
        max_stage = 3)
    corbit_plot(vts = sim3, net = fiveNet, max_lag = 10, 
        max_stage = 3, partial = "yes")

\end{verbatim}
\begin{figure}[H]
        \centering
        \begin{subfigure}{0.48\textwidth}
            \includegraphics[width=\textwidth]{images/nacf_corbit_sim2.pdf}
            \caption{}
            \label{fig: nacf sim2}
        \end{subfigure}
        \hfill
        \begin{subfigure}{0.48\textwidth}
            \includegraphics[width=\textwidth]{images/nacf_corbit_sim3.pdf}
            \caption{}
            \label{fig: nacf sim3}
        \end{subfigure}
        \hfill
        \begin{subfigure}{0.48\textwidth}
            \includegraphics[width=\textwidth]{images/pnacf_corbit_sim2.pdf}
            \caption{}
            \label{fig: pnacf sim2}
        \end{subfigure}
        \hfill
        \begin{subfigure}{0.48\textwidth}
            \includegraphics[width=\textwidth]{images/pnacf_corbit_sim3.pdf}
            \caption{}
            \label{fig: pnacf sim3}
        \end{subfigure}
        \caption{
            Figure \ref{fig: nacf sim2} NACF Corbit plot for the simulation given by \eqref{eq: simulation two}, \ref{fig: nacf sim3} NACF Corbit plot for the simulation given by \eqref{eq: simulaiton three}. Figure \ref{fig: pnacf sim2} PNACF Corbit plot for the simulation given by \eqref{eq: simulation two}, \ref{fig: pnacf sim3} PNACF Corbit plot for the simulation given by \eqref{eq: simulaiton three}
        }
        \label{fig: pnacf level plots}
\end{figure}
The Corbit plots in Figure \ref{fig: pnacf level plots} highlight the different correlation decay between the simulations given by \eqref{eq: simulation two} and \eqref{eq: simulaiton three}.
\subsection{R-Corbit plots}
We propose the R-Corbit plot as a exploratory tool for studying whether the network autocorrelation structure remains similar throughout different time slices when analysing network time series data, our aim is to diagnose whether or not there is evidence of stationarity in the observed time series. 
\par
Furthermore, we might also want to study if the network autocorrelation structure changes with respect to different covariates, this is particularly useful when analysing nodal time series for which we know there might be exogenous variables that might aid us in understanding the similarities and differences between clusters on the network. For example, we can use the population size (urban, rural and hub town) when studying economic and epidemiological data.
\par
For the purposes of illustration and simplicity we use R-Corbit plots for analysing the different network autocorrelation structure present in the models given by \eqref{eq: simulation one}, \eqref{eq: simulation two} and \eqref{eq: simulaiton three}.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.35]{images/nacf_R-Corbit.pdf}
    \caption{R-Corbit plot for the simulations coming from the models given by \eqref{eq: simulation one}, \eqref{eq: simulation two} and \eqref{eq: simulaiton three}. This R-Corbit plot is for NACF values with a maximum lag of twenty and a maximum $r$-stage depth of three. This plot compares the three simulations as proxies for different time slices, we split the realisations every 50 time-steps. See Figure \ref{fig:sim_one} for a plot of the realisations.}
    \label{fig:R-Corbit_nacf}
\end{figure}

Figure \ref{fig:R-Corbit_nacf} shows a R-Corbit plot for NACF values. It allows us to compare the network autocorrelation structure between the three GNAR realisations coming from \eqref{eq: simulation one}, \eqref{eq: simulation two} and \eqref{eq: simulaiton three}.
As for the Corbit plot, the numbers next to the last inner ring correspond to the lag, and each inner ring corresponds to stage adjacency depth, furthermore, we compare the NACF between the three simulations at each ($h$, $r$) pair by plotting each value around the mean value for the time slices. The R-Corbit plot includes a legend to its right where we can see the time slices and/or covariates we use for that specific R-Corbit plot. Figure \ref{fig:R-Corbit_nacf} suggests that the autocorrelation structure is not equivalent among these time slices of the network time series realisation. Further, it also shows that autocorrelation decays for all three models with respect to lag and $r$-stage depth. We can also produce the R-Corbit plot for the PNACF values. This allows us to compare the partial network autocorrelation structure between different time slices and/or covariates. In this work we are only comparing time slices.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.35]{images/pnacf_R-Corbit.pdf}
    \caption{R-Corbit plot for the simulations coming from the models with parameter values given by \eqref{eq: simulation one}, \eqref{eq: simulation two} and \eqref{eq: simulaiton three}. This R-Corbit plot is for PNACF values with a maximum lag of twenty and a maximum stage depth of three. On this plot we compare the three simulations as proxies for different time slices, we split the realisations every 50 time-steps. See Figure \ref{fig:sim_one} for a plot of the realisations.}
    \label{fig:R-Corbit_pnacf}
\end{figure}

Figure \ref{fig:R-Corbit_pnacf} shows that the R-Corbit plot highlights the strong autoregressive model given by \eqref{eq: simulation two}, this point has the largest PNACF value among the three models and sharpest decay with respect to lag and stage depth. Moreover, we see that the third point at each ($h$, $r$) pair highlights the strong neighbourhood regression present in the model given by \eqref{eq: simulaiton three}, these points show that decay is sharp after the second stage, and not as sharp just for lag one. Finally, the R-Corbit plot exhibits that the model given  by \eqref{eq: simulation one} is more complex as the PNACF decay is not as sharp as in the other particular cases. The code for producing these R-Corbit plots is:
\begin{verbatim}
    sim_frames <- list(sim1, sim2, sim3)

    R_corbit_plot(vts_frames = sim_frames, 
        network_list = list(fiveNet), max_lag = 10,
        max_stage = 3, weight_matrices = list(W), 
        same_net = "yes", 
        frame_names = c("sim1", 'sim2', 'sim3'))
    
    R_corbit_plot(vts_frames = sim_frames, 
        network_list = list(fiveNet), max_lag = 10, 
        max_stage = 3, weight_matrices = list(W), 
        same_net = "yes", 
        frame_names = c("sim1", 'sim2', 'sim3'), 
        partial = "yes")
\end{verbatim}