\section*{Proof of Theorem \ref{th: pop_sparsity gnar}}
Assume that $\boldsymbol{X}_t$ is a stationary $\GNAR(p, [s_1, \cdots, s_p])$ process for which the underlying network structure $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ is static and with full rank spectral matrix. Throughout this section the following notation and definitions are used 
\begin{itemize} \label{def: notation and sets}
    \item $\boldsymbol{X}_t$ multivariate time series.
    \item $r^* = \max(s_1, s_2, \cdots, s_p)$ the maximum $\rstage$ active depth.
    \item $\cup_{r = 1}^{r^*} \mathcal{N}_{r} (i) =: \mathcal{N} (i)$ the set of active $\rstage$ neighbors with respect to node $i$.
    \item $\mathcal{A}_r := (i, j) \in \mathcal{A}_r \iff \exists l \in \mathcal{K}$ such that $\delta(i, l) \leq r$ and $\delta(l, j) \leq r;$ 
    the set of active $\rstage$ pairs in neighborhood regressions.
    \item $\mathcal{P}_{r} := \mathcal{A}_{r} - \mathcal{A}_{r - 1}$ the set of pairs which activate exactly at a distance equal to $r$ 
    \newline (i.e., $(i, j) \in \mathcal{P}_r \iff \exists l \in \mathcal{K}$ such that $\delta(i, l) \leq r$, $\delta(l, j) \leq r$ and $(i, j) \notin \mathcal{A}_{r - 1}$).
    \item $\bmat{f}(\omega)$ spectral matrix of $\boldsymbol{X}_t$.
    \item $\bmat{S}(\omega) = \{\bmat{f}(\omega)\}^{-1}$ inverse spectral matrix of $\boldsymbol{X}_t$. 
    \item $f(\cdot) \equiv 0$ the function $f$ is identically zero in its domain.
    \item $f_{i, j | \mathcal{A}^{(ij)}}(\omega)$ conditional cross-spectrum between $X_{i, t}$ and $X_{j, t}$.
    \item $X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)}$ the two time series are independent given all the other ones in $\mathcal{A}^{(ij)}$.
    \item $\{\bmat{S}_{\mathcal{A}^{(ij)}} (\omega) \}$ inverse spectral matrix constrained to the set $\mathcal{A}^{(ij)}$.
    \item $ \boldsymbol{f}_{i, \mathcal{A}^{(ij)}} (\omega)$ $i$th column of the spectral matrix without the $i$th entry.
\end{itemize}
%Proof starts here
\newpage
\begin{proposition} \label{prop: distance for conditional independence}
    Assume that $\boldsymbol{X}_t$ is a stationary $\GNAR(p, [s_1, \cdots, s_p])$ process for which the underlying network structure $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ is static and with the same notation as above, then $(i, j) \notin \mathcal{A}_{r^*}$ if and only if $\delta(i, j) \geq 2 r^* + 1$. And, that $(i, j) \in \mathcal{P}_{r^*}$ if and only if $\delta(i, j) \in \{2 r^* - 1, 2 r^* \}$.
\end{proposition}
\begin{proof}
    Let $(i, j) \notin \mathcal{A}_{r^*}$, then there does not exist $l \in \mathcal{K}$ such that
    \begin{equation} \label{eq: borough active set condition}
        \delta(i,l) \leq r^* \textit{ and  } \delta(j, l) \leq r^* .
    \end{equation}
    Using (\ref{eq: borough active set condition}) note that if there is at least one $l \in  \mathcal{N} (i) \cap  \mathcal{N} (j)$, then such $l$ satisfies $\delta(i, l) \leq r^*$ and $\delta(j, l) \leq r^*$, and by definition we have that $(i, j) \in \mathcal{A}_{r^*}$, thus, $(i, j) \notin \mathcal{A}_{r^*}$ if and only if $\mathcal{N} (i) \cap  \mathcal{N} (j) = \varnothing$.
    \par
    Now, suppose that $\delta(i, j) \leq 2 r^*$, then we can choose $l$ such that 
        $\delta(i, l) = r$,
    for some $r \in \{1, \cdots, r^* \}$, hence, $l \in \mathcal{N} (i)$. Note that since $\delta(i, j) > r^*$, we can write
    \begin{align*}
        \delta(i, j) &= \delta(i, l) + \delta(j, l) \\
        &= r + \delta(j, l) \leq 2 r^* \\
        &\iff \delta(j, l) \leq r^*,
    \end{align*}
    so we have that $l \in \mathcal{N} (j)$, and 
    \begin{equation*}
        l \in  \mathcal{N} (i) \cap  \mathcal{N} (j) \dagger
    \end{equation*}
    a contradiction since $(i, j) \notin \mathcal{A}_{r^*}$, therefore, 
        $\delta(i, j) \geq 2 r^* + 1$.
    Now, decrease $r^*$ to $(r^* -1)$, since $(i, j) \notin \mathcal{A}_{r^* - 1}$ if and only if $\delta(i, j) \geq 2 r^* - 1$, we have that $\mathcal{A}_{r^* - 1} \subseteq \mathcal{A}_{r^*}$, and 
    $$2 r^* - 1 \leq \delta(i, j) < 2 r^* + 1 \iff (i, j) \in \mathcal{A}_{r^*} - \mathcal{A}_{r^* - 1},$$
    hence, $(i, j) \in \mathcal{P}_{r^*}$ if and only if $\delta(i, j) \in \{2 r^* - 1, 2 r^* \}$.
