Bayesian Inference for Relational Graph in a Causal Vector Autoregressive Time Series

The computation of the likelihood is difficult for a multivariate time series. For a Gaussian VAR process, the likelihood can be computed relatively fast using the Durbin-Levison (DL) or innovations algorithm. However, under different parameterizations, the computational burden can increase significantly depending on the complexity of the parameterization. Moreover, the DL-type algorithm can still be challenging if the process dimension is high.

In models with a high number of parameters, a common approach is to seek low-dimensional sub-models that could adequately describe the data and answer basic inferential questions of interest.

Using the proposed parameterization of $\bm{\Omega}$ and the differences $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$, we provide a low-rank formulation of a causal VAR process that leads to an efficient algorithm for computing the likelihood based on a Gaussian VAR($p$) sample. The algorithm achieves computational efficiency by avoiding inversion of large dimensional matrices. The sparse precision matrix $\bm{\Omega}=\bm{C}_{0}^{-1}$ is modeled as a separate parameter, allowing direct inference about the graphical structure under sparsity constraints . The likelihood computation requires the inversion of only $r_{j}\times r_{j}$ matrices instead of $d\times d$ matrices. The last fact substantially decreases the computational complexity. In particular, when $r_{j}=r=1$, i.e. when the conditional precision updates, $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$, are all rank-one, all components of likelihood computation can be done without matrix inversion or factorization.

Before we proceed, we define some notations that are used throughout the article. The stationary variance matrix $\bm{\Upsilon}_{j}$ of $(\bm{X}_{t},\bm{X}_{t-1},\ldots,\bm{X}_{t-j})^{\mathrm{T}},$ as defined in (2.3), will be written in the following nested structures

$\displaystyle\bm{\Upsilon}_{j}=\begin{bmatrix}\bm{\Gamma}(0)&\bm{\xi}_{j}^{\mathrm{T}}\\ \bm{\xi}_{j}&\bm{\Upsilon}_{j-1}\end{bmatrix}=\begin{bmatrix}\bm{\Upsilon}_{j-1}&\bm{\kappa}_{j}\\ \bm{\kappa}_{j}^{\mathrm{T}}&\bm{\Gamma}(0)\end{bmatrix}$

(3.2)

where

$\displaystyle\bm{\xi}_{j}^{\mathrm{T}}=(\bm{\Gamma}(1),\,\bm{\Gamma}(2),\ldots,\bm{\Gamma}(j)),\quad\bm{\kappa}_{j}^{\mathrm{T}}=(\bm{\Gamma}(j)^{\mathrm{T}},\ldots,\bm{\Gamma}(2)^{\mathrm{T}},\,\bm{\Gamma}(1)^{\mathrm{T}}).$

(3.3)

Denote the Schur-complements of $\bm{\Upsilon}_{j-1}$ in the two representations as

$\displaystyle\bm{C}_{j}=\bm{\Gamma}(0)-\bm{\xi}_{j}^{\mathrm{T}}\bm{\Upsilon}_{j-1}^{-1}\bm{\xi}_{j},\quad\bm{D}_{j}=\bm{\Gamma}(0)-\bm{\kappa}_{j}^{\mathrm{T}}\bm{\Upsilon}_{j-1}^{-1}\bm{\kappa}_{j}.$

(3.4)

Also, let $\phi_{d}(\cdot|\bm{\mu},\bm{\Sigma})$ and $\Phi_{d}(\cdot|\bm{\mu},\bm{\Sigma})$ denote the probability density and cumulative distribution function of the $d$-dimensional normal with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$. The basic computation will be to successively compute the likelihood contribution of the conditional densities $f(\bm{X}_{j}|\bm{X}_{j-1},\ldots,\bm{X}_{1}),\,j=2,\ldots,T$. Let $\bm{Y}_{j}=(\bm{X}_{j}^{\mathrm{T}},\ldots,\bm{X}_{\max(1,j-p)}^{\mathrm{T}})^{\mathrm{T}}$. Under the assumption that the errors are Gaussian, i.e., $\bm{Z}_{j}\sim\mathrm{N}_{d}(\bm{0},\bm{\Sigma})$, and that they are independent, and writing $f(\bm{X}_{1}|\bm{Y}_{0})=f(\bm{X}_{1})$, $\mu_{1}=\bm{0}$, $\bm{\Sigma}_{1}=\bm{C}_{0}$ we have,

