\documentclass{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}

\title{{\bf CARar GNAR comparison}}
\author{Daniel Salnikov}
\date{June 2024}

\begin{document}

\maketitle

\section{Problem setting and notation}

We are interested in modelling data that can be described as either a network time series, i.e., $\mathcal{X} = (\boldsymbol{X}_t, \mathcal{G})$, or as a spatio-temporal process characterised by a parsimonious covariance structure. In this document, we compare a conditional autoregressive temporal autoregressive (CARar) model to a global-$\alpha$ generalised network autoregressive (GNAR) model. Both models assume knowledge of an underlying adjacency structure, in the case of GNAR it is the network's adjacency matrix, whereas, for CARar it is the neighbourhood adjacency matrix for areal unit data. We will use the following notation:

\begin{itemize}
    \item $\mathcal{K} := [d] = \{1, \dots, d\}$, set of nodes, i.e., dimension of the multivariate process.
    \item $\mathcal{E} \subseteq \mathcal{K} \times \mathcal{K}$, edge set, i.e., if $(i, j) \in \mathcal{E}$, then there is an edge between nodes $i$ and $j$.
    \item $\mathcal{G} := (\mathcal{K}, \mathcal{E)}$ underlying network structure of the process, e.g., land borders between counties.
    \item $\boldsymbol{X}_t:$ a multivariate time series of dimension $d$, such that each univariate time series $X_{i, t}$ is linked to node $i$ in the network $\mathcal{G}$.
    \item $\mathbf{S}_r:$ the $r$-stage adjancency matrix of $\mathcal{G}$, i.e., $[\mathbf{S}_r]_{i j} = \mathbb{I} \{ d(i, j) = r \}$, where $\mathbb{I}$ is the indicator function and $d(i, j)$ is the shortest path between nodes $i$ and $j$ in $\mathcal{G}$.
    \item $\mathbf{W}:$ the weights matrix for network $\mathcal{G}$, i.e., $[\mathbf{W}]_{ij} = w_{ij} \in [0, 1]$ is the weight between nodes $i$ and $j$, e.g., $1 / (\sum_{\ell = 1}^d [\mathbf{S}_r]_{i \ell} )$.
    \item $\mathbf{\Sigma_u}:$ covariance matrix of the $d$-dimensional residual process $\boldsymbol{u}_t$.
    \item $\alpha_{k} \in \mathbb{R}$ autoregressive coefficient at the $k$th lag, $\beta_{kr} \in \mathbb{R}$ neighbourhood regression coefficient at network distance $r$ and $k$th lag.
\end{itemize}

\subsection{CAR model}
CARar models are hierarchical models that describe spatio-temporal interactions based on a neighbourhood structure for random effects, the model is given by

\begin{align*}
    &X_{i, t} = \mu + \phi_{i, t} + u_{i, t} \\
    &\boldsymbol{\phi}_t | \boldsymbol{\phi}_{t - 1} \sim \mathrm{N} (\rho \boldsymbol{\phi}_{t - 1}, \tau^2 \mathbf{Q}^{-1} ) \\
    &\mathbf{Q} = \xi \left \{ \mathrm{diag} \left (\sum_{j = 1}^d [\mathbf{S}_1]_{i j} \right ) - \mathbf{S}_1 \right \} + (1 - \xi ) \mathbf{I_d} \\
    &\tau^2 \sim \mathrm{inverseGamma} (a, b) \\
    &\rho, \xi \sim \mathrm{Unif} (0, 1),
\end{align*}
where $\mathbf{I_d}$ is the $d \times d$ identity matrix, $u_{i, t} \sim \mathrm{N} (0, \sigma_u^2)$ at all times $t$ is independent of $\mu$ and $\phi_{i, t}$ for all $i \in [d]$, and $a, b$ are hyperparameters. Note that given $\mu$, $\rho$, $\xi$ and $\tau^2$, we can write the model in autoregressive form as
\begin{equation} \label{eq: car model}
    \boldsymbol{X}_t = \boldsymbol{\mu} + \rho (\boldsymbol{X}_{t - 1} - \boldsymbol{\mu} ) + \boldsymbol{v}_t