\end{proof}
\newpage
Recall from time series analysis that the best linear predictor up to a certain lag restricted to a subset of active nodes $\mathcal{A} \subseteq \mathcal{K}$ solves for the optimal linear filter with respect to mean squared error. Set $\mathcal{A}^{(ij)} := \mathcal{K} - \{i, j\}$, then the error of the optimal filter conditioning on $\mathcal{A}^{(ij)}$ is given by

\begin{equation} \label{eq: opt filter}
   u_{i, t}^{(ij)} := u_{i, t | \mathcal{A}^{(ij)} } = X_{i, t} - a_{i} - \sum_{k = 1}^{p} \sum_{l \in \mathcal{A}_{r^*}^{(j)}} 
   \gamma_{k}^{il} X_{l, t - k},
\end{equation}
where $\mathcal{A}_{r^*}^{(j)} = \mathcal{A}^{(ij)} \cap \mathcal{N} (i)$  and $a_i$ is the optimal mean constant; see \cite{brockwell_davis}. 
Following \cite{graph_series}, define the conditional dependence structure based on the cross-spectrum between the conditional time series $X_{i, t} | \mathcal{A}^{(ij)}$ and $X_{j, t} | \mathcal{A}^{(ij)}$ by noting that these two time series are uncorrelated at all lags if and only if the conditional cross-spectrum $f_{i, j | \mathcal{A}^{(ij)}}(\omega)$ is identically zero for all frequencies $\omega \in [-\pi, \pi]$. We will make use of the following result proved in \cite{graph_series}.
\begin{theorem}{[\cite{graph_series}]} \label{th: inv spec cross indentity}
Let $\boldsymbol{X}_t$ be a stationary multivariate time series, if the spectral matrix $\bmat{f}(\omega)$ has full rank, then
\begin{align*}
    f_{i, j | \mathcal{A}^{(ij)}}(\omega) &:= f_{u_{i}^{(ij)} u_{j}^{(ij)}} (\omega), \\
    f_{i, j | \mathcal{A}^{(ij)}}(\omega) &= - [\bmat{S} (\omega)]_{ij}.
\end{align*}
\end{theorem}
The relevance of Theorem \ref{th: inv spec cross indentity} is that it enables us to identify conditional independence by looking at the inverse spectral matrix $\bmat{S} (\omega)$, rather than by computing all possible conditional cross-spectrum values. Hence, the inverse spectral matrix is related to multivariate time series models as the concentration (inverse-covariance) matrix is to Gaussian Graphical Models; see \cite{graph_lasso}. 
This allows to state results from \cite{graph_series} and \cite{brillinger} in the following Lemma. 
\newpage
\begin{lemma}{[\cite{graph_series}, and \cite{brillinger} ]} \label{lma: conditional equivalences}
Let $\boldsymbol{X}_{t} \in \mathbb{R}^d$ be a stationary multivariate process, with the same notation as above, then 
\begin{align} \label{eq: conditional independence equivalence}
    X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ for all } h  &\iff \cov (u_{i, t}^{(ij)}, u_{j, t + h}^{(ij)} ) \equiv 0, \nonumber \\
    &\iff f_{i, j | \mathcal{A}^{(ij)}}(\omega) \equiv 0, \nonumber \\
    &\iff [\bmat{f}(\omega)]_{ij} - \boldsymbol{f}_{i, \mathcal{A}^{(ij)}} (\omega) \{ \bmat{S}_{\mathcal{A}^{(ij)}} (\omega) \}^{-1}  
    \boldsymbol{f}_{\mathcal{A}^{(ij)}, j} (\omega), \nonumber \\
    &\iff [\bmat{S} (\omega)]_{ij} \equiv 0.
\end{align}
\end{lemma}
A partial correlation graph satisfies the Markov properties if we define the absence of an edge as the relation $ X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all } h$, which makes explicit the relationship between nodes $i$ and $j$ (i.e., $i$ and $j$ are uncorrelated at all lags given all the other nodes in the graph); see \cite{graph_series}. GNAR generalizes this idea by noting that the conditional independence can be extracted from paths in the graph rather than just from existing edges. 

