Bayesian Tensor Autoregression

• Bayesian Tensor Autoregressive Framework is a model that extends VAR to array-valued time series using tensor decomposition for dimensionality reduction.
• It leverages CP and Tucker decompositions with hierarchical shrinkage priors to address over-parameterization in complex, high-dimensional data.
• Scalable Bayesian inference methods such as block-Gibbs and adaptive rank selection enable effective forecasting in applications like international trade, neuroimaging, and network analysis.

A Bayesian Tensor Autoregressive Framework generalizes classical vector autoregressive (VAR) models to array-valued time series, addressing the over-parameterization challenges inherent to high-dimensional, multiway data such as international trade flows, neuroimaging, or multilayer networks. The approach leverages tensor decomposition—typically CANDECOMP/PARAFAC (CP) or Tucker methods—to achieve parsimony, introduces hierarchical shrinkage priors for adaptive regularization, and applies sophisticated Bayesian inference strategies to enable scalable estimation, factor identification, and uncertainty quantification in settings with thousands or millions of dynamic parameters.

1. Formal Model Specification and Tensor Decomposition

Consider a $p$-lag tensor autoregressive process for a time series $\{Y_t\}$ of $N$-order responses, each $Y_t \in \mathbb{R}^{I_1 \times \cdots \times I_N}$. A general Bayesian Tensor Autoregression (TAR) is specified as: $$Y_t = A_0 + \sum_{h=1}^{p}\mathcal{A}_h \times_{N+1} \mathrm{vec}(Y_{t-h}) + E_t$$ where:

• $\mathcal{A}_h$ is a $(N+1)$-order coefficient tensor.
• $\times_{N+1}$ denotes mode-$N+1$ contraction with lagged, vectorized data.
• $E_t$ is an $N$-way tensor-normal innovation, with

$$\mathrm{vec}(E_t) \sim \mathcal{N}\bigl(0,\,\Sigma_N \otimes \cdots \otimes \Sigma_1\bigr)$$

For practical parsimony, $\mathcal{A}_h$ is replaced by a low-rank tensor decomposition:

• CP/Canonical PARAFAC decomposition: $$\mathcal{A}_h = \sum_{r=1}^R \beta^{(r)}_1 \circ \cdots \circ \beta^{(r)}_{N+1}$$
• Tucker decomposition (for higher-order tensors or panel VARs): $$\mathcal{A} = \mathcal{G} \times_1 U^{(1)} \times_2 U^{(2)} \times_3 U^{(3)}$$, reducing parameters from $O(I^{2N})$ to $O(R \sum_{i} I_i)$ or $O(R_1 R_2 R_3 + \sum_j I_j R_j)$ (Luo et al., 2022, Fan et al., 2022, Qi, 5 Nov 2025).

By vectorizing the model, one can prove that the tensor structure induces Kronecker (mode-wise) noise covariance and allows mapping to a VAR process with block structure, facilitating analytic developments and computationally efficient inference (Billio et al., 2017).

2. Hierarchical Bayesian Prior Specification

To achieve regularization and rank/lag selection, Bayesian TAR frameworks employ structured priors:

• Multiplicative Gamma Process (MGP) prior for CP margins: For each factor,

$$\beta^{(r)}_{j,i_j} \sim \mathcal{N}\bigl(0, (\phi_{r,j,i_j}\,\tau_r)^{-1}\bigr)$$

$\phi_{r,j,i_j} \sim \mathrm{Gamma}(\nu/2, \nu/2)$; $\tau_r = \prod_{l=1}^r \delta_l$, with increasing shrinkage on higher-order components (Luo et al., 2022).

• Global-local shrinkage via hierarchical horseshoe or Laplace priors, placed both on tensor factors and core tensors in Tucker decompositions (Fan et al., 2022, Qi, 5 Nov 2025, Billio et al., 2017).
• Stick-breaking/cumulative shrinkage on lag factors, employing spike-and-slab for data-driven lag selection (Zhang et al., 2021).
• Covariance priors: Inverse-Wishart on mode marginals, optionally linked by a global hyperparameter (Billio et al., 2017, Qi, 5 Nov 2025).
• Stochastic volatility: Innovations scale factors $\omega_t = e^{h_t}$, $h_t$ AR(1), are modeled to capture time-varying volatility (Qi, 5 Nov 2025).

Priors are chosen for adaptivity, enabling the data to select effective rank ($R$), lag order, and sparsity structure, with explicit Gamma, Dirichlet, and Beta tuning parameters governing the degree of regularization.

3. Posterior Computation and Adaptive Rank Selection

Posterior inference is accomplished via block-Gibbs, Hamiltonian Monte Carlo, or Metropolis–Hastings methods, with tailored steps for core tensor components, margin factors, shrinkage hyperparameters, and covariance updates.

