\documentclass[12pt]{article}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{enumerate}
\usepackage{natbib}
\usepackage{url} % not crucial - just used below for the URL 
\usepackage{amssymb}
\usepackage{amsthm}
\DeclareMathOperator*{\argmin}{arg\, min}
\usepackage{wrapfig}
\usepackage{mathtools}
\usepackage{float}
% \usepackage{setspace}\doublespacing
\usepackage[hidelinks]{hyperref}
\usepackage{mathtools}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{scalerel,stackengine}
\usepackage{relsize}
\usepackage{mathrsfs}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{xcolor}
\usepackage{changebar}
\pagenumbering{arabic}



\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{proposition}{Proposition}
\newtheorem{remark}{Remark}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newcommand{\size}[2]{{\fontsize{#1}{0}\selectfont#2}}
\newcommand\normx[1]{\left\Vert#1\right\Vert}
\newcommand\bmat[1]{\mathbf{#1}}
\newcommand\myprod[2]{#1^{T} #2}
\newcommand\tildesym[1]{\boldsymbol{\tilde{#1}}}
\newcommand{\nacf}[0]{\operatorname{nacf}}
\newcommand{\colsums}[0]{\operatorname{colsums}}
\newcommand{\pnacf}[0]{\operatorname{pnacf}}
\newcommand{\GNAR}[0]{\operatorname{GNAR}}
\newcommand{\diag}[0]{\operatorname{diag}}
\newcommand{\bigchi}[0]{\mathlarger{\mathlarger{ \mathlarger{ \chi}}}}
\newcommand{\rstage}[0]{\mbox{$r$-stage}}
\newcommand{\hlag}[0]{\mbox{$h$-lag}}
\newcommand{\cov}[0]{\operatorname{cov}}
\newcommand{\sft}[0]{\operatorname{sft}}
\newcommand{\rscc}[0]{\operatorname{rscc}}
\newcommand{\local}[0]{\operatorname{local}}
\newcommand{\globindex}[0]{\operatorname{globindex}}


\renewcommand{\qedsymbol}{$\blacksquare$}


%\pdfminorversion=4
% NOTE: To produce blinded version, replace "0" with "1" below.
\newcommand{\blind}{1}

% DON'T change margins - should be 1 inch all around.
\addtolength{\oddsidemargin}{-.5in}%
\addtolength{\evensidemargin}{-1in}%
\addtolength{\textwidth}{1in}%
\addtolength{\textheight}{1.7in}%
\addtolength{\topmargin}{-1in}%




\begin{document}
\nocite{*}



\def\spacingset#1{\renewcommand{\baselinestretch}%
{#1}\small\normalsize} \spacingset{1}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\if1\blind
{
  \title{\bf New tools for network time series with an application
    to COVID-19  hospitalisations}
  \author{G. P. Nason\footnotemark[1] \footnotemark[2],
    D. Salnikov\thanks{Dept.\ Mathematics, Huxley Building, Imperial College, 180 Queen's Gate, South Kensington, London, SW7 2AZ, UK.} \footnotemark[2]\\
    Imperial College London\\
    and \\
    M. Cortina--Borja\thanks{Population, Policy and Practice Research and Teaching Department, Great Ormond Street Institute of Child Health, 30 Guilford Street, London, WC1N 1EH, UK.} \footnotemark[1] \\
    Great Ormond Street Institute of Child Health\\
    University College London}
  \maketitle
} \fi

\if0\blind
{
  \bigskip
  \bigskip
  \bigskip
  \begin{center}
    {\LARGE\bf New tools for network time series with an application
    to COVID--19  hospitalisations}
\end{center}
  \medskip
} \fi


%\bigskip
\begin{abstract}
Network time series models are increasingly important across many areas, involving known or inferred underlying network structure, which can be exploited to make sense of high--dimensional dynamic phenomena. We introduce two new association measures: the network and partial network autocorrelation functions and define Corbit (correlation--orbit) visualisation plots. Corbit plots permit interpretation of underlying correlation structures and, crucially, aid model selection more rapidly than general tools such as information criteria. We introduce interpretations of generalised network autoregressive (GNAR) processes as generalised graphical models.
We shine new light on how incorporating prior information is related to variable selection and shrinkage in the GNAR context. We illustrate the usefulness of GNAR models, network autocorrelations and Corbit plots for
 a novel network time series modelling of COVID--19 mechanical ventilation bed occupancies at 140 NHS Trusts. We also introduce the R--Corbit plot that shows correlations over different time periods or with respect to external covariates and plots that quantify the relevance of individual nodes. Our analysis provides insight on the COVID--19 series’ underlying dynamics, highlights two groups of geographically co--located `relevant' NHS Trusts, and demonstrates excellent predictive performance.
\end{abstract}

\noindent%
{\it Keywords:}  multivariate time series, network autocorrelation function, partial network autocorrelation function, 
Corbit plot, R--Corbit  plot.

% Comment   

%\spacingset{1.9} % DON'T change the spacing!
\section{Introduction}
\label{sec:intro}
Network time series in many areas benefit from
the development of statistical methods that enable interpretation and inference of relationships present in dynamic phenomena. Modelling such network time series necessitates studying a constant flux of complex data characterised by a temporal component among large numbers of interacting variables.
The interacting variables are often associated with a network structure, for example, networks in
neuroscience, biology, medicine and business to name a few. In the absence of a network,
techniques
are available to learn possible network structure, which enables efficient modelling, analysis and forecasting
of a large number of interactions; see, e.g., \cite{lauritzen, graph_lasso, Songsiri200989}. A complete overview of network science and multivariate time series is beyond the scope of this 
work; see \cite{network_time_series_overview} for a review.
For a non-exhaustive overview of network methods in statistics {\em and}\ time series see, e.g.,
\cite{graph_series, lauritzen, luth, brockwell_davis, Songsiri200989, sna_book, ts_forecast_book} and
\cite{DALLAKYAN2022107557}.
\par
Recently, the generalized network autoregressive (GNAR) model has been developed
\citep{gnar_org, Zhu17, gnar_paper}, which provides a parsimonious interpretable model that has been often shown to
have both simpler interpretability {\em and} superior forecasting performance in a number of scenarios.
Model extensions in this rapidly
developing area include, for example, \cite{Zhu19} for quantiles, \cite{Zhou20} for Network GARCH models,
\cite{NasonWei22} to admit time-changing exogeneous variables,
\cite{ArmillottaFokianos21, ArmillottaFokianos24} for Poisson/count data, \cite{Mantziou23} for GNAR processes on
the edges of networks,
and \cite{nar_local_local} develop a model with local per--node network parameters along with a test for equality and subsequent asymptotic theory. 
Such models have proven useful for many (network) time series where the characteristics of the series are similar from variable (node) to variable (node)
in a network, although GNAR processes' utility is not limited to this situation.

Crucial elements of any statistical modelling exercise are model elicitation and specification.
For many time series models, and for GNAR in particular, this involves
choosing
quantities such as the order $p, q$ of any autoregressive and moving average terms, respectively, and the differencing
parameter, $d$.
GNAR models also involve $p$,
and, in addition, the \mbox{$p\times 1$} vector, $s$:
the number of 
`stage-neighbours' per lag. Until now, the
main tools for GNAR model
order choice have been Akaike's and the Bayesian information criteria (AIC
and BIC),
GNAR versions
of which appear in the \texttt{GNAR} CRAN package developed by~\cite{gnar_package}
for \texttt{R}, for example.

In regular univariate time 
series modelling, users benefit
from both AIC and BIC, but are
aware of their shortcomings,
especially for short series. For example,
several different models
that yield similar AIC/BIC values make it
difficult to choose between them.
In addition, AIC and BIC computation
can be time consuming for models with
geometric parameter growth.
For example, for an order $p$
global-$\alpha$ GNAR model
that admits $s$ stages of neighbours per lag this
means estimating $p\,s+p$ parameters.
Typically, the `final' model will
be parsimonious, but $2^{p\,(s+1)}$ models might
need to be investigated with AIC/BIC
until that final model is found.
For instance, for monthly data, we
might want to start with $p=12$ to
enable detection of annual cycles and
with, e.g., $s=4$ this would result in
an unwieldy
$2^{12\times 5} \approx 10^{18}$ models to
be investigated.
\newpage
With regular univariate
time series,
analysts additionally 
have access 
to the well--known
tools of
the autocorrelation (acf) and partial
autocorrelation (pacf) functions,
which aid model order determination
and can also detect other behaviours
such as trend, seasonality and non-stationarity.
Indeed, acf and
pacf plots are often inspected before any formal modelling 
is attempted.
We introduce
a new  network--enabled version of these tools:
Section~\ref{sec:corbit} introduces both
the network autocorrelation function
(NACF) and the partial NACF (PNACF).
Section~\ref{sec:corbit} also introduces
a new graphical tool, the 
Corbit plot, to clearly convey
information in either
the NACF or PNACF
for a network time series
(Figure~\ref{fig: corbit sim} below shows
a preview).
Corbit plots help
users to quickly and
directly identify model
order and other characteristics of
network time series, just
as the acf and pacf
plots do for 
univariate series.
Section~\ref{sec:corbit} also exhibits
 Corbit plots on simulated data,
 clearly showing their advantage for model
 interpretation and selection.

Section~\ref{sec: graph model connection} proposes a  novel  interpretation of GNAR processes as generalised graphical models, which constrains the processes' autoregressive structure and proves
some interesting connections to graphical models by incorporating higher--order interactions.
Section~\ref{sec: graph model connection} 
also proposes a  new  interpretation of how including prior information into our analysis is related to performing variable selection and shrinkage for GNAR models. 
In particular, we prove a new result
explaining the connection
between our multi--stage GNAR
neighbourhood
structure and a hierarchy imposed
on the process'
inverse cross--spectrum matrix.
Section~\ref{sec: graph model connection} also
 shows  how classical graphical time
 series models can be seen as a
 special case of GNAR
 processes.
 

Section~\ref{sec: application}  proposes a new analysis of the problem of modelling, analysis and prediction of
the number of patients occupying mechanical ventilation beds
 (prevalence)  during the COVID--19 pandemic in $140$ 
 National Health Service (NHS) Trusts in England.
Corbit plots conveniently and rapidly give strong guidance as to choice of GNAR process order that suggests
a very parsimonious model.
We demonstrate the excellent predictive
performance for our
fitted GNAR models 
compared to established time series
models, but also their close connections to conditional
autoregressive spatial models, but only for the GNAR$(1, [1])$ model.

We further introduce an extension of the Corbit plot, namely
the R--Corbit plot, which permits
analysts
to understand (i) the 
effect of covariates on the
network time series correlation
structure and  (ii) how the
correlation structure can change
over different time periods.
The latter can provide a
clear and immediate indication
of nonstationarity, where it exists. For the COVID--19 mechanical  ventilation beds time series we find that a 
R--Corbit
plot can
show how the dynamics of the process change
during different waves of the pandemic.
The `R--' in the name
R--Corbit plot highlights that such plots are comprised of rings of circles.  We also present two new plots that show the local and global relevance
of individual NHS Trusts  within the network. 

GNAR is an alternative method for  time series modelling, which is especially useful when the data satisfy certain conditions, which we make explicit in Sections~\ref{sec:gnar methodology} and~\ref{sec: graph model connection}.
Thus, GNAR should be thought of as an addition to the existing toolbox for
multivariate time series or spatio--temporal data. 

Next, Section~\ref{sec:gnar methodology} formulates a hierarchical representation for GNAR processes, which allows us to write the model in compact matrix notation and later efficiently
define the GNAR NACF and PNACF and associated Corbit and R--Corbit plots.




%Body of paper
\input{GNAR_methods}

\input{gnar_model_selection}

\input{gnar_connection_graphical_models}

\input{application_covid_19}

\input{conclusion}

\section*{Acknowledgments}
We gratefully acknowledge the following support: Nason from EPSRC NeST Programme grant EP/X002195/1; 
Salnikov from the UCL Great Ormond Street Institute of Child Health, NeST
contribution from Imperial College London, 
the Great Ormond Street Hospital DRIVE Informatics Programme, and the Bank of Mexico.
Cortina--Borja supported by the NIHR Great Ormond Street Hospital Biomedical Research Centre. The views expressed are those of the authors and not necessarily those of the EPSRC, National Health Service (NHS), the NIHR or the UK Department of Health.


\newpage

\begin{center}
{\large\bf SUPPLEMENTARY MATERIAL IN APPENDICES}
\end{center}


\begin{description}

\item [A: Model Predictions] Tables of model prediction results.

\item [B: Further plots] Corbit plot with dashed line separating lag zero
    from highest lag. Rectangular Corbit plots.

\item [C: CARar and GNAR] Detailed comparison and notes.
\end{description}

\begin{center}
{\large\bf SUPPLEMENTARY MATERIAL IN SEPARATE DOCUMENT}
\end{center}

\begin{description}

\item[Correlation Structure Proof:]
Proof and illustration of Theorem \ref{th: pop_sparsity gnar}. (PDF/LaTeX)
\item[Network Autocorrelation Function Properties:]
Derivation and properties of the NACF given by Definition \ref{def: nacf}. (PDF/LaTeX)
\item[Algorithms and Properties of $\rstage$ Adjacency Matrices:]
Brief document exhibiting the properties and possible computation techniques for $\rstage$ adjacency matrices given by Definition \ref{def: rstage adjacency}. (PDF/LaTeX)
%\input{gnar_pop_sparse_proof}
\item[R-package for  Graphical Aids routine:] R-package \textbf{GNAR} and scripts containing code to perform the diagnostic methods described in the article. The package also contains all data sets used as examples in the article. (GNU zipped tar file). The GNAR package can
be accessed \href{https://cran.r-project.org/web/packages/GNAR/index.html}{\underline{HERE}}.
Further examples, further details on some of the
paper's examples and some new examples are presented in the following RMarkdown document:
\href{https://dansal182.github.io/new_tools_for_network_time_series_with_application_to_COVID_19_hospitalisations/network_autocorrelation_notebook.html}{\underline{R Markdown link: click here.}}

\item[COVID-19 Series data set:] Data set analysed in Section 5. The data set is available in the \texttt{GNAR} package, the file is \texttt{logMVbedMVC.vts}. A help file contains a detailed description of the data set and how it was obtained from the UK Government Coronavirus Dashboard
({\tt coronavirus.data.gov.uk}) between April 2020 and July 2021.
\end{description}

\appendix
\setcounter{equation}{0}
\counterwithin{equation}{section}

\input{model_prediction_appendix}


\section{Further Plots, Alternatives and Options}

\label{app: furtherplots}


\subsection{R-Corbit plot}
Figure~\ref{fig: R-Corbit pnacf} shows the PNACF R-Corbit plot referred to in the text
    above in Section~\ref{sec: model interpretation application}.
\begin{figure}
    \centering
   \hspace{-3.5cm} \includegraphics[scale=0.43]{images/R_Corbit_pnacf.pdf}
        \caption{COVID-19 series PNACF R-Corbit plot. The plot compares the PNACF values between the two COVID-19 outbreaks and the gap between them.
        The maximum lag is ten and the maximum $\rstage$ depth is equal to six.  See Section \ref{sec: model interpretation application} for description.}
    \label{fig: R-Corbit pnacf}
\end{figure}

\subsection{Corbit plot with `line gap'}
Figure~\ref{fig: corbit + line gap} shows a Corbit plot with a `line gap'. This reflects the
fact that lag~1 and lag~20 are not actually close, so a dotted line is drawn to indicate that
lag~20 is the last and highest lag and that the next lag in the circle (lag~1) is not a continuation.
Producing this is an option in the software and is the result of a suggestion of
Professor Berthold Lausen, University of Essex.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.60]{images/Corbit_gdpVTS_line_separation.pdf}
    \caption{NACF Corbit plot of the \textbf{gdpVTS} (included in the \texttt{GNAR} package). The plot includes the optional {\em line gap} argument, which highlights that lags are not a periodic quantity, and separates the first and last lags across all $r$-stages explicitly. }
    \label{fig: corbit + line gap}
\end{figure}

\subsection{Rectangular Corbit plots}
It is possible to format the NACF and PNACF in rectangular form, which might be more comfortable
for some users. Figures~\ref{fig:nacfrect} is the rectangular version of the GNAR$(2, [1, 1])$ simulation whose Corbit plot is shown in Figure~2 of the main article.
%Figure~\ref{fig:pnacfrect} is the rectangular version of a PNACF plot (not shown) from the same
%simulation.
Producing the rectangular versions is an option in the software.

\begin{figure}[H]
    \centering
    \resizebox{0.9\textwidth}{!}{\includegraphics{images/nacf_sim1_rectangular_version_two.pdf}}
    \caption{Rectangular NACF plot.}
    \label{fig:nacfrect}
\end{figure}

%\begin{figure}
%    \centering
%    \resizebox{0.8\textwidth}{!}{\includegraphics{images/pnacf_sim1_rectangular_version_two.pdf}}
%    \caption{Rectangular PNACF plot.}
%    \label{fig:pnacfrect}
%\end{figure}

\input{CARar_GNAR_comparison/CARar_GNAR_latex}

\input{glob_node_relevance_edit}

%\bibliographystyle{Chicago}
\bibliographystyle{guy3}
\bibliography{Bibliography-MM-MC}


\end{document}