\begin{definition} \label{def: conditional indep definition}
Let $\boldsymbol{X}_t$ be a stationary multivariate process with an underlying $\GNAR$ structure $\mathcal{G} = (\mathcal{K}, \mathcal{E})$, the $\GNAR$ generalized correlation graph $G = (\mathcal{K}, E)$ is an undirected graph which satisfies $d_{E} (i, j) > r^*$ if and only if nodes $i$ and $j$ are independent given all the other nodes, which we write as
\begin{equation*}
     X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all } h \iff d_{E} (i, j) > r^*.
\end{equation*}
Furthermore, the strength of dependence can be recovered from $G$ in the sense that if $d_{E} (i, j) = r$ for some $r \in \{1, \cdots, r^* \}$, then constraining the $\GNAR$ model to maximum $\rstage$ depth of $(r - 1)$ will uncorrelate the two nodal time series (i.e., $ X_{i, t} \perp X_{j, t + h}$ at all lags $h$).
\end{definition}
\newpage
Intuitively, Definition \ref{def: conditional indep definition} proposes that two nodes in the underlying network $\mathcal{G}$ are conditionally independent if and only if the shortest path between their active $\rstage$ neighbors has a length of at least one (i.e., they do not have common active $\rstage$ neighbors). This condition expresses mathematically that no active projection path is long enough to contain both $i$ and $j$.
\par
Before continuing we mention the following result from \cite{graph_series}. 
\begin{proposition}{[\cite{graph_series}]} \label{prop: conditional residual}
Let $\boldsymbol{X}_t$ be a multivariate time series with correlation structure given by $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ and spectral matrix $\bmat{f}(\omega)$ which has full rank for all frequencies $\omega \in [-\pi, \pi]$, then
    $$    (i, j) \notin \mathcal{E}  \iff [\bmat{S}(\omega)]_{ij} \equiv 0.
    $$
\end{proposition}
The following proposition is necessary as well.
\begin{proposition} \label{prop: var constrains}
Let $\boldsymbol{X}_t$ be a stationary $\GNAR$ process with full rank spectral matrix $\bmat{f} (\omega)$, static underlying network $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ structure and maximum active $\rstage$ neighborhood regression $r^*$, then we have that
\begin{equation*}
    (i, j) \notin \mathcal{A}_{r^*} \textit{ if and only if } X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all lags } h.
\end{equation*}
\end{proposition}
\begin{proof}
Recall the VAR representation of  $\GNAR(p, [s_1, \cdots, s_p])$ process given by (\ref{eq: sparse var}),
\begin{equation*}
    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*}
for which the autoregressive matrices are given by 
$$\bmat{\Phi}_k = \{\diag(\alpha_{ik}) + \sum_{r = 1}^{s_k} \beta_{kr} (\bmat{W} \odot \bmat{S}_r ) \},$$
and the entries satisfy
\begin{equation*}
    [\bmat{\Phi}_k]_{ij} = \alpha_{ik} \mathbb{I}(i = j) + \beta_{k r} w_{i j} \mathbb{I}(i \neq j),
\end{equation*}
\newpage
moreover, since $\phi_{k}^{ij} = \beta_{kr} w_{ij}$ and $\phi_{k}^{ji} = \beta_{kr} w_{ji}$ we have that
\begin{equation} \label{eq: ar coeffss equiv}
    \phi_{k}^{ij} \equiv 0 \iff \beta_{kr} \equiv 0 \iff \phi_{k}^{ji} \equiv 0,
\end{equation}
regardless of the connection weights between nodes, it is active $\rstage$ depth which activates a connection between two nodes. And that the characteristic matrix polynomial for $z \in \mathbb{C}$ is given by
$$ \bmat{\Phi}(z) := \bmat{I_d} - \sum_{k = 1}^{p} \bmat{\Phi}_k z^k, $$
and recall that $\bmat{S} (\omega)$ satisfies; see \cite{brockwell_davis},
\begin{equation} \label{eq: inverse spec coeffs}
    \bmat{S} (\omega) = \frac{2 \pi}{\sigma^2} \{\bmat{\Phi}(e^{\imath \omega})\}^T \bmat{\Phi}(e^{-\imath \omega}).
\end{equation}
Now combine (\ref{eq: inverse spec coeffs}) and (\ref{eq: ar coeffss equiv}), so we have that
\begin{equation*}
    [\bmat{S} (\omega)]_{ij} = \frac{2 \pi}{\sigma^2} \sum_{l = 1}^{d} [\{\bmat{\Phi}(e^{\imath \omega})\}^T]_{il} [\bmat{\Phi}(e^{-\imath \omega})]_{lj},
\end{equation*}
hence, we also have that
\begin{equation} \label{eq: inv spec and ar coefs}
    [\bmat{S} (\omega)]_{ij} \equiv 0 \iff  \sum_{l = 1}^{d} [\{\bmat{\Phi}(e^{\imath \omega})\}^T]_{il} [\bmat{\Phi}(e^{-\imath \omega})]_{lj} \equiv 0.
