\section{NACF derivation}
\label{sec: derivation}
We are interested on estimating a generalised notion of autocorrelation, hence, we start by estimating a generalised notion of correlation. To do this we begin our analysis with respect to a multivariate Gaussian with common variation scale $\sigma^2$. We assume that there is a known positive definite matrix $\bmat{M} \in \mathbb{R}^{d \times d}$ that acts as a concentration matrix (inverse covariance) matrix for the vector $\boldsymbol{X} \in \mathbb{R}^d$ with respect to a mean $\boldsymbol{\mu}  \in \mathbb{R}^d$, however, we also assume that there is an unknown common variance scale $\sigma^2 > 0$  such that 
\begin{equation} \label{eq: deriv eq one}
    \boldsymbol{X} \sim \mathcal{N} (\boldsymbol{\mu}, \sigma^2 \bmat{M}^{-1} ) .
\end{equation}
From  \eqref{eq: deriv eq one} we have that the log-probability density function for $\boldsymbol{X}$ is 
\begin{equation} \label{eq: loglikelihood chi square}
    \log p(\boldsymbol{x}) = \log | \sigma^2 \bmat{M} | - \frac{1}{2 \sigma^2} \{ (\boldsymbol{x} - \boldsymbol{\mu})^{T} \bmat{M} (\boldsymbol{x} - \boldsymbol{\mu}) \} .
\end{equation}
Now, we focus on the Mahalanobis distance $\{ \sigma^{-2} (\boldsymbol{x} - \boldsymbol{\mu})^{T}  \bmat{M}  (\boldsymbol{x} - \boldsymbol{\mu}) \}^{1/2}$ that exploits the assumed structure for a known matrix $\bmat{M}$ and an unknown $\sigma^2$.  We will use the following proposition.
\begin{proposition} \label{prop: chi distr.}
Let $\boldsymbol{X} \in \mathbb{R}^d$ , such that $\boldsymbol{X} \sim \mathcal{N} (\boldsymbol{\mu}, \sigma^2 \bmat{M}^{-1} )$ and assume that $\bmat{M}$ is known. Then, using the residual random variables
$ \boldsymbol{\epsilon} := (\boldsymbol{X} - \boldsymbol{\mu}) $
we can define the scalar random variable $\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} $ and we have that 
\begin{equation*}
    \frac{\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon}}{\sigma^2} \sim \bigchi_{_{(d)}}^{2} .
\end{equation*}
\end{proposition}
\begin{proof}
Since $\bmat{M}^{- 1}$ is a symmetric positive definite matrix, we can find its Cholesky decomposition and factorise its inverse,
\begin{align} \label{eq: cholesky}
    &\bmat{M}^{-1} = \bmat{L} \bmat{L}^{T} , \nonumber \\
    &\bmat{M} = \big (\bmat{L}^{-1} \big )^{T} \bmat{L}^{-1} .
\end{align}
Moreover, there exists a random vector $\boldsymbol{u} \sim \mathcal{N}(0, \sigma^2 \bmat{I} )$ such that
\begin{equation} \label{eq: noise chol repar}
    \boldsymbol{\epsilon} = \bmat{L} \boldsymbol{u},
\end{equation}
and using \eqref{eq: cholesky} we have that
\begin{align} \label{eq: eps to chi square}
    \boldsymbol{u}^{T} \boldsymbol{u} &= \big ( \bmat{L}^{-1} \boldsymbol{\epsilon} \big )^{T} \big ( \bmat{L}^{-1} \boldsymbol{\epsilon} \big ), \nonumber \\
    &= \boldsymbol{\epsilon}^{T} (\bmat{L}^{-1} \big )^{T} \bmat{L}^{-1} \boldsymbol{\epsilon}, \nonumber \\
    &= \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} .
\end{align}
Now, we know that $\boldsymbol{u}^{T} \boldsymbol{u} / \sigma^2$ is a sum of $d$ independent squared standard Gaussian random variables, hence, $\boldsymbol{u}^{T} \boldsymbol{u} / \sigma^2$ has a chi-squared distribution with \mbox{$d$-degrees} of freedom. Therefore,  by substituting \eqref{eq: eps to chi square} into $\boldsymbol{u}^{T} \boldsymbol{u} / \sigma^2 \sim \bigchi_{_{(d)}}^{2}$ we have that
\begin{equation} \label{eq: end chi square}
     \frac{\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon}}{\sigma^2} \sim \bigchi_{_{(d)}}^{2} .
\end{equation}
\end{proof}

Proposition \ref{prop: chi distr.} shows that the generalised dot product for the residuals $\boldsymbol{\epsilon}$ has a chi-squared distribution with \mbox{$d$-degrees }of freedom, thus, we can form the likelihood function with the chi-squared density function,
\begin{equation} \label{eq: log-likelihood sigma }
    p(\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} | \sigma^2 ) = {\bigg \{ (2 \sigma^2)^{d / 2} \Gamma \bigg (\frac{d}{2} \bigg ) \bigg \} }^{-1}
    \bigg \{ \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} \bigg \}^{\frac{d}{2} - 1} 
    \exp \bigg( - \frac{\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} }{2 \sigma^2} \bigg ) .
\end{equation}
Now, we can think of $\boldsymbol{\epsilon}$ as $d$ observations for estimating the quantity $\sigma^2$, so we compute the maximum likelihood estimator using \eqref{eq: log-likelihood sigma }. 
\begin{equation*}
   \dfrac{d}{d \sigma^2} \ln L ( \sigma^2| \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} ) = -\frac{d}{2} \ln (2 \sigma^2 ) - 
   \bigg ( \frac{d}{2} - 1 \bigg ) \ln  \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} - 
   \frac{1}{2 \sigma^2}  \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} + k ,
