PaperlayerView exact revision on arXiv

Source-linked AI summary

A Direct Estimation of High Dimensional Stationary Vector Autoregressions

Fang Han, Huanran Lu, Han Liu

arXiv:1307.0293v3stat.ML

TL;DR

High-dimensional VAR estimation is difficult when the parameter dimension exceeds the observed series length, while temporal dependence complicates theory. The paper uses a Yule–Walker-based linear program for transition-matrix estimation. It establishes estimation, sign-recovery, and prediction consistency results, with analysis that highlights the transition matrix’s spectral norm.

  • Problem

    In high-dimensional VARs, the number of parameters can exceed the time-series length, and temporal dependence adds theoretical challenges to estimation.

    Source: §1 · pp. 2, 4

  • Method

    The method uses the Yule–Walker relationship to formulate transition-matrix estimation as a linear program with an elementwise maximum-norm constraint.

    Source: §1 · p. 3

  • Results

    For lag-1 series, the paper gives explicit estimation rates under sparsity and establishes sign-recovery and prediction consistency.

    Source: §1 · p. 4

  • Takeaways & Limitations

    The analysis highlights the transition matrix’s spectral norm in estimation accuracy and does not require the irrepresentable condition for exact sign recovery.

    Source: §1 · p. 4

    One-step prediction retains an unpredictable noise term that cannot be canceled even with almost infinite data.

    Source: §4 · p. 12

Abstract

from arXiv · show

The vector autoregressive (VAR) model is a powerful tool in modeling complex time series and has been exploited in many fields. However, fitting high dimensional VAR model poses some unique challenges: On one hand, the dimensionality, caused by modeling a large number of time series and higher order autoregressive processes, is usually much higher than the time series length; On the other hand, the temporal dependence structure in the VAR model gives rise to extra theoretical challenges. In high dimensions, one popular approach is to assume the transition matrix is sparse and fit the VAR model using the "least squares" method with a lasso-type penalty. In this manuscript, we propose an alternative way in estimating the VAR model. The main idea is, via exploiting the temporal dependence structure, to formulate the estimating problem into a linear program. There is instant advantage for the proposed approach over the lasso-type estimators: The estimation equation can be decomposed into multiple sub-equations and accordingly can be efficiently solved in a parallel fashion. In addition, our method brings new theoretical insights into the VAR model analysis. So far the theoretical results developed in high dimensions (e.g., Song and Bickel (2011) and Kock and Callot (2012)) mainly pose assumptions on the design matrix of the formulated regression problems. Such conditions are indirect about the transition matrices and not transparent. In contrast, our results show that the operator norm of the transition matrices plays an important role in estimation accuracy. We provide explicit rates of convergence for both estimation and prediction. In addition, we provide thorough experiments on both synthetic and real-world equity data to show that there are empirical advantages of our method over the lasso-type estimators in both parameter estimation and forecasting.

1 Introduction

High-dimensional VAR estimation is challenging when the number of coefficients exceeds the time-series length, and temporal dependence complicates theoretical analysis. The paper proposes a Yule–Walker-based linear-programming estimator and develops estimation, sign-recovery, and prediction results.

1 Introduction

  • The paper focuses on transition-matrix estimation, which is a foundation for forecasting and also informs Granger-causality analysis.VAR analysis also targets noise covariance and data covariance, though the paper’s stated focus is transition matrices.Source: §1 · p. 2
  • High-dimensional VARs pose a dimensionality problem because the product of lag order p and series count d can exceed the time-series length T, motivating constraints such as sparsity.The introduction describes transition-matrix zeros and nonzeros as encoding Granger non-causalities and causalities.Source: §1 · pp. 2–3
  • The proposed estimator uses the relationship between transition matrices and marginal and lag-1 autocovariances to formulate estimation as a linear program, allowing weak rather than exact sparsity.The formulation adds an elementwise maximum-norm constraint for identifiability and generalizes the Dantzig selector to multivariate-response regression.Source: §1 · p. 3
  • The paper provides nonasymptotic convergence analysis highlighting the transition matrix’s spectral norm as important to estimation accuracy, without requiring the irrepresentable condition for exact sign recovery.Earlier lasso-type analyses rely on less transparent design-matrix regularity conditions.Source: §1 · p. 4
  • The estimation problem can be split into d linear programs and solved in parallel, unlike the lasso-type procedures discussed in the introduction.The method also permits transition matrices that are only weakly sparse.Source: §1 · p. 4
  • For lag-1 series, the paper gives explicit estimation rates under sparsity and establishes sign-recovery and prediction consistency results.The sign-recovery result uses maximum-norm consistency, while prediction consistency follows from L1 consistency; weak-sparsity results are deferred to Section 4.Source: §1 · p. 4

2 Background

The background specifies the stationary VAR setting, its covariance relationships, and notation used in the analysis. It relates lag-p models to lag-1 models and introduces the norms and submatrix notation used later.