\end{equation}
Assume that $[\{\bmat{\Phi}(e^{\imath \omega})\}^T]_{il} [\bmat{\Phi}(e^{-\imath \omega})]_{lj} \equiv 0$, then the sum in the right hand side of (\ref{eq: inv spec and ar coefs}) is equal to zero, and by combining Lemma \ref{lma: conditional equivalences} and (\ref{eq: inv spec and ar coefs}) we have that
\begin{align} \label{eq: inactive to condtitional independende}
    [\{\bmat{\Phi}(e^{\imath \omega})\}^T]_{il} [\bmat{\Phi}(e^{-\imath \omega})]_{lj} \equiv 0 &\implies  [\bmat{S} (\omega)]_{ij} \equiv 0, \nonumber \\
    &\iff X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all lags } h.
\end{align}
Now, suppose that $(i, j) \notin \mathcal{A}_{r^*}$, then $\delta(i, j) > r^*$ and there does not exist $l \in \mathcal{K}$ such that $\delta(i, l) \leq r^*$ and 
$\delta(j, l) \leq r^*$, where the distances are computed on the underlying network $\mathcal{G}$. 
\newpage
Thus, for all $l \in \mathcal{A}^{(ij)}$ we have that
\begin{align*} \label{eq: first prop cond independence}
    \delta(i, l) \leq r^* &\implies \delta(j, l) > r^*, \nonumber \\
    &\implies \phi_{k}^{lj} \equiv 0,
\end{align*}
then since $l$ and $j$, using the above and (\ref{eq: ar coeffss equiv}), do not influence each other we have that
\begin{equation*}
    [\bmat{\Phi}(e^{-\imath \omega})]_{lj} \equiv 0 \implies [\{\bmat{\Phi}(e^{\imath \omega})\}^T]_{il} [\bmat{\Phi}(e^{-\imath \omega})]_{lj} \equiv 0.
\end{equation*}
Note that the argument above holds similarly if $ \delta(j, l) \leq r^*$, and that if $\delta(i, l) > r^*$ and $\delta(j, l) > r^*$, then we have that 
$\phi_{k}^{lj} \equiv 0 \equiv \phi_{k}^{li},$
so by (\ref{eq: inactive to condtitional independende}) and the arguments above we have that
\begin{equation}  \label{eq: first prop cond independence end}
    (i, j) \notin \mathcal{A}_{r^*} \implies X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all lags } h.
\end{equation}
We proceed to prove the other direction, assume that $(i, j) \in \mathcal{A}_{r^*}$, then there exists at least one other nodal time series such that its corresponding node $l \in \mathcal{K}$, $\delta(i, l) \leq r^*$ and  $\delta(j, l) \leq r^*$. Hence, by (\ref{eq: inv spec and ar coefs}, \ref{eq: inverse spec coeffs}, \ref{eq: ar coeffss equiv}) we have that 
\begin{equation*}
    (i, j) \in \mathcal{A}_{r^*} \implies  [\bmat{S} (\omega)]_{ij} \not\equiv 0,
\end{equation*}
which by Lemma \ref{lma: conditional equivalences} gives
\begin{equation} \label{eq: second prop cond independence end}
     (i, j) \in \mathcal{A}_{r^*} \implies X_{i, t} \not\perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all lags } h.
\end{equation}
Finally, combining (\ref{eq: first prop cond independence end}, \ref{eq: second prop cond independence end}) finishes the proof.
\end{proof}
\newpage
Now it is possible to prove Theorem \ref{th: pop_sparsity gnar}, we start by proving the first part (i.e., $a$).
\begin{proof}
Consider the stationary GNAR process $\boldsymbol{X}_t$, we can embed a latent graphical process $\tilde{\boldsymbol{X}}_t$ to this GNAR process by connecting all nodes which fail to satisfy 
\newline $X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all lags } h$. And, 
%And, note that the optimal filters given by (\ref{eq: opt filter}) are uncorrelated at all lags $h$ if and only if nodes $i$ and $j$ do not collide in an active projection path, then 
that by Lemma \ref{lma: conditional equivalences} and Theorem \ref{th: inv spec cross indentity} we have that 
\begin{equation} \label{eq: path length equivalence}
    X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ for all } h \iff [\bmat{S}(\omega)]_{ij} \equiv 0 .
\end{equation}

Now, denote by $\tilde{\mathcal{G}} = (\mathcal{K}, \tilde{\mathcal{E}})$ the corresponding partial correlation graph. And, connecting (\ref{eq: path length equivalence}), Theorem \ref{th: inv spec cross indentity} and Proposition \ref{prop: conditional residual} we have that 
\begin{align} \label{eq: end first part}
     (i, j) \notin \tilde{\mathcal{E}} &\iff  [\bmat{S}(\omega)]_{ij} \equiv 0 \nonumber \\
     &\iff X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)} \textit{ at all } h \nonumber \\
     &\iff \cov (u_{i, t}^{(ij)}, u_{j, t + h}^{(ij)} ) \equiv 0,