$\displaystyle f(\bm{X}_{j}|\bm{Y}_{j-1})=\phi_{d}(\bm{X}_{j}|\bm{\mu}_{j},\bm{\Sigma}_{j}),$

(3.5)

where the conditional mean and variance are given by $\bm{\mu}_{j}=\bm{\xi}_{j-1}^{\mathrm{T}}\bm{\Upsilon}_{j-2}^{-1}\bm{Y}_{j-1}$, $\bm{\Sigma}_{j}=\bm{C}_{j-1}$ for any $2\leq j\leq p,$ and $\bm{\mu}_{j}=\bm{\xi}_{p}^{\mathrm{T}}\bm{\Upsilon}_{p-1}^{-1}\bm{Y}_{j}$, $\bm{\Sigma}_{j}=\bm{C}_{p}$ for $p+1\leq j\leq T$. Thus, the full Gaussian likelihood

$\displaystyle L=f(\bm{X}_{1})\prod_{j=2}^{\mathrm{T}}f(\bm{X}_{j}|\bm{Y}_{j-1})$

(3.6)

can be obtained by recursively deriving the conditional means and variances from the constraint-free parameters describing the sparse precision matrix $\bm{\Omega}$ and the reduced-rank conditional variance differences, $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1},j=1,\ldots,p$. For parameterization of the conditional precision updates, we use a low-rank parameterization. Specifically, let

$\displaystyle\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}=\bm{L}_{j}\bm{L}_{j}^{\mathrm{T}},\mbox{ where }\bm{L}_{j}\mbox{ are }d\times r_{j}\mbox{ matrices with }r_{j}\ll d.$

(3.7)

The rank factors $\bm{L}_{j}$ are not directly solvable from the precision updates $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$. To complete the parameterization, we need to define a bijection. The bijection can be defined by augmenting $\bm{L}_{j}$ with $d\times r_{j}$ semi-orthogonal matrices $\bm{Q}_{j}$ and define $\bm{L}_{j}\bm{Q}_{j}^{T}$ as a unique square-root of $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$, such as the unique symmetric square-root. This would lead to an identifiable parameterization. However, given that the $d\times r_{j}$ entries in $\bm{Q}_{j}$ are constrained by the orthogonality requirement, we will use a slightly over-identified system of parameters to describe the reduced rank formulation of $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$, $j=1,\ldots,p$. The over-identification facilitates computation enormously without creating any challenges in inference for the parameters of interest.

Specifically, for each $j=1,\ldots,p,$ along with $\bm{L}_{j}$ we will use $d\times r_{j}$ parameters arranged in a $d\times r_{j}$ matrix, $\bm{K}_{j}$, to describe the reduced-rank updates $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$. The exact definition of the matrices $\bm{K}_{j}$ along with the steps for computing the causal VAR likelihood based on the basic constraint-free parameters $\bm{\Omega}$, $\bm{L}_{j}$, and $\bm{K}_{j}$, $j=1,\ldots,p$ are given below. We assume the that the ranks $r_{1},\ldots,r_{p}$ are specified and fixed. Also, unless otherwise specified, we will assume the matrix square roots for pd matrices to be the unique symmetric square-root. Then, given the parameters $\bm{\Omega},\bm{L}_{1},\ldots,\bm{L}_{p}$ and the initialization $\bm{C}_{0}^{-1}=\bm{\Omega},$ the following steps describe the remaining parameters $\bm{K}_{1},\ldots,\bm{K}_{p}$ recursively along with recursive computation of the components $f(\bm{X}_{j}|\bm{Y}_{j-1})$ of the VAR likelihood.

For $j=1,\ldots,p,$