2 Background

  • The background defines indexed subvectors and submatrices, along with vector and matrix norms including Lq, Lmax, and Frobenius norms.It also denotes the singular values of a matrix in descending order.Source: §2 · pp. 4–5
  • A lag-p VAR is defined through a stationary mean-zero vector sequence satisfying the autoregressive relation, with Gaussian innovations independent across time and of the past.The stated determinant condition ensures stationarity.Source: §2 · p. 5
  • For a lag-1 VAR, the Yule–Walker equation relates the transition matrix to the lag-1 and marginal covariance matrices.The paper extends this relationship to lag-p models by reformulating them as lag-1 models with augmented vectors and transition matrices.Source: §2 · pp. 5–6

3 Methods and Algorithms

The method estimates VAR transition matrices through a covariance-based optimization that uses empirical marginal and lag-1 covariances. Its subproblems can be formulated as linear programs and solved in parallel.

3 Methods and Algorithms

  • Empirical marginal and lag-1 covariances replace population covariances in the Yule-Walker-based estimate, producing a convex optimization program.The program uses the Lmax norm.Source: §3 · p. 6
  • The constraint seeks an estimate close to the true parameter, while the minimization component induces sparsity.The method uses an ℓ1 objective.Source: §3 · p. 6
  • The optimization decomposes into subproblems for each β_j, enabling parallel solution and greater computational efficiency than lasso-type procedures.Source: §3 · p. 7
  • For lag p, each transition matrix A_k is assembled from rows J_k of the estimated Ω, with J_k selecting the corresponding block.Source: §3 · p. 7
  • Positive/negative variable decomposition reformulates each optimization as a linear program that can be solved with the simplex algorithm.Source: §3 · p. 7

4 Theoretical Properties

The paper establishes estimation, sign-recovery, and prediction results for the proposed estimator, with temporal dependence affecting estimation accuracy. It also extends the results to lag-p VAR models and considers covariance structures and symmetric transition matrices.

4 Theoretical Properties

  • For lag 1, the paper presents nonasymptotic parameter-estimation convergence rates under matrix L1 and Lmax norms, assuming A1 ∈ M(q,s,Md) for some 0 ≤ q < 1.For T ≥ 6 log d + 1 and d ≥ 8, the guarantee has probability no smaller than 1 − 14d−1.Source: §4 · p. 8
  • The estimation rates depend on Σ and A_1, with ∥A_1∥_2 characterizing temporal dependence and the amount of information available for estimation.The paper attributes this dependence to the effect of ∥A_1∥_2 on sample covariance convergence.Source: §4 · pp. 8–9
  • The rate examples cover entrywise weak dependence under diagonal dominance and blockwise weak dependence under block diagonal dominance, with stated covariance conditions for simplifying the rates.The blockwise example also invokes an upper-bounded condition number.Source: §4 · pp. 9–10
  • Under exact sparsity and the corollary's conditions, thresholding recovers sign(A_1) with probability at least 1 − 14d^-1, supporting consistent estimation of the Granger-causality network.Source: §4 · p. 11
  • Prediction error includes an unpredictable new-observation noise term that cannot be canceled by more data, while the estimator-error term can converge to zero under conditions.Source: §4 · pp. 11–12
  • With a symmetric A_1, the estimator selects the smaller-magnitude entry for each symmetric pair and retains the stated L1 and L∞ error bounds.The paper then gives an L2 prediction result for this symmetric setting.Source: §4 · p. 12
  • The lag-1 results extend to lag p through the augmented VAR formulation, including analogous estimation, sign-recovery, and prediction results.Source: §4 · pp. 8, 12, 13

5 Experiments

The experiments compare the proposed estimator with ridge and lasso methods on synthetic and equity data, using cross-validation to select tuning parameters. They assess estimation and prediction accuracy, the role of transition-matrix spectral norms, computation time, and memory usage.

5.1 Cross-Validation Procedure

  • Cross-validation grids lag and tuning-parameter choices, averages prediction errors from training windows, and selects the pair with the smallest average error.When the lag is predetermined, the procedure can instead select only the tuning parameter.Source: §5.1 · p. 14

5.2.1 Performance Comparison: Lag p = 1

  • The proposed method nearly uniformly outperforms Hsu et al. (2008) and Hamilton (1994) across Frobenius, L2, and L1 estimation errors.The improvement over Hsu et al. tends to be more significant at larger dimensions; Hamilton’s method performs worse than the other two.Source: §5.2.1 · p. 15, Table 1 · p. 16

5.2.2 Synthetic Data: Lag p ≥1

  • With multiple lags, the proposed method still outperforms both competing methods in the reported synthetic experiments.The experiments use lags 1–9, dimension 50, and time-series length 100, with hub or scale-free transition patterns.Source: §5.2.2 · pp. 15, 21, Table 11 · p. 21, Table 12 · p. 22

5.2.3 Synthetic Data: Impact of Transition Matrices’ Spectral Norms

  • As the spectral norm of A1 increases, the estimation error increases, indicating that transition-matrix spectral norms matter for estimation accuracy.The experiments vary ∥A1∥2 from 0.05 to 0.9; Figure 2 plots estimation error against the spectral norm.Source: §5.2.3 · pp. 21, 23, Figure 2 · p. 22