\end{align}
Then by Proposition \ref{prop: var constrains} and (\ref{eq: end first part}) we have that
\begin{equation*}
      (i, j) \notin \tilde{\mathcal{E}} \iff (i, j) \notin \mathcal{A}_{r^*},
\end{equation*}
and by Proposition \ref{prop: distance for conditional independence}, we have that
\begin{equation*}
     (i, j) \notin \tilde{\mathcal{E}} \iff \delta(i, j) \geq 2 r^* + 1,    
\end{equation*}
which establishes the first part of Theorem \ref{th: pop_sparsity gnar}.
\end{proof}

Before proving the second part we mention the following notions, and recall that our underlying motivation and assumption is that the larger $r$ is the weaker the dependence between $\rstage$ neighbors is. First, we introduce the complex soft-threshold operator.
\newpage
\begin{definition} \label{def: complex sft-thrsh}
    The complex soft-threshold operator is given by
    \begin{align} \label{eq: cmplx sft}
        \sft(z; \rho) : \mathbb{C} &\rightarrow \mathbb{C} \nonumber \\
        z &\mapsto (|z| - \rho)_{+} e^{\imath \phi},
    \end{align}
    where $\phi = \arg(z) \in [0, 2 \pi]$, $\rho \geq 0$, $\imath$ denotes the imaginary unit (i.e., $\imath^2 = -1$), $(\cdot)_{+}$ denotes the positive part operator, and recall that $|z| e^{\imath \phi} = z$. 
\end{definition}
Note that if the argument $\phi \in \{0, \pi \}$, then $\sft(z; \rho)$ recovers the real-valued soft-threshold operator.
Next, we state a mathematical form of distinguishing between different degrees of relevance between nodes. Motivated by distance as a measure of dissimilarity, assume that there exists a cross-correlation hierarchy and define $\rstage$ neighborhood inclusion, with maximum active stage $r^*$, as follows.
\begin{definition} \label{def: rstage hierarchy}
    Let $\boldsymbol{X}_t$ be a stationary $\GNAR$ process with static underlying structure $\mathcal{G} = (\mathcal{K}, \mathcal{E})$ and conditional cross-correlation hierarchy 
    $0 < \tau^{(r^*)}_h < \cdots < \tau^{(r)}_h < \tau^{(r - 1)}_h < \cdots < \tau^{(1)}_h$, then  $(i, j) \in \mathcal{A}_r \textit{ if and only if }$ the following conditions hold
    \begin{align*}
      \sft (| \cov (u_{i, t}^{(ij)}, u_{j, t + h}^{(ij)} )|;  \tau^{(r - 1)}_h) &\equiv 0, \\
      \sft (| \cov (u_{i, t}^{(ij)}, u_{j, t + h}^{(ij)} )|;  \tau^{(r)}_h) &\not\equiv 0.
    \end{align*}
\end{definition}
Definition \ref{def: rstage hierarchy} makes explicit our assumption that the further nodes are in the graph, the less dependent they are. It enables us to reparameterize a GNAR process by punching out weak dependencies just by decreasing the length of maximum $\rstage$ neighborhood regression (i.e., if we increase the conditional cross-correlation threshold, then nodes further from each other become conditionally independent). 
\newpage
A possible way of interpreting Definition \ref{def: rstage hierarchy} is that it makes explicit the idea of a correlation index by distinguishing between different thresholds $\tau^{(r)}_h$ of relevance for which the underlying GNAR structure $\mathcal{G}$ captures the different strengths.
\par
A useful consequence of Definition \ref{def: rstage hierarchy} is that if the cross-correlation hierarchy exists, then by Theorem \ref{th: inv spec cross indentity} there exists a cross-spectrum hierarchy 
$\infty > \rho^{(1)} >  \rho^{(2)} > \cdots > \rho^{(r^*)} >  \rho^{(r^* + 1)} = 0 $ such that
\begin{equation*}
    \tau^{(r)}_h \leq |\cov (u_{i, t}^{(ij)}, u_{j, t + h}^{(ij)} )| < \tau^{(r - 1)}_h \iff 
     \rho^{(r)} \leq |[\bmat{S}(\omega)]_{ij}| < \rho^{(r - 1)},
\end{equation*}
for all $h$ and at all $\omega \in [-\pi, \pi]$, thus, defining $\rstage$ neighborhood inclusion as in Definition \ref{def: rstage hierarchy} guarantees the existence of a cross-spectrum hierarchy.
\par
Now, it is possible to define an $\rstage$ sequence of inverse spectral matrices as
\begin{equation} \label{eq: r stage inv matrix}
    \bmat{S}^{(r)}(\omega) := \sft \big( \bmat{S} (\omega) + \rho^{(r)} \diag\{e^{\imath \phi_{ii}}\};  \rho^{(r)} \big),