1. 1.

   Since $\bm{C}^{-1}_{j}-\bm{C}_{j-1}^{-1}=\bm{L}_{j}\bm{L}_{j}^{T},$ we have $\bm{C}_{j}^{-1}=\bm{C}_{j-1}^{-1}+\bm{L}_{j}\bm{L}_{j}^{\mathrm{T}}.$

2. 2.

   Using the Sherman-Woodbury-Morrison (SWM) formula for partitioned matrices, $\bm{C}_{j}=\bm{C}_{j-1}-\bm{U}_{j}\bm{U}_{j}^{\mathrm{T}},$ where

   $\displaystyle\bm{U}_{j}=\bm{C}_{j-1}\bm{L}_{j}(\bm{I}+\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j})^{-1/2}$

   (3.8)

   Note that,

   	$\displaystyle\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j}$	$\displaystyle=$	$\displaystyle(\bm{I}+\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j})^{-1/2}\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}\bm{C}_{j-1}^{-1}\bm{C}_{j-1}(\bm{I}+\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j})^{-1/2}$
   		$\displaystyle=$	$\displaystyle(\bm{I}+\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j})^{-1/2}\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j}(\bm{I}+\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j})^{-1/2}.$

   Since for a positive definite matrix $\bm{A}$, $\|(\bm{I}+\bm{A})^{-1/2}\bm{A}(\bm{I}+\bm{A})^{-1/2}\|=\frac{\|\bm{A}\|}{1+\|\bm{A}\|}<1,$ $\bm{U}_{j}$ satisfies the restriction $\|\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j}\|<1.$

3. 3.

   Define

   $\displaystyle\bm{W}_{j}=\bm{\Gamma}(j)^{\mathrm{T}}-\bm{\xi}_{j-1}\bm{\Upsilon}_{j-2}^{-1}\bm{\kappa}_{j-1},$

   (3.9)

   with $\bm{W}_{1}=\bm{\Gamma}(1)^{T}.$ From Roy et al. (2019), $\bm{U}_{j}\bm{U}_{j}^{\mathrm{T}}=\bm{W}_{j}\bm{D}_{j-1}^{-1}\bm{W}_{j}^{\mathrm{T}}.$ and hence for any $\bm{V}_{j}$ such that $\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{V}_{j}=\bm{I}$ we have

   $\displaystyle\bm{W}_{j}=\bm{U}_{j}\bm{V}_{j}^{\mathrm{T}}.$

   (3.10)

   While $\bm{U}_{j}$ are determined by $\bm{L}_{j}$, the $\bm{V}_{j}$ matrices are determined by the other parameters $\bm{K}_{j}$. We construct $\bm{V}_{j}$ from the basic parameters $\bm{K}_{j}$ as,

   $\displaystyle\bm{V}_{j}=\bm{K}_{j}(\bm{K}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{K}_{j})^{-1/2}.$

   (3.11)

   Note that $\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{V}_{j}=\bm{I}$.

4. 4.

   Thus, entries of the covariance matrices can be updated as

   $$\bm{\Gamma}(j)^{\mathrm{T}}=\bm{U}_{j}\bm{V}_{j}^{\mathrm{T}}+\bm{\xi}_{j-1}^{\mathrm{T}}\bm{\Upsilon}_{j-2}^{-1}\bm{\kappa}_{j-1},\quad\bm{\xi}_{j}^{\mathrm{T}}=(\bm{\xi}_{j-1}^{\mathrm{T}},\;\;\bm{\Gamma}(j)),\qquad\bm{\kappa}_{j}^{\mathrm{T}}=(\bm{\Gamma}(j)^{\mathrm{T}},\;\;\bm{\kappa}_{j-1}^{\mathrm{T}}).$$

