Technical Summary: Bayesian Inference of Sparsity in Stable Vector Autoregressive Processes

Problem Statement

The paper addresses the challenge of performing Bayesian inference on stable graphical vector autoregressive (VAR) processes where both sparsity and stationarity are required. In high-dimensional time-series applications (e.g., genetics, neuroscience, macroeconomics), identifying directed relationships (Granger causality) and contemporaneous conditional independence is crucial. These relationships are encoded in the sparsity patterns of the autoregressive coefficient matrices ($\Phi$) and the error precision matrix ($K = \Sigma^{-1}$).

While sparsity can be induced via standard priors (e.g., spike-and-slab), enforcing stationarity is computationally difficult. A VAR process is stable if and only if the roots of its characteristic polynomial lie outside the unit circle, a condition that restricts the parameters to a complex, non-convex geometric region known as the stationary region ($C_{p,m}$).

• The Conflict: Standard sparse priors are defined on unconstrained Euclidean space. Truncating these priors to $C_{p,m}$ introduces intractable normalizing constants, making posterior inference doubly-intractable.
• The Consequence: Without explicit enforcement, posterior mass often spills outside the stationary region, particularly when the time-series length is short relative to the dimension or when sparsity is low. This leads to unbounded forecast variances and poor predictive performance.

Methodology

The authors propose a novel Bayesian framework that simultaneously enforces stationarity and induces sparsity through parameter expansion and a reparameterization strategy.

1. Reparameterization for Stationarity

To avoid the geometric complexity of $C_{p,m}$, the authors introduce an auxiliary variable $\nu \in \mathbb{R}^+$ and map the constrained autoregressive coefficients $\Phi$ to an unconstrained set of square matrices $Z = (Z_1, \dots, Z_p)$ and a scaling factor $u$.

• Mapping: The transformation is defined as $Z_s = \nu \phi_s$ and $u = \rho(C_\phi)$, where $\rho(\cdot)$ is the spectral radius and $C_\phi$ is the companion matrix.
• Inverse Mapping: The original coefficients are recovered via $\phi_s = u Z_s / \rho(C_z)^s$, where $C_z$ is the companion matrix of $Z$.
• Key Property: This bijection preserves zeros: $z_{s,ij} = 0$ if and only if $\phi_{s,ij} = 0$ (for off-diagonal elements). This allows sparsity to be induced in the unconstrained space $Z$ while guaranteeing the resulting $\Phi$ lies within the stationary region.

2. Prior Specification

• Autoregressive Coefficients ($\Phi$): A spike-and-slab prior is placed on the unconstrained matrices $Z$.
  • Off-diagonal elements $z_{s,ij}$ ($i \neq j$) are modeled as $z_{s,ij} = \gamma_{s,ij} \tilde{z}_{s,ij}$, where $\gamma_{s,ij}$ is a binary indicator and $\tilde{z}_{s,ij}$ is a continuous effect size.
  • The slab is a normal distribution, and the spike is an atom at zero (or a narrow normal).
  • Hyperparameters control the mixing probability and shrinkage.
• Error Precision Matrix ($K$): A mixture of G-Wishart distributions is used.
  • The prior is specified as $\pi(G_2, K) = \pi(G_2)\pi(K|G_2)$, where $G_2$ is the undirected graph of contemporaneous dependencies and $K|G_2 \sim W_{G_2}(d, D)$.
  • This approach, common in Gaussian graphical models, allows for sparse $K$ while maintaining conjugacy with the Gaussian likelihood, though the normalizing constant is handled via exchange algorithms.

3. Computational Inference (MCMC)

The authors develop a Metropolis-within-Gibbs algorithm to sample from the joint posterior:

1. Continuous Parameters ($\theta_{1,c}$): The vector containing the continuous parts of $Z$, the scaling factor $u$, and hyperparameters is sampled using the No-U-Turn Sampler (NUTS). This gradient-based method is essential for navigating the complex posterior dependence induced by the reparameterization.
2. Discrete Indicators ($\theta_{1,d}$): The binary indicators for sparsity in $Z$ are updated via systematic scan using their full conditional Bernoulli distributions.
3. Graph and Precision Matrix ($\theta_2$): The undirected graph $G_2$ and precision matrix $K$ are sampled using the graphical model determination algorithm (Wang and Li, 2012), which employs an exchange algorithm to bypass the intractable G-Wishart normalizing constant.

Key Results

Simulation Experiments

The authors conducted extensive simulations varying time-series length ($n$), dimension ($m$), order ($p$), and sparsity levels ($r$).

• Sparsity Identification: Enforcing stationarity significantly reduced misclassification rates for sparsity indicators, particularly when the true sparsity was low ($r=0.1$) or moderate ($r=0.5$). Without stationarity constraints, the posterior often failed to identify the correct zero patterns, leading to high Hamming distances between the true and estimated graphs.
• Stationarity Violation: In analyses where stationarity was not enforced, 10% of the simulations resulted in posterior distributions where up to 16% of the mass lay outside the stationary region, confirming the risk of unbounded variance in unconstrained models.
• Comparison: The proposed method outperformed both unconstrained Bayesian sparse VARs and Lasso-based regularized likelihood estimators, particularly in high-dimensional, low-sample regimes.

Real-World Applications

1. U.S. Macroeconomic Data:
   • Setup: A VAR(4) model on 20 variables (including GDP, CPI, and Federal Funds Rate) from 1959–2007.
   • Finding: The proposed prior (sparse + stationary) yielded the best out-of-sample predictive performance across all horizons ($h=1, 2, 4, 8$).
   • Insight: Models that did not enforce stationarity produced unbounded forecast variances. The sparse stationary model identified a strong feedback loop between the Consumer Price Index (CPI) and the Federal Funds Rate, consistent with economic theory.
2. Neuroscience (Brain Connectivity):
   • Setup: Analysis of intracranial EEG (iEEG) data from epilepsy patients, focusing on band power in beta and delta frequencies.
   • Finding: The model revealed distinct patterns of directed connectivity. In the beta band, there was a significantly higher posterior probability for left-to-right (LR) directed edges compared to other directions.
   • Insight: This directional asymmetry aligns with clinical history where the patient's seizures originated in the left hemisphere (amygdala/hippocampus) and were resolved after left temporal lobe resection. The results suggest that interictal (non-seizure) brain activity contains information about the functional network generating seizures.

Significance and Claims

The paper claims to provide the first methodology for constructing a prior for vector autoregressive coefficients that simultaneously constrains inference to the stationary region and allows for sparsity in the parameter space.

• Theoretical Contribution: The work resolves the tension between the geometric complexity of the stationary region and the need for sparse parameter estimation by introducing a zero-preserving bijection via parameter expansion.
• Practical Impact: The authors demonstrate that enforcing stationarity is not merely a theoretical constraint but a practical necessity for improving predictive accuracy and scientific interpretability, especially in high-dimensional settings where data is limited.
• Modesty: The authors note that while their approach is computationally more intensive than unconstrained methods (due to the auxiliary variable and NUTS sampling), the cost is justified by the guarantee of stationarity and the resulting improvement in inference quality. They also acknowledge that the method could be adapted to other spaces with norm-like constraints (e.g., the simplex).

The paper concludes that ignoring stationarity in graphical VARs can lead to severe inferential errors and that their proposed framework offers a robust solution for learning causal structures in stable dynamic systems.