\end{equation}
where $\phi_{ii}$ is the argument of the diagonal entries in $\bmat{S} (\omega)$ and $\sft(z; \rho)$ is applied entry-wise to $\bmat{S} (\omega)$. Note that the diagonal of $\bmat{S} (\omega)$ does not change for any $r \in \{1, \cdots, r^*\}$.
And, that the sequence satisfies 
$$\normx{\bmat{S}(\omega) }_0 < \normx{\bmat{S}^{(r)}(\omega) }_0,$$ 
for all active $r$, and since $\rho^{(r + 1)} < \rho^{(r)}$, also
$$\normx{\bmat{S}^{(r + 1)}(\omega) }_0 \leq \normx{\bmat{S}^{(r)}(\omega) }_0,$$
where $\normx{ \cdot }_0$ counts the number of non-zero entries in the matrix. The sequence above makes explicit the different strengths of dependence between nodes, as $r$ decreases the $\bmat{S}^{(r)}(\omega)$ become sparser which results in a model with fewer nodes influencing each other; see Section \ref{sec: graph model connection}.
\newpage
We are ready to prove the second part (i.e., $b$) of Theorem \ref{th: pop_sparsity gnar}.
\begin{proof}
Assume that $\bmat{f} (\omega)$ has full rank and that the sequence of $\{\bmat{S}^{(r)} (\omega) \}$ are non-singular. Define the sequence of partial correlation graphs 
$\mathcal{G}^{(r)} = (\mathcal{K}, \mathcal{E}^{(r)})$ as the one which satisfies $(i, j) \notin \mathcal{E}^{(r)}$ if and only if $X_{i, t} \perp X_{j, t + h} | \mathcal{A}^{(ij)}$ at all lags $h$ when constraining $\rstage$ depth to $r \in \{1, \cdots, r^*\}$.
By Theorem \ref{th: inv spec cross indentity} and Proposition \ref{prop: var constrains}
we have that
\begin{equation} \label{eq: r sequence equivalence}
    (i, j) \notin \mathcal{E}^{(r)}  \iff [\bmat{S}^{(r)}(\omega)]_{ij} \equiv 0, 
\end{equation}
Using (\ref{eq: r sequence equivalence}) we have that the sequence of edge sets satisfies
\begin{equation} \label{eq: edge set sequence}
     \mathcal{E}^{(1)} \subset \cdots \subset \mathcal{E}^{(r)} \subset \mathcal{E}^{(r + 1)} \subset \cdots  \subset \tilde{\mathcal{E}},
\end{equation}
where $\tilde{\mathcal{E}}$ is defined as in (\ref{eq: end first part}). Note that once a pair of edges become conditionally independent they stay so since the soft-threshold operator maps that entry to zero for all larger $\rho^{(r)}$ in the cross-spectrum hierarchy, this idea allows us to prune the underlying graph by removing edges based on the strength (i.e., magnitude) of the conditional cross-spectrum $f_{i, j | \mathcal{A}^{(ij)}}(\omega)$; see Lemma \ref{lma: conditional equivalences}. Now consider a pair of nodes $i$ and $j$ which satisfy
\begin{equation} \label{eq: second part heart}
    \rho^{(r)} \leq |[\bmat{S}(\omega)]_{ij}| < \rho^{(r - 1)},
\end{equation}
for all $\omega \in [-\pi, \pi]$, then, applying the soft-threshold given by Definition \ref{def: complex sft-thrsh} and (\ref{eq: cmplx sft}, \ref{eq: r sequence equivalence}) to this pair results in 
$(i, j) \in \mathcal{E}^{(r)}$, and $(i, j) \notin \mathcal{E}^{(r - 1)}$, thus, we have that 
\begin{align} \label{eq: second part inv equivalences}
    &(\ref{eq: second part heart}) \implies (i, j) \in \mathcal{E}^{(r)} - \mathcal{E}^{(r - 1)}, \nonumber \\
    &(i, j) \in \mathcal{E}^{(1)} \implies (i, j) \in \bigcap_{r = r^*}^{1} \mathcal{E}^{(r)}, \nonumber \\
    &(i, j) \in \mathcal{E}^{(r)} \implies (i, j) \in \bigcap_{s = r^*}^{r} \mathcal{E}^{(s)} .