\end{equation}
where $\boldsymbol{\mu} = \mu \boldsymbol{1}$, $\boldsymbol{1}$ is the $d$-dimensional vector of ones, i.e., $\boldsymbol{1} = (1, \dots, 1) \in \mathbb{R}^d$, and
\begin{equation*}
    v_{i, t} = u_{i, t} + \nu_{i, t},
\end{equation*}
where $\boldsymbol{\nu}_t \sim \mathrm{N} (\boldsymbol{0}, \tau^2 \mathbf{Q}^{-1} )$ are spatially correlated effects on the mean. Note that 
\begin{equation*}
    \boldsymbol{v}_t \sim \mathrm{N} \{ 0, (\tau^2 + \sigma_u^2) \mathbf{Q}^{-1} \},
\end{equation*}
thus, 
\begin{equation*}
    \boldsymbol{\tilde{v}}_t = \mathbf{Q}^{1 / 2} \boldsymbol{v}_t \sim \mathrm{N} \{\boldsymbol{0}, (\tau^2 + \sigma_u^2) \mathbf{I_d} \} .
\end{equation*}
\subsection{GNAR model}
The global-$\alpha$ GNAR$(p, [s_k])$ is given by
\begin{align} \label{eq: gnar model}
    \boldsymbol{Z}_{r, t} &= (\mathbf{S}_r \odot \mathbf{W} ) \boldsymbol{X}_{t} \nonumber \\
    \boldsymbol{X}_t &= \sum_{k = 1}^p \left (\alpha_k \boldsymbol{X}_{t - k} + \sum_{r = 1}^{s_k} \beta_{kr} \boldsymbol{Z}_{r, t - k} \right ) + \boldsymbol{u}_t,
\end{align}
where $\boldsymbol{u}_t \sim \mathrm{N}_d (0, \mathbf{\Sigma_u})$, $s_k$ is the maximum $r$-stage depth at lag $k$ (i.e., the number of nonzero $\beta_{kr}$ parameters at lag $k$), and $p$ is the maximum lag. 


\section{CAR as a GNAR model}
Notice that by constraining \eqref{eq: gnar model} to no neighbourhood regressions and one-lag, i.e, model order $(1, [0])$, we can estimate the spatial random effects by generalised least-squares by setting the covariance matrix equivalence $\mathbf{\Sigma_u} = (\tau^2 + \sigma_u^2) \mathbf{Q}^{-1}$, and by including an intercept term $\boldsymbol{\mu}$, we can write \eqref{eq: gnar model} as
\begin{equation*}
     \boldsymbol{X}_t = \boldsymbol{\mu} + \alpha_1 (\boldsymbol{X}_{t - 1} - \boldsymbol{\mu} ) + \boldsymbol{v}_t,
\end{equation*}
where $\boldsymbol{v}_t \sim \mathrm{N} ( 0, \mathbf{\Sigma_u})$, which is identical to the CAR model given by \eqref{eq: car model}. 
\par
Further, note that by the partitioned matrix inverse formula we have that

\begin{equation*}
    [\mathbf{Q}]_{i \cdot} = - \mathbf{Q}_{-i} [\mathbf{Q}^{- 1}]_{i \cdot},
\end{equation*}
thus, non-zero entries correspond to $1$-stage neighbours. Let $w_{i \cdot} =  [\mathbf{Q}]_{i \cdot}$, i.e., the $i$th row of $\mathbf{Q}$ without the diagonal entry. Then, by properties of the multivariate normal distribution, given all other random effects, i.e., $\phi_{j , t}$ for $j \neq i$, we can write
\begin{equation*}
    \phi_{i, t} = \rho \phi_{i, t - 1} + \beta \sum_{j = 1}^d q_{i j} (\phi_{j, t} - \rho \phi_{j, t -1} ) + u_{i, t},