5. 5.

   Applying SWM successively, we have

   	$\displaystyle\bm{D}_{j}$	$\displaystyle=$	$\displaystyle\bm{D}_{j-1}-\bm{R}_{j}\bm{R}_{j}^{\mathrm{T}},$
   	$\displaystyle\bm{D}_{j}^{-1}$	$\displaystyle=$	$\displaystyle\bm{D}_{j-1}^{-1}+\bm{P}_{j}\bm{P}_{j}^{\mathrm{T}},$

   where $\bm{R}_{j}=\bm{V}_{j}(\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j})^{1/2}$ and $\bm{P}_{j}=\bm{D}_{j-1}^{-1}\bm{R}_{j}(\bm{I}-\bm{R}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{R}_{j})^{-1/2}$. For the last update to hold one needs $(\bm{I}-\bm{R}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{R}_{j})$ to be positive definite. This follows from the fact that :

   $$\bm{R}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{R}_{j}=(\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j})^{1/2}\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{V}_{j}(\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j})^{1/2}=\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j},$$

   where the final equation holds because $\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{V}_{j}=\bm{I}.$ Recalling $\|\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j}\|<1,$ the result follows. Thus, $\bm{D}_{j}^{-1}=\bm{D}_{j-1}^{-1}+\bm{D}_{j-1}\bm{V}_{j}(\bm{I}-\bm{U}_{j}^{\mathrm{T}}\bm{C}_{j-1}^{-1}\bm{U}_{j})^{-1}\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}.$

6. 6.

   Finally, the $j$th conditional density in the likelihood is updated as

   $$f(\bm{X}_{j}\;|\;\bm{Y}_{j-1})=\phi_{d}(\bm{X}_{j}|\bm{\xi}_{j-1}^{\mathrm{T}}\bm{\Upsilon}_{j-2}^{-1}\bm{Y}_{j-1},\bm{C}_{j-1}).$$

   The determinant term can be updated recursively as $\det(\bm{C}_{j}^{-1})=\det(\bm{C}_{j-2}^{-1})\det(\bm{I}+\bm{L}_{j-1}^{\mathrm{T}}\bm{C}_{j-2}\bm{L}_{j-1}).$

The subsequent updates $j=(p+1),\dots,T$, required to compute the rest of the factors of the likelihood, can be obtained simply noting $\bm{\mu}_{j}=\bm{\xi}_{p}\bm{\Upsilon}_{p-1}^{-1}\bm{Y}_{j-1,p}$ and $\bm{C}_{j-1}^{-1}=\bm{C}_{p}^{-1}$ where $\bm{Y}_{j-1,p}=(\bm{X}_{j-1}^{T},\ldots,\bm{X}_{j-p}^{T})^{T}$

In the recursive algorithm described above, the parameters $\bm{\Omega},(\bm{L}_{1},\bm{K}_{1}),\ldots,(\bm{L}_{p},\bm{K}_{p})$ are mapped to the modified parameters $\bm{\Omega},(\bm{U}_{1},\bm{V}_{1}),\ldots,(\bm{U}_{p},\bm{V}_{p})$. The mapping is not one-to-one since each $\bm{V}_{j}$ has the restriction $\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{V}_{j}=\bm{I}$ while the basic parameters $\bm{K}_{j}$ do not have any restrictions. The intermediate parameters in the mapping, $\bm{\Omega},(\bm{U}_{1},\bm{V}_{1}),\ldots,(\bm{U}_{p},\bm{V}_{p})$, are identifiable only up to rotations. Consider the equivalence classes of pair of $d\times r_{j}$ matrices $(\bm{U},\bm{V})$ defined through the relation $(\bm{U},\bm{V})\equiv(\bm{U}\bm{Q},\bm{V}\bm{Q})$ for any $r_{j}\times r_{j}$ orthogonal matrix $\bm{Q}$. Let for each $r_{j}$, $\mathbb{C}_{r_{j}}$ be the set of such equivalence classes. Also let $\mathcal{S}^{++}_{d,s}$ be a subset of the positive definite cone $\mathcal{S}^{++}_{d}$ with specified sparsity level $s$. The following proposition shows that there is a bijection from the set $\mathcal{S}^{++}_{d,s}\times\mathbb{C}_{r_{1}}\times\cdots\times\mathbb{C}_{r_{p}}$ to the causal VAR$(p)$ parameter space defined by the VAR coefficients $\bm{A}_{1},\ldots,\bm{A}_{p}$ and innovation variance $\bm{\Sigma}$, where polynomials with coefficients $\bm{A}_{1},\ldots,\bm{A}_{p}$ satisfy low-rank restrictions and the stationary precision satisfies $s$-sparsity. The definition of sparsity is made more clear in the prior specification section. For simplicity, we use a mean zero VAR($p$) process. Before we formally state the identification result, we state and prove a short Lemma that is used to show identification.