\end{align}
\newpage
Following the ideas above, consider the graph $\mathcal{G}^{(1)} = (\mathcal{K}, \mathcal{E}^{(1)})$ associated to the strongest connections with edge set $\mathcal{E}^{(1)}$, which is the one with the strongest dependence between nodes as proposed in Definition \ref{def: rstage hierarchy}. 
\par
Using (\ref{eq: second part heart}, \ref{eq: second part inv equivalences}) we have that if $(i, j) \in \mathcal{E}^{(1)}$, then $i$ and $j$ are connected throughout the whole sequence in (\ref{eq: edge set sequence}), and that the edge set sequence satisfies
\begin{align} \label{eq: second part inverse end}
    \mathcal{A}_{1} &= \mathcal{E}^{(1)}, \nonumber \\
    \mathcal{A}_{2} - \mathcal{A}_{1} &= \mathcal{E}^{(2)} - \mathcal{E}^{(1)}, \nonumber \\
    &\cdot \nonumber \\
    \mathcal{A}_{r} - \mathcal{A}_{r - 1} &= \mathcal{E}^{(r)} - \mathcal{E}^{(r - 1)}, \nonumber \\
    &\cdot \nonumber \\
    \mathcal{A}_{r^*} - \mathcal{A}_{r^* - 1} &= \mathcal{E}^{(r^*)} - \mathcal{E}^{(r^* - 1)}.
\end{align}
Moreover, computing the shortest path between nodes $i$ and $j$ on the graph $\mathcal{G}^{(1)}$ can be done as follows. 
\par
If $(i,j) \in \mathcal{E}^{(1)}$, then $\delta(i,j) = 1$, else if $(i,j) \in \mathcal{E}^{(2)} - \mathcal{E}^{(1)}$, then since $(i,j)$ are $1-$stage neighbors in $\mathcal{G}^{(2)}$ there exists another node $l$ such that $(i, j, l)$ form a triangle in $\mathcal{G}^{(2)}$ and $(i, l), (l, j) \in \mathcal{E}^{(1)}$, hence, $\delta(i, j) = 2$ in $\mathcal{G}^{(1)}$. Sequentially carrying this one-step computation on the whole sequence $\{ \mathcal{G}^{(r)}\}$ results in 
\begin{equation*}
    \delta(i, j) = r \mathbb{I}\{(i, j) \in \mathcal{E}^{(r)} - \mathcal{E}^{(r - 1)} \}, 
\end{equation*}
where $\mathbb{I}$ is the indicator function and $\mathcal{E}^{(0)} = \varnothing$, note that these distances correspond to $\mathcal{G}^{(1)}$, and that since all active neighborhood paths remain active if the active length increases by one (i.e., $\mathcal{A}_{r - 1} \subseteq \mathcal{A}_{r}$). By (\ref{eq: second part inverse end}) we have that 
\begin{equation} \label{eq: end first side part b}
    \rho^{(r)} \leq |[\bmat{S}(\omega)]_{ij}| < \rho^{(r - 1)} \implies (i, j) \in \mathcal{A}_r - \mathcal{A}_{r - 1} = \mathcal{P}_{r},
\end{equation}
%which shows the first direction of the second part of Theorem \ref{th: pop_sparsity gnar}.
\end{proof}
\newpage
Furthermore, this shows that defining $\rstage$ neighborhood inclusion based on 
\newline Definition \ref{def: rstage hierarchy} results in the same $\mathcal{A}_{r}$ as defining inclusion based on the underlying network $\mathcal{G}$ under the assumption that larger distance between nodes implies a weaker dependence.
\par
We proceed to prove the other direction of the second part.
\begin{proof}
Consider reparameterizing the process for some $s \leq r^*$ such that $\beta_{kr} \equiv 0$ if and only if $r > s$. 
By Proposition \ref{prop: var constrains} the new linear filters given by (\ref{eq: opt filter}) satisfy that if $(i, j) \notin \mathcal{A}_s,$ then
$$\cov (u_{i, t}^{(ij)|s}, u_{j, t + k}^{(ij)|s}) \equiv 0,$$ 
where $u_{i, t}^{(ij)|s}$ denotes the optimal conditional residuals when restricting the GNAR process to maximum active $\rstage$ depth $s$.
\par
Following the ideas above, for $r \in \{1, \cdots, r^*\}$ restricting the length of $\rstage$ neighborhood regression to a value smaller than $r^*$ results in the optimal filters
\begin{equation*}
     u_{i, t}^{(ij)|r} = X_{i, t} - a_{i} - \sum_{k = 1}^{p} \sum_{l \in \mathcal{A}_r^{(j)}} 
        \gamma_{k}^{il} X_{l, t - k},
\end{equation*}
where $\mathcal{A}_r^{(j)} = \mathcal{A}^{(ij)} \cap \mathcal{A}_r$. And, using Theorem \ref{th: inv spec cross indentity}, (\ref{eq: end first part}, \ref{eq: r stage inv matrix}) we have that
\begin{align} \label{eq: second part end equivalence neigh sets r}
    (i, j ) \in \mathcal{P}_r &\implies (i, j) \in \mathcal{A}_{r}, \nonumber \\
    &\iff |\cov (u_{i, t}^{(ij)|r}, u_{j, t + h}^{(ij)|r} )| \not\equiv 0, \nonumber \\
    &\iff |[\bmat{g}^{(r)} (\omega)]_{ij}| \not\equiv 0. 