\end{equation*}
note that $k$ is constant with respect to $\sigma^2$, hence, 
\begin{equation} \label{eq: mle sigma}
    \hat{\sigma}^2_{ML} = \frac{ \boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon} }{d} .
\end{equation}
Above, \eqref{eq: mle sigma} proposes an estimator for the unknown variance $\sigma^2$ if we have knowledge of the matrix $\bmat{M}$. This suggests comparing the covariance matrix $\sigma^2 \bmat{M}^{-1}$ with the isotropic case $\sigma^2 \bmat{I_d}$. Suppose that we can find a constant $\nu > 0$ such that for all $\boldsymbol{x} \in \mathbb{R}^d$ we have that 
\begin{equation*}
    \normx{\sigma^2 \bmat{M} \boldsymbol{x}}_2 \leq \sigma \nu \normx{\boldsymbol{x}}_2 .
\end{equation*}
Thus, we can define the generalised correlation as 
\begin{equation} \label{eq: network correlation}
    \rho_{\nu, \bmat{M}} := \frac{\boldsymbol{\epsilon}^{T} \bmat{M} \boldsymbol{\epsilon}}{
                \nu \normx{\boldsymbol{\epsilon}}^2_2
    } .
\end{equation}
Above, \eqref{eq: network correlation} computes the ratio between the estimator for $\sigma^2$ given by \eqref{eq: mle sigma} and the classic estimator $\normx{\boldsymbol{\epsilon}}^2_2 / d$. 
Furthermore, we note that if there is no partial covariance between two variables (i.e., $[\bmat{M}]_{ij} = 0$), then we do not have to bound said interaction when finding $\nu$. Thus, we see that a possible value for $\nu$ is an upper bound on the interactions given by $\bmat{M}$.  This can be done by bounding the pair-wise covariances between all possible pairs.  Note that defining $\nu$ as this results in a quantity lower bounded by the operator norm associated to the matrix $\bmat{M}$.
\par
Now, assume that $\boldsymbol{X}_t$ is a zero mean network time series that arises from a known network with graph $\mathcal{G}$ and that $\bmat{W}_{t}$ is the corresponding weight matrix. Note that the off-diagonal entries are the connection weights which quantify the influence each node has on other nodes in the network. Hence, if $[\bmat{W}_{t}]_{ij} = 0$, then $X_{t + h, j}$ is not an \mbox{$r$-stage} neighbour considered for neighbourhood regression at any lag with respect to node $X_{i, t}$. Thus, following the ideas from above we can use the following concentration matrix
\begin{equation} \label{eq: network structure matrix}
    \bmat{M}_t = \big ( \bmat{W}_t \odot \bmat{S}_r + \bmat{I_d} \big ) .
\end{equation}
This allows us to compute the variation that is common across the network by considering the linear path given by the matrix in \eqref{eq: network structure matrix} for each node. 
Furthermore, assume that we can find $\lambda_t$ such that for any choice of \mbox{$h$-lag} and \mbox{$r$-stage} neighbour the network covariance is upper bounded by the autocovariances
\begin{equation*}
    E \{ X_{i, t + h} \langle [\mathbf{W}_{t}]_{i\cdot}, \boldsymbol{X}_{t} \rangle \} \leq
    \lambda_t \sum_{i = 1}^{d} E ( X_{i, t + h} X_{i, t} ) \leq  \lambda_t \sum_{i = 1}^{d} E ( X_{i, t}^2 ).
\end{equation*}
This suggests that we can think of a common scale $\sigma^2$ for the network time series with matrix $\big ( \bmat{W}_t \odot \bmat{S_r} + \bmat{I_d} \big )$, and get a formula motivated by \eqref{eq: network correlation}. Now, assume that $E ( \boldsymbol{X}_{t} ) = \boldsymbol{\mu} < +\infty$ for all $t$ , and consider the quantity
\begin{equation} \label{eq: bounded covariance}
    E \{ (  \boldsymbol{X}_{t + h} - \boldsymbol{\mu})^{T} \big ( \bmat{W}_t \odot \bmat{S}_r + \bmat{I_d} \big )
    ( \boldsymbol{X}_{t} - \boldsymbol{\mu}) \} .
\end{equation}
Above, \eqref{eq: bounded covariance} calculates the sum of the expected values of covariation between each nodal time series with respect to itself and its \mbox{$r$-stage} neighbours at the $k$th-lag. Thus, we can use $\lambda_t$ to bound the expectation in \eqref{eq: bounded covariance}, and for a network time series $\boldsymbol{X}_t$ paired to the network $\mathcal{G}$ we have that
\begin{equation} \label{eq: nacf theoretical}
     - 1 \leq \frac{E \{ (  \boldsymbol{X}_{t + h} - \boldsymbol{\mu})^{T} \big ( \bmat{W}_t \odot \bmat{S}_r + \bmat{I}_t \big )
    ( \boldsymbol{X}_{t} - \boldsymbol{\mu}) \}}{ E \{ (  \boldsymbol{X}_{t} - \boldsymbol{\mu})^{T}
    (1 + \lambda_t)\bmat{I_d} ( \boldsymbol{X}_{t} - \boldsymbol{\mu}) \}} \leq 1 .
\end{equation}
We can estimate $\boldsymbol{\mu}$ with the sample mean $\overline{\boldsymbol{X}}$. Combining these ideas, we define the sample network autocorrelation function by substituting the observations $\{ \boldsymbol{X}_t \} $ into \eqref{eq: nacf theoretical} and estimating the expectations with sample means (i.e., averaging the observed values). Said procedure results in the NACF given by \eqref{eq:nacf definition}.