The over-parameterization of the reduced rank matrices is intentional; it reduces computational complexity. Consistent posterior inference on the set of identifiable parameters of interest is still possible, and the details are given in subsequent sections. For a general $d\times d$ matrix with rank $r$, the number of free parameters is $d^{2}-(d-r)^{2}=2dr-r^{2}$ where $(d-r)^{2}$ is the co-dimension of the rank manifold of $d\times d$ rank $r$ matrices. In the proposed parameterization, we are parameterizing the rank $r_{j}$ updates $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$ using $2dr_{j}$ parameters in the matrix pair $(\bm{L}_{j},\bm{K}_{j}).$ Thus, we have $r_{j}^{2}$ extra parameters that are being introduced for computational convenience. Typically, $r_{j}$ will be small and hence the number of extra parameters will be small as well.

The case when the increments $\bm{C}_{j}^{-1}-\bm{C}_{j-1}^{-1}$ are rank one matrices, i.e. $r_{j}=1$ for all $j$, is of special interest. Then, the parameters $\bm{U}_{j},\bm{V}_{j}$ are fully identifiable if one fixes the sign of the first entry in $\bm{U}_{j}$. Moreover, in this case, the likelihood can be computed without having to invert any matrices or computing square roots.

Specifically, the quantities from the original algorithm can be simplified and recursively computed. The key quantities to be computed at each iteration are the matrices $\bm{C}_{j},\bm{C}_{j}^{-1},\bm{D}_{j},\bm{D}_{j}^{-1},\bm{U}_{j},\bm{V}_{j}.$ At the $j$th stage of recursion, Let $l_{j}$ and $k_{j}$ be scalar quantities defined as $l_{j}=\bm{L}^{\mathrm{T}}_{j}\bm{C}_{j-1}\bm{L}_{j}$ and $k_{j}=\bm{K}_{j}^{\mathrm{T}}\bm{D}_{j-1}^{-1}\bm{K}_{j}$. Then, given the vectors $L_{j},K_{j}$

	$\displaystyle\bm{C}_{j}^{-1}$	$\displaystyle=$	$\displaystyle\bm{C}_{j-1}^{-1}+\bm{L}_{j}\bm{L}_{j}^{\mathrm{T}},$
	$\displaystyle\bm{U}_{j}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{(1+z_{j})}}\bm{C}_{j-1}\bm{L}_{j},$
	$\displaystyle\bm{C}_{j}$	$\displaystyle=$	$\displaystyle\bm{C}_{j-1}-\bm{U}_{j}\bm{U}_{j}^{\mathrm{T}},$
	$\displaystyle\bm{V}_{j}$	$\displaystyle=$	$\displaystyle\frac{1}{k_{j}}\bm{K}_{j},$
	$\displaystyle\bm{D}_{j}$	$\displaystyle=$	$\displaystyle\bm{D}_{j-1}-\frac{l_{j}}{(1+l_{j})}\bm{V}_{j}\bm{V}_{j}^{\mathrm{T}},$
	$\displaystyle\bm{D}_{j}^{-1}$	$\displaystyle=$	$\displaystyle\bm{D}_{j-1}^{-1}+\frac{1}{1+l_{j}}\bm{D}_{j-1}\bm{V}_{j}\bm{V}_{j}^{\mathrm{T}}\bm{D}_{j-1}.$

When the dimension is large, the rank-one update parameterization can effectively capture a reasonable dependence structure while providing extremely fast updates in the likelihood computation. In the simulation section, we demonstrate the effectiveness of the rank-one update method via numerical illustrations.