\end{equation*}
where $\beta = \xi  (\tau^2 + \sigma_u^2)^{-1}$. Notice that if we set $\boldsymbol{\varphi}_t = \{ \mathbf{Q} - \mathrm{diag} (\mathbf{Q}) \} \boldsymbol{\phi}_t$, then we have that
\begin{equation*}
    \boldsymbol{\phi}_t = \rho \boldsymbol{\phi}_{t - 1} + \beta (\boldsymbol{\varphi}_t - \rho \boldsymbol{\varphi}_{t- 1} ) + \boldsymbol{u}_t,
\end{equation*}
which should be thought as conditional given all other entries for each random effect $\phi_{i, t}$. This is a global-$\alpha$ GNAR process, albeit, using conditional distributions for $1$-stage neighbourhood regressions, which might result in unstable estimators. The non-normalised weights matrix is $\{ \mathbf{Q} - \mathrm{diag} (\mathbf{Q}) \}$.


\subsection{One-step forecast}

Assume that we observe a realisation of length $T$. Then, we estimate the unknown parameters in \eqref{eq: car model} by Markov Chain Monte Carlo and compute one-step forecasts by inputting $\hat{\boldsymbol{\mu}}$ and $\hat{\rho}$, i.e, 
\begin{equation*}
    \hat{\boldsymbol{X}}_{T + 1} =  \hat{\boldsymbol{\mu}} + \hat{\rho} (\boldsymbol{X}_{T} - \hat{\boldsymbol{\mu}}),
\end{equation*}
note that we can also compute the above for each set of parameters sampled from the posterior, and subsequently compute the predictive one-step forecast, which integrates the parameters. However, computation is cumbersome and we could not find a forecast option in \texttt{CARBayesST}. 
\par
Alternatively, the {\em generalised least squares} estimator is given by 
\begin{equation*}
    \hat{\boldsymbol{\theta}} = (\mathbf{R}^T \mathbf{\Sigma}_T^{-1} \mathbf{R} )^{-1} \mathbf{R}^T  \mathbf{\Sigma}_T^{-1} \boldsymbol{y}_T,
\end{equation*}
where we have written a length $T$ realisation in linear model form, i.e., $\boldsymbol{y}_T =  \mathbf{R} \boldsymbol{\theta} + \boldsymbol{u}_T$, $\boldsymbol{u}_T \sim  \mathrm{N} (0, \mathbf{\Sigma}_T)$, $\mathbf{\Sigma}_T =  \mathbf{I_d} \otimes \mathbf{\Sigma_u}$, and $\boldsymbol{\theta} = (\mu, \rho) \in \mathbb{R}^2$ is the vector of unknown autoregressive parameters. The one-step forecast is given by 
\begin{equation*}
     \hat{\boldsymbol{X}}_{T + 1} =  \hat{\boldsymbol{\theta}}_1 + \hat{\theta}_2 (\boldsymbol{X}_{T} -  \hat{\boldsymbol{\theta}}_1),
\end{equation*}
where $ \hat{\boldsymbol{\theta}}_1 =  \hat{\theta}_1 \boldsymbol{1}$. Notice that it strongly resembles the one computed using Markov Chain Monte Carlo. Further, if we assume that $\boldsymbol{X}_t$ is a global-$\alpha$ GNAR process with (spatial)network correlated errors, then the estimator is
\begin{equation*}
     \hat{\boldsymbol{\theta}} = (\mathbf{R}^T \mathbf{\Sigma}_T^{-1} \mathbf{R} )^{-1} \mathbf{R}^T  \mathbf{\Sigma}_T^{-1} \boldsymbol{y}_T,
\end{equation*}
which is identical to the generalised least squares estimator for CARar models. Thus, one-step ahead forecasts computed by inputting $\boldsymbol{\hat{\theta}}$ are equal for GNAR$(1, [0])$ with correlated errors and CARar$(1)$ models.

\end{document}