\end{align}
\newpage
And, that the reparameterized GNAR process for distances no longer than $(r - 1)$ satisfies
\begin{align} \label{eq: second part end equivalence neigh sets r - 1}
    (i, j) \in \mathcal{P}_{r} &\implies (i, j) \notin \mathcal{A}_{r - 1} \nonumber \\
    &\iff \cov (u_{i, t}^{(ij)|r - 1}, u_{j, t + h}^{(ij)|r - 1}) \equiv 0 \nonumber \\
    &\iff [\bmat{g}^{(r - 1)} (\omega)]_{ij} \equiv 0,
\end{align}
where $\bmat{g}^{(r)} (\omega)$ and $\bmat{g}^{(r - 1)} (\omega)$ are the inverse spectral matrices of the constrained GNAR process. Also note that since
$(i, j) \in \mathcal{A}_{r}$ we have that 
$\bmat{S} (\omega) \not\equiv 0, $
which holds in conjunction with (\ref{eq: second part end equivalence neigh sets r}, \ref{eq: second part end equivalence neigh sets r - 1}). Hence, there exist two threshold values $\rho^{(r)}, \rho^{(r - 1)}$ for each $r \in \{1, \cdots, r^*\}$ in the cross-spectrum hierarchy such that
\begin{align*}
    &\bmat{S}^{(r)}(\omega) = \sft \big( \bmat{S} (\omega) + \rho^{(r)} \diag\{e^{\imath \phi_{ii}}\};  \rho^{(r)} \big) = \bmat{g}^{(r)} (\omega), \\
     &\bmat{S}^{(r - 1)}(\omega) = \sft \big( \bmat{S} (\omega) + \rho^{(r - 1)} \diag\{e^{\imath \phi_{ii}}\};  \rho^{(r - 1)} \big) = \bmat{g}^{(r - 1)} (\omega),
\end{align*}
which gives $\rho^{(r)} \leq |[\bmat{S}(\omega)]_{ij}| < \rho^{(r - 1)}$, and we have that
\begin{equation} \label{eq: end second side b}
    (i, j) \in \mathcal{P}_{r} \implies \rho^{(r)} \leq |[\bmat{S}(\omega)]_{ij}| < \rho^{(r - 1)},
\end{equation}
for all $\omega \in [-\pi, \pi]$.
\end{proof}
%Proof ends here
Finally, combining (\ref{eq: end first side part b}, \ref{eq: end second side b}) and Proposition \ref{prop: distance for conditional independence} finishes the proof.
\par
The ideas in this proof show how stronger dependence between closer nodes can be exhibited either by pruning the latent graph $\tilde{\mathcal{G}}$ or by applying the soft-threshold given by (\ref{eq: cmplx sft}) to $\bmat{S} (\omega)$ (i.e., a form of reparameterization), which are equivalent. Furthermore, they establish the equivalence between using the generalized GNAR correlation structure $\mathcal{G}$ and the notions of different degrees of dependence in Definitions \ref{def: conditional indep definition} and \ref{def: rstage hierarchy} with the notions in Section \ref{sec:gnar methodology}.
\newpage
\begin{remark} \label{rma: graph generalization}
Let $\mathcal{G}^{(1)} = (\mathcal{K}, \mathcal{E}^{(1)})$ be the graph associated to the strongest connections with edge set $\mathcal{E}^{(1)}$ given by (\ref{eq: r sequence equivalence}), then the graph $G = (\mathcal{K}, E)$ given by Definition \ref{def: conditional indep definition} is the same graph (i.e., $\mathcal{G}^{(1)} = G$).
\end{remark}

Remark \ref{rma: graph generalization} shows how GNAR processes generalize graphical multivariate time series by constraining VAR models with the introduction of higher-order interactions, furthermore, it shows that the neighborhood regression coefficients $\beta_{kr}$ exploit the correlation structure induced by $\mathcal{G}$ since they account for node-wise interactions common across the network and manage inclusion in active $\rstage$ conditional dependence paths.
\begin{corollary} \label{cor: simplify to graphical models}
If one constrains a stationary $\GNAR(p, [s_1, \cdots, s_p])$ $\boldsymbol{X}_t$ process to $1-$stage neighborhood regression, then the correlation hierarchy in Definition \ref{def: rstage hierarchy} simplifies to only distinguishing between $[\bmat{S} (\omega)]_{ij} \equiv 0$ and $[\bmat{S} (\omega)]_{ij} \not\equiv 0$, hence, Theorem \ref{th: pop_sparsity gnar} recovers the results in \cite{graph_series}.
\end{corollary}
Corollary \ref{cor: simplify to graphical models} shows that $\GNAR(p, [1, \cdots, 1])$ is equivalent to parsimonious partial correlation graphical models that satisfy the Markov properties of graphical models; see \cite{graph_series} and \cite{lauritzen}.