• Block–Gibbs on tensor factors: Reshape the likelihood as regression over each margin to enable conjugate updates; each factor vector is regressed conditional on others (exploiting multiple equivalent forms of the likelihood) (Luo et al., 2022, Fan et al., 2022).
• Adaptive rank selection: Begin with an over-fitted $R_{\text{max}}$, periodically drop or add columns based on thresholded inactivity (e.g., component mass below $1\mathrm{e}{-3}$ in $>90\%$ entries). Stop adaptation post burn-in to guarantee ergodicity (Luo et al., 2022, Zhang et al., 2021).
• Interweaving algorithms: Three-way ASIS interweaving addresses identifiability issues, improving MCMC mixing and allowing correct scale allocation to margin factors (Luo et al., 2022).
• Cumulative shrinkage and MC³: For time-varying connectivity, auxiliary-variable and Metropolis-coupled MCMC address multimodal posteriors and slow mixing (Zhang et al., 2021).
• Multimodal panel inference: For Bayesian panel tensor VARs, subject-specific random effects inherit shared tensor structure but allow cross-sectional heterogeneity without dimensional explosion (Fan et al., 2022).

Hyperparameters (e.g., adaptation rates, shrinkage thresholds) are set based on simulation calibration and mixed-model diagnostics.

4. Identification and Factor Interpretation

A major concern in tensor decompositions is non-identifiability under scale, sign, and permutation of margin factors. Bayesian TAR frameworks address this via:

• Post-processing (pivot-matching/sign-matching): After posterior sampling, match columns to a pivot sample by minimizing Euclidean distance; perform sign-flip corrections to align factor directionality (Luo et al., 2022).
• Channel separation in Tucker decomposition: Distinct margin factor matrices allow identification of key modes in economic or brain network data—e.g., import/export/goods dimensions or spatial/temporal/subject factors (Qi, 5 Nov 2025, Fan et al., 2022).
• Dynamic connectivity motifs: For neuroscientific applications, binary selector processes ($\delta_{r,t}$) allow identification of latent connectivity patterns, activating them only when supported by evidence (Zhang et al., 2021).

Interpretation of factors is context-sensitive; in trade data, margin columns correspond to major economic hubs, commodity groups, or policy blocks, while in neuroimaging data, they reveal latent functional networks and their dynamic engagement.

5. Empirical Validation and Computational Properties

Empirical studies demonstrate significant gains in forecasting efficiency, interpretability, and scalability:

• Forecast superiority: Bayesian Tensor VARs (CP or Tucker) outperform standard VARs with Minnesota, Normal-Gamma, and Horseshoe priors in MSFE, MAE, and ALPL across macroeconomic series and horizons (Luo et al., 2022, Qi, 5 Nov 2025).
• Dimensionality reduction: For $N=40$ series, tensor models infer $\sim$450 parameters versus 8,000 for dense VARs, reducing runtime from $\sim$10 hours to $\sim$3 hours (Luo et al., 2022).
• Factor alignment: Tucker factor trajectories in large trade datasets closely mirror global macro indicators (industrial production, trade volatility) and detect exogenous events (COVID-19, geopolitical conflict spikes) (Qi, 5 Nov 2025).
• Sparsity and selection: Shrinkage priors successfully drive superfluous components and lags to zero, enabling robust estimation from over-parameterized initializations (Zhang et al., 2021).
• Impulse response functions: Tensor autoregressive models generalize VAR IRF tools to shocks in multilayer networks, revealing cross-layer substitution, persistent propagation, and dependence on network topology, not size (Billio et al., 2017).

A plausible implication is that the combination of tensor decomposition, adaptive shrinkage, and stochastic volatility allows Bayesian TARs to unify dynamic factor analysis, causal connectivity, and scalable high-dimensional inference for modern panel and multivariate time series.

6. Extensions, Limitations, and Open Directions

Observed limitations and recommended practices include:

• Hyperparameter choices: Simulation studies support $\nu=3$, Gamma(5,1) for MGP, small Dirichlet concentration for global shrinkage, and empirical tuning for inactivity thresholds (Luo et al., 2022, Zhang et al., 2021).
• Ergodicity and mixing: Adaptive rank selection must be terminated before final recording to ensure valid posterior samples. Interweaving is necessary for mixing of CP margins (Luo et al., 2022).
• Post-processing necessity: Identification requires post hoc pivot/sign-matching; naively fixing factor orientations can interfere with shrinkage dynamics (Luo et al., 2022).
• Model specification: Alternatives such as Tucker decomposition can offer better forecast performance if true data-generating processes exhibit mode-wise low-rank structure (Qi, 5 Nov 2025).
• Ordering and covariance: Cholesky SV imposes variable ordering; ordering-invariant decompositions may improve robustness, albeit at higher computational cost (Luo et al., 2022).
• Further extensions: Time-varying ranks, component-wise shrinkage per margin-row, and application to higher-order arrays are areas of ongoing methodological development (Qi, 5 Nov 2025, Fan et al., 2022).

In conclusion, the Bayesian Tensor Autoregressive Framework provides a rigorously founded methodology for scalable, interpretable dynamic modeling in high-dimensional, multiway time series, validated across econometric and neuroscientific domains and offering robust handling of uncertainty, regularization, and complex dependencies.