If the model is assumed to be causal, the prior should charge only autoregression coefficients complying with causality restrictions while putting prior probability directly on the stationary covariance matrix. Presently, the literature lacks a prior that charges only causal processes with the prior on the precision $\bm{\Omega}$ as an independent component of the prior. Using our parametrizations, we achieve that. Recall that in our setup, $p$ is considered given; hence, no prior distribution is assigned to $p$. Also, the ranks $r_{1},\ldots,r_{p}$ are specified. While all of the procedures described in this paper go through where the $p$ ranks are potentially different from each other, we develop the methodology of low-rank updates by fixing $r_{1}=\cdots=r_{p}=r.$

We consider the modified Cholesky decomposition of $\bm{\Omega}=(\bm{I}_{p}-\bm{E})\bm{F}(\bm{I}_{p}-\bm{E})^{\mathrm{T}}$, where $\bm{E}$ is strictly lower triangular, $\bm{F}=\hbox{diag}(\bm{f})$ is diagonal with $\bm{f}$ the vector of entries. For sparse estimation of $\bm{\Sigma}^{-1}$, we put sparsity inducing prior on $\bm{E}$. We first define a hard thresholding operator $H_{\lambda}(\bm{H})=\{h_{i,j}\mathbbm{1}(|h_{i,j}|>\lambda)\}$. Then we set, $\bm{E}=H_{\lambda}(\bm{E}_{1})$, where $\bm{E}_{1}$ is a strictly lower triangular matrix. For off-diagonal entries $\bm{E}_{1}$, we let $e_{ij,1}\sim\mathrm{N}(0,\sigma^{2}_{e})$ and $\sigma^{2}_{e}\sim\hbox{Inv-Ga}(c_{1},c_{1})$ for $i<j$. The components $f_{1},\ldots,f_{p}$ of $\bm{f}$ are independently distributed according to the inverse Gaussian distributions with density function $\pi_{d}(t)\propto t^{-3/2}e^{-(t-\xi)^{2}/(2t)}$, $t>0$, for some $\xi>0$; see (Chhikara, 1988). This prior has an exponential-like tail near both zero and infinity. We put a weakly informative mean-zero normal prior with a large variance on $\xi$. While employing sparsity using hard thresholding, the modified Cholesky form is easier to work with. Lastly, we put a Uniform prior on $\lambda$.

• •

  Conditional precision updates $\bm{L}_{j}$ of dimension $d\times r$: We first build the matrix $\bm{\Lambda}=[\bm{\lambda}_{1};\bm{\lambda}_{2};\ldots\bm{\lambda}_{p}]$, where $\bm{\lambda}_{j}=\hbox{vec}(\bm{L}_{j})$ placing them one after the other and $p$ is a pre-specified maximum order of the estimated VAR model. Subsequently, we impose a cumulative shrinkage prior on $\bm{\Lambda}$ to ensure that the higher-order columns are shrunk to zero. Furthermore, for each individual $\bm{L}_{j}$ too, we impose another cumulative shrinkage prior to shrinking higher-order columns in $\bm{L}_{j}$. Our prior follows the cumulative shrinkage architecture from Bhattacharya and Dunson (2011), but with two layers of cumulative shrinkage on $\Lambda$ to shrink its higher-order columns as well as the entries corresponding to $\bm{L}_{j}$’s with higher $j$:

  $$\lambda_{\ell k}|\phi_{\ell k},\tau_{k}\psi_{k,\lceil{\ell}/{r}\rceil}\sim\mathrm{N}(0,\phi_{\ell k}^{-1}\tau_{k}^{-1}\psi_{k,\lceil\frac{\ell}{r}\rceil}^{-1}),$$

  $$\phi_{lk}\sim\mathrm{Ga}(\nu_{1},\nu_{1}),\quad\tau_{k}=\prod_{i=1}^{k}\delta_{i},\quad\psi_{k,m}=\prod_{i=1}^{m}\delta^{(k)}_{i}$$

  $$\delta_{1}\sim\mathrm{Ga}(\kappa_{1},1),\quad\delta_{i}\sim\mathrm{Ga}(\kappa_{2},1)\textrm{ for }i\geq 2.$$

  $$\delta^{(k)}_{1}\sim\mathrm{Ga}(\kappa_{1},1),\quad\delta^{(k)}_{i}\sim\mathrm{Ga}(\kappa_{2},1)\textrm{ for }i\geq 2,$$

  where $\lceil x\rceil$ stands for the ceiling function that maps to the smallest integer, greater than or equal to $x$, and Ga stands for the gamma distribution. The parameters $\phi_{lk}$ control local shrinkage of the elements in $\Lambda$, whereas $\tau_{k}$ controls column shrinkage of the $k$-$th$ column. Similarly, $\psi_{k,\lceil{\ell}/{r}\rceil}^{-1})$ helps to shrink higher-order columns in $\bm{L}_{k}$. Following the well-known shrinkage properties of multiplicative gamma prior, we let $\kappa_{1}=2.1$ and $\kappa_{2}=3.1$, which work well in all of our simulation experiments. However, in the case of rank-one updates, the $\psi_{k,m}$’s are omitted. Computationally, it seems reasonable to keep $\psi_{k,m}=1$ as long as $r$ is pre-specified to a small number.