5.2.4 Computation Time and Memory Usage

  • Under the reported setting, memory usage increases approximately linearly with the number of lags p.The experiment varies p from 1 to 9 with T = 100 and d = 50.Source: §5.2.4 · p. 23
  • In the tested setting, the proposed method is faster than ridge without parallel computation and faster than lasso with at least 16 cores.The authors attribute ridge’s slower computation to the need to calculate an inverse of a large matrix.Source: §5.2.4 · p. 23, Figure 3 · p. 24

5.3 Real Data

  • The method outperforms two competitors in minimized averaged prediction accuracy; at day t = 1258, using its cross-validation-selected p = 1 and corresponding λ, error is 7.62, versus lasso’s 11.11 and ridge’s 11.94.Source: §5.3 · p. 25

6 Discussions

The discussion presents the linear-program formulation as an alternative for VAR estimation and highlights the role of the transition-matrix operator norm in accuracy. It also marks limits of the treatment and identifies further directions.

6 Discussions

  • The proposed linear-program approach offers an alternative to penalized least squares, and the analysis highlights ∥A1∥2 as relevant to estimation accuracy.The discussion also notes that the formulation brings theoretical insights into transition-matrix estimation.Source: §6 · p. 25
  • The main text leaves estimation of Σ and Ψ untreated, although the appendix outlines possible covariance-estimation routes, including assuming Σ is sparse.The discussion says an estimator of Ψ can be obtained when Σ and A are estimated.Source: §6 · p. 26
  • For lag p, a similar estimator can be built using the augmented formulation in Equation (2.3).The equity-data table reports optimized averaged prediction errors for three methods across lags p=1 to 9, but this passage gives no values or winner.Source: §6 · p. 26, Table 14 · p. 26
  • The method is designed for stationary VAR processes; handling unstable VARs and determining lag order beyond cross-validation remain future directions.The discussion notes prior work on unstable VAR models and lag-order selection, without extending the present method to those tasks.Source: §6 · p. 26

A Proofs of Main Results

This section provides proofs of the manuscript’s main results.

A Proofs of Main Results

A.1 Proof of Theorem 4.1

The proof of Theorem 4.1 develops covariance-concentration lemmas and uses them to establish feasibility and a norm bound for the estimator. The appendix also relates sample and population covariance quantities in the max norm.

A.1 Proof of Theorem 4.1

  • The proof prepares for Theorem 4.1 with lemmas on sample covariance concentration, including an Lmax convergence bound for S to Σ.The stated bound generalizes an upper bound derived for independently generated data.Source: §A.1 · p. 27
  • A separate lemma quantifies the Lmax difference between lag-1 sample covariance S1 and population covariance Σ1.Lemma A.2 gives a probability condition when T ≥ max(6 log d + 1, 2), though the supplied passage does not include its bound.Source: §A.1 · p. 28
  • With probability at least 1 − 14d^-1, A1 is feasible for the optimization problem, yielding ∥bΩ∥1 ≤ ∥A1∥1.The proof checks feasibility and then uses the equivalence of the optimization formulations.Source: §A.1 · p. 30
  • For lag one, the proof identifies bΩ with bA1, completing the estimator correspondence used in the argument.Source: §A.1 · p. 31

A.2 Proof of the Rest Results

The remaining results follow from Theorem 4.1 or from the connection between lag-one and lag-p VAR models. The proofs also use Gaussian concentration to bound prediction error.

A.2 Proof of the Rest Results

  • Corollary 4.4 follows directly from Theorem 4.1.Source: §A.2 · p. 31
  • A prediction-error bound decomposes the error into a transition-matrix estimation term and the next noise vector’s norm.The displayed inequality bounds the prediction error by ∥A1 − bA1∥2∥XT∥2 + ∥ZT+1∥2.Source: §A.2 · p. 32
  • The proof applies a concentration inequality for Lipschitz functions of a standard Gaussian vector.The passage introduces this step for a Gaussian random vector Y ∼ Nd(0,Q).Source: §A.2 · p. 32
  • Theorem 4.7 uses the connection between lag-one and lag-p autoregressive models established in Equation (2.3).Its proof technique is similar to that of Theorem 4.1.Source: §A.2 · p. 32

B Supporting Lemmas

Appendix B presents two Gaussian concentration lemmas. Their proof sketches use Lipschitz concentration of a Gaussian norm and Gaussian tail bounds with a dimension-dependent threshold.

B Supporting Lemmas

  • Lemma B.1 considers a Gaussian vector Y ∼NT (0, Q) and states a condition η > 2T −1/2.The supplied passage does not include the lemma's resulting inequality.Source: B Supporting Lemmas · p. 33
  • The proof of Lemma B.1 applies concentration for the Lipschitz function ∥Y ∥2 and factors the difference of squared norms.It then combines that factorization with the preceding concentration result.Source: B Supporting Lemmas · p. 33
  • Lemma B.2 defines Qmax as the largest diagonal entry of Q for a Gaussian vector Z, then uses Gaussian tails with t = (Qmax · α log d)1/2.The passage states this choice yields the desired result but does not include the bound itself.Source: B Supporting Lemmas · p. 33
Loading 1307.0293v3…