• •

  $\bm{K}_{j}$ of dimension $d\times r$: We consider non-informative flat prior for the entries in $\bm{K}_{j}$. The rank of $\bm{K}_{j}$ is specified simultaneously with the rank of $\bm{L}_{j}$. Hence, the cumulative shrinkage prior for $\bm{K}_{j}$ is not needed.

We use Markov chain Monte Carlo algorithms for posterior inference. The individual sampling strategies are described below. Due to the non-smooth and non-linear mapping between the parameters and the likelihood, it is convenient to use the Metropolis-within-Gibbs samplers for different parameters.

Adaptive Metropolis-Hastings (M-H) moves are used to update the lower-triangular entries in the latent $\bm{E}_{1}$ and the entries in $\bm{f}$. Due to the positivity constraint on $\bm{f}$, we update this parameter in the log scale with the necessary Jacobian adjustment. Specifically, we generate the updates from a multivariate normal, where the associated covariance matrix is computed based on the generated posterior samples. Our algorithm is similar to Haario et al. (2001) with some modifications, as discussed below. The initial part of the chain relies on random-walk updates, as no information is available to compute the covariance. However, after the 3500th iteration, we start computing the covariance matrices based on the last $S$ accepted samples. The choice of 3500 is based on our extensive simulation experimentation. Instead of updating these matrices at each iteration, we update them once at the end of each 100-$th$ iteration using the last $S$ accepted samples. The value of $S$ is increased gradually. Furthermore, the constant variance, multiplied with the covariance matrix as in Haario et al. (2001) is tuned to maintain a pre-specified level of acceptance. The thresholding parameter $\lambda$ is also updated in the log scale with a Jacobian adjustment.

1. 1.

   Adaptive Metropolis-Hastings update for each column in $\bm{\Lambda}$ and full conditional Gibbs updates for the other hyperparameters.

2. 2.

   Adaptive M-H update for $\bm{K}_{j}$.

To speed up the computation, we initialize $\bm{\Omega}$ to the graphical lasso Friedman et al. (2008) output using glasso R package based on the marginal distribution of multivariate time series at every cross-section, ignoring the dependence for a warm start. From this $\bm{\Omega}$, the modified Cholesky parameters $\bm{E}$ and $\bm{F}$ are initialized. We start the chain setting $p=\min(p_{m},T/2)$, where $p_{m}$ is a pre-specified lower bound, and initialize the entries in $\bm{L}_{j}$’s and $\bm{K}_{j}$’s from $\hbox{Normal}(0,1/j)$. At the 1000th iteration, we discard the $j$’s if the entries in $\bm{L}_{j}$ have very little contribution. In the case of the M-H algorithm, the acceptance rate is maintained between 25% and 50% to ensure adequate mixing of posterior samples.

Recently, Heaps (2023) applied the parametrization from Roy et al. (2019) with flat priors on the new set of parameters and developed HMC-based MCMC computation using rstan. Our priors, however, involve several layers of sparsity structures and, hence, are inherently more complex. We primarily rely on M-H sampling, as direct computation of the gradients that are required for HMC is difficult. We leave the exploration of more efficient implementation using rstan to implement gradient-based samplers as part of possible future investigations.

For the simulation experiment, we primarily study the impact of the sample size and the sparsity level on estimation accuracy for the stationary precision matrix of a first order VAR process. We compare the accuracy of estimating the stationary precision matrix and the associated partial correlation graph for the proposed method along with two other popular graphical model estimation methods, the Gaussian Graphical Model (GGM) and the Gaussian Copula Graphical Model (GCGM).

We generate the data from a VAR(1) model with a fixed marginal precision matrix, $\bm{\Omega}_{1}$. Specifically, the data $\bm{X}_{1},\ldots,\bm{X}_{T}$ are generated from the model

$\displaystyle\bm{X}_{t}=\bm{A}_{1}\bm{X}_{t-1}+\bm{Z}_{t}$

(5.1)

with the stationary precision matrix chosen as $\bm{\Omega}=\bm{\Omega}_{1}$ and the parameters associated with the update of the conditional precision, $(\bm{L}_{1},\bm{K}_{1})$, are chosen as random $30\times 30$ with entries generated independently from $N(0,2.5^{2}).$ The VAR coefficient $\bm{A}_{1}$ and the innovation variance $\bm{\Sigma}_{1}$ are solved from the specified parameters, $\bm{\Omega}_{1},\bm{L}_{1},\bm{K}_{1},$ following the steps defined in section 3.1. The initial observation, $\bm{X}_{1}$, is generated from the stationary distribution and subsequent observations are generated via the iteration in (5.1). Two different sample sizes are used, $T=40$ and $T=60.$ For the sparse precision matrix, we used two different levels of sparsity, $15\%$ and $25\%$. The sparse precision $\bm{\Omega}_{1}$ is generated using the following method.

1. (1)

   Adjacency matrix: we generate three small-world networks, each containing 10 disjoint sets of nodes with two different choices for nei variable in sample_smallworld of igraph (Csárdi et al., 2024) as 10 and 5. Then, we randomly connect some nodes across the small worlds with probability $q$. The parameters in the small world distribution and $q$ are adjusted to attain the desired sparsity levels in the precision matrix.

2. (2)

   Precision matrix: Using the above adjacency matrix, we apply G-Wishart with scale 6 and truncate entries smaller than 1 in magnitude.

The simulated model is a full-rank 30-dimensional first-order VAR. Thus, the proposed method of fitting low-rank matrices for the conditional precision updates is merely an approximation and is fitting a possibly misspecified model. Posterior computation for the proposed method is done using the steps given in section 3.5, where a higher-order VAR with low-rank updates is used to fit the data. The upper bound on the order of the progress is chosen to be $p_{m}=10.$ Thus, we are fitting a mis-specified model, and evaluating the robustness characteristics of the method along with estimation accuracy. Table 1 shows the median MSE for estimation of the 30$\times$30 precision matrix for model (5.1). The proposed method has higher estimation accuracy than the GCGM and GGM methods that ignore the temporal dependence. Thus, explicitly modeling the dependence, even under incorrect specification of the rank, provides substantial gain in terms of MSE. The comparative performance is better for the proposed method when the sparsity level is $25\%.$ The Gaussian Graphical Model seemed sensitive to the specification of the simulation model. The MSE for he GGM increases with increasing sample size, a phenomenon presumably an artifact of ignoring the dependence in the sample. We simulated the parameters of the conditional precision update, $\bm{L}_{1},\bm{K}_{1},$ independently from $N(0,1)$ (not reported here) as well and found that the MSE for the GGM to be decreasing with the sample size. This could be because the $N(0,1)$ is more concentrated around zero, making the simulated model closer to the independent model. The MSE for all the methods goes up with a decreased sparsity level. A greater number of nonzero entries in the precision matrix makes the problem harder for sparse estimation methods and leads to higher estimation errors. Our results above are all based on rank-1 updates. We also obtain results for rank-3 updates (not shown here) but the results show marginal improvement over rank-1 updates for precision matrix estimation.