High-dimensional Simultaneous Inference on Non-Gaussian VAR Model via De-biased Estimator

Consider a VAR(d) model:

where $X_{i}=(X_{i1},X_{i2},\dots,X_{ip})\in\mathbb{R}^{p}$ is the random process of interests, $A_{i}\in\mathbb{R}^{p\times p}$, $i=1,\dots,d$, are the transition matrices and $\varepsilon_{i}$, $i\in\mathbb{Z}$, are i.i.d. innovation vectors with zero mean and symmetric distribution, i.e. $\varepsilon_{i}=-\varepsilon_{i}$ in distribution, for all $i\in\mathbb{Z}$. By a rearrangement of variables, VAR(d) models can be formulated as VAR(1) models (see Liu and Zhang, (2021)). Therefore, without loss of generality, we shall work with VAR(1) models:

This type of random process has a wide range of application, such as finance development (Shan, (2005)), economy (Juselius, (2006)) and exchange rate dynamics (Wu and Zhou, (2010)).

To ensure model stationarity, we assume that the spectral radius $\rho(A)<1$ throughout the paper, which is also the sufficient and necessary condition for a VAR(1) model to be stationary. However, a more restrictive condition that $\|A\|<1$ is always assumed in most of the earlier work. See for example, Han et al., (2015), Loh and Wainwright, (2012) and Negahban and Wainwright, (2011). For a non-symmetric matrix $A$, it could happen that $\|A\|\geq 1$ while $\rho(A)<1$. To fill the gap between $\rho(A)$ and $\|A\|$, Basu and Michailidis, (2015) proposed stability measures for high-dimensional time series to capture temporal and cross-section dependence via the spectral density function

where $\Gamma_{X}(\ell)=\text{Cov}(X_{i},X_{i+\ell}),~{}i,\ell\in\mathbb{Z}$ is the autocovariance function of the process $\{X_{i}\}$. In a more recent work, Liu and Zhang, (2021) defined spectral decay index to connect $\rho(A)$ with $\|A\|$ from a different point of view. In this paper, we will adopt the framework of spectral decay index in Liu and Zhang, (2021).

Next, we are interested in building some estimators of $A$ for which we could establish asymptotic distribution theory. This allows one to conduct statistical inference, such as finding simultaneous confidence interval. There have been some work on the robust estimation only. Liu and Zhang, (2021) provides both a Lasso-type estimator and a Dantzig-type estimator to consistently estimate the transition coefficient $A$ given $\{X_{i}\}$, under very mild moment condition on $X_{i}$ and $\epsilon_{i}$. It turns out that both Lasso-type and Dantzig-type estimators are not unbiased for estimating the transition matrix, thus insufficient for tasks like statistical inference. Therefore, one needs to develop more refined method to establish results in terms of asymptotic distributional theory. In the following sections, we will construct a de-biased estimator based on the existing one and derive the limiting distribution for the de-biased estimator.

Note that unlike many other existing work (Han et al., (2015), Basu and Michailidis, (2015), etc.), we do not assume $\varepsilon_{i}$ to be Gaussian or sub-Gaussian. Instead, it could happen that the innovations $\varepsilon_{i}$ only have some finite moments, which makes the standard techniques for estimation and inference invalid.

In this section, we construct a de-biased estimator using the techniques introduced in Bickel, (1975). To fix the idea, let $a_{j}^{\top}$ be the $j$-th row of $A$ and $\beta^{*}=\text{Vec}(A)=(a_{1}^{\top},a_{2}^{\top},\dots,a_{p}^{\top})^{\top}\in\mathbb{R}^{p^{2}}$. Suppose we are given a consistent, possibly biased, estimator $\widehat{\beta}$ of $\beta^{*}$, i.e. $|\widehat{\beta}-\beta^{*}|=o(1)$ (for example, Lasso-type or Dantzig-type estimators in Liu and Zhang, (2021)). Define a loss function $L:\mathbb{R}^{p^{2}}\to\mathbb{R}$ as

where $\beta=(\beta_{1}^{\top},\dots,\beta_{p}^{\top})^{\top}$ with $\beta_{k}\in\mathbb{R}^{p}$ for $1\leq k\leq p$, the weight function

for some threshold $T>0$ to be determined later, and the robust loss function $\ell(x)$ satisfies:

• (i)

  $\ell(x)$ is a thrice differentiable convex and even function.

• (ii)

  For some constant $C>0$, $|\ell^{\prime}|,|\ell^{\prime\prime}|,|\ell^{(3)}|\leq C$.

We give two examples of such loss functions from Pan et al., (2021) that satisfy the above conditions.

Direct calculation shows that everywhere twice differentiable and almost everywhere thrice differentiable. Also, the derivative of first three orders are bounded in magnitude. We mention that generalization to other loss functions that does not satisfy the differentiability conditions (for example, huber loss) may be derived under more refined arguments, but will be omitted in this paper.

Denote by $\psi(x)=\ell^{\prime}(x)$ the derivative of $\ell(x)$, then $\psi(x)$ is twice differentiable by condition (i) and $|\psi(x)|\leq C$ for all $x\in\mathbb{R}$ by condition (ii). Let $\mu=(\mu_{1},\dots,\mu_{p})^{\top}\in\mathbb{R}^{p}$ with $\mu_{k}=\mathbb{E}[\psi^{\prime}(\varepsilon_{ik})]$ and $\mu^{-1}=(\mu_{1}^{-1},\dots,\mu_{p}^{-1})^{\top}$. Let $\widehat{\mu}=(\widehat{\mu}_{1},\dots,\widehat{\mu}_{p})$ be the estimate of $\mu$ with $\widehat{\mu}_{k}=\frac{1}{n}\sum_{i=1}^{n}\psi^{\prime}(\widehat{\varepsilon}_{ik})$, where $\widehat{\varepsilon}_{ik}=X_{ik}-X_{i-1}^{\top}\widehat{\beta}_{k}$. Let $\Sigma_{x}=\mathbb{E}[X_{i}X_{i}^{\top}w(X_{i})]\in\mathbb{R}^{p\times p}$ be the weighted covariance matrix and $\Omega_{x}=\Sigma_{x}^{-1}\in\mathbb{R}^{p\times p}$ be the weighted precision matrix. Denote by $\widehat{\Sigma}_{x}=n^{-1}\sum_{i=1}^{n}X_{i-1}X_{i-1}^{\top}w(X_{i-1})$ the weighted sample covariance. Furthermore, suppose that $\widehat{\Omega}_{x}$ is a suitable approximation of the weighted precision matrix $\Omega_{x}$ (e.g., CLIME estimator introduced by Cai et al., (2011)), as will be discussed in section 2.3. To ensure the validity of such estimator, the sparsity of each row of $\Omega_{x}$ is always assumed due to high dimensionality. Now we introduce a few more notations:

and analogously,

Following the one-step estimator in Bickel, (1975), we de-bias $\widehat{\beta}$ by adding an additional term involving the gradient of the loss function $L$:

To briefly explain the presence of $\widehat{\Omega}$, consider Taylor expansion of $\nabla L_{n}(\widehat{\beta})$ around $\nabla L_{n}(\beta^{*})$. Write

where the remainder term $\sqrt{n}R=o(1)$ under certain conditions. Moreover, we also hope $\Delta$ to be negligible. As will be shown in the following sections,

To this end, $\widehat{\Omega}$ needs to be a good approximation of the precision matrix $\Omega$, which inspires the construction of such $\widehat{\Omega}$. More rigorous arguments will be presented in the subsequent sections.

Note that the estimator $\check{\beta}$ is closely related to the de-sparsifying Lasso estimator (Van de Geer et al., (2014) and Zhang and Zhang, (2014)), which is employed to conduct simultaneous inference for linear regression models in Zhang and Cheng, (2017). $\check{\beta}$ will reduce to de-sparsifying Lasso estimator if the loss $\ell(x)$ in (2.4) is squared error loss and the weight $w(x)\equiv 1$. Moreover, Loh, (2018) uses this one-step estimator to build the limiting distribution of high-dimensional vector restricted to a fixed number of coordinates, and delivers a result that agrees with Bickel, (1975) for low-dimensional robust M-estimators. Different from that, we will derive such conclusions simultaneously for all $p^{2}$ coordinates of $\beta^{*}$.

In the subsequent sections, we aim at obtaining a limiting distribution for $\check{\beta}$.

In this section, we mainly discuss the validity of having $\widehat{\Omega}$ as an approximation of $\Omega$. By the structure of $\Omega$, we need to first find a suitable estimator of the weighted precision $\Omega_{x}$.

The estimation of the sparse inverse covariance matrix based on a collection of observations $\{X_{i}\}$ plays a crucial role in establishing the asymptotic distribution. In high-dimensional regime, one cannot obtain a suitable estimator for the precision matrix by simply inverting the sample covariance, as the sample covariance is not invertible when the number of features exceeds the number of observations. Depending on the purposes, various methodology have been proposed to solve problem of estimating the precision. See for example, graphical Lasso (Yuan and Lin, (2007) and Friedman et al., (2008)) and nodewise regression (Meinshausen and Bühlmann, (2006)). From a different perspective, Cai et al., (2011) proposed a CLIME approach to sparse precision estimation, which shall be applied in this paper. For completeness, we reproduce the CLIME estimator in the following.

Suppose that the sparsity of each row of $\Omega_{x}$ is at most $s$, i.e., $s=\max_{1\leq i\leq p}|\{j:\Omega_{x,ij}\neq 0\}|.$ We first obtain $\widehat{\Theta}$ by solving

for some regularization parameter $\lambda_{n}>0$. Note that the solution $\Theta$ may not symmetric. To account for symmetry, the CLIME estimator $\widehat{\Omega}_{x}$ is defined as

For more analysis of CLIME estimator, see Cai et al., (2011). Next, we present the convergence theorem for CLIME estimator.

The following theorem shows that $\|\Omega-\widehat{\Omega}\|$ enjoys the same convergence rate as in the previous theorem.

The above theorem is built upon two facts: $\widehat{\Omega}_{x}$ approximates $\Omega_{x}$ and $\widehat{\mu}$ approximates $\mu$. The result regarding the latter approximation will be given in Lemma 3.7.

In this section, consider the following hypothesis testing problem:

versus the alternative hypothesis $H_{1}:A_{ij}\neq A^{0}_{ij}$ for some $i,j=1,\dots,p$. Instead of projecting the explanatory variables onto a subspace of fixed dimension (Javanmard and Montanari, (2014), Zhang and Zhang, (2014), Van de Geer et al., (2014) and Loh, (2018)), we allow the number of testings to grow as fast as an exponential order of the sample size $n$. Zhang and Cheng, (2017) presented a more related work, where it’s also allowed that the testing size to grow as a function of $p$. However, they conducted such simultaneous inference procedure under linear regression setting with independent random variables.

Employing the de-biased estimator $\check{\beta}$ defined in (2.7), we propose to use the test statistics

where $\check{\beta}$ is defined in (2.7). In the next several theorems, we elaborate a multiplier bootstrap method to obtain the critical value of the test statistics, which requires a few scaling and moment assumptions. Recall definition 2.1 for $\tau$ and theorem 2.5 for the definition of $\gamma$. Also recall that $s=\max_{1\leq i\leq p}|\{j:\Omega_{x,ij}\neq 0\}|.$

Assumptions

• (A1)

  $\sqrt{n}T^{3}|\widehat{\beta}-\beta^{*}|_{1}^{2}=o(1)$.

• (A2)

  $\|\Omega_{x}\|_{1}^{2}s\gamma^{2}\tau^{4}T^{4}(\log p)^{3}/\sqrt{n}=o(1)$.

• (A3)

  $s\gamma\tau^{2}T^{2}(\log p)^{3/2}|\widehat{\beta}-\beta^{*}|_{1}=o(1)$.

• (A4)

  $sT^{2}(\log(pn))^{7}/n\lesssim n^{-c}$.

• (A5)

  $(\log p)^{3/2}(\log n)^{1/2}T\sqrt{s\tau}\gamma/n^{1/4}=o(1).$

Additionally, throughout the paper we assume that for some constant $C>0$, $\mathbb{E}[X_{ik}^{2}]\leq C$ and $\mathbb{E}[\varepsilon_{ik}^{2}]\leq C$ for all $1\leq k\leq p$. We also suppose that $\|\Sigma_{x}\|_{\max}=O(1)$ and $0<c\leq\lambda_{\min}(\Sigma_{x})\leq\lambda_{\max}(\Sigma_{x})\leq C$. Thus, $\|\Omega_{x}\|_{2}\leq 1/\lambda_{\min}(\Sigma_{x})=O(1)$ and $\|\Omega_{x}\|_{1}=O(\sqrt{s})$, where the row sparsity $s=\max_{1\leq i\leq p}|\{j:\Omega_{x,ij}\neq 0\}|$.

Theorem 2.8 rigorously verifies that $\sqrt{n}R=o(1)$ and $\Delta=o(1)$ in (2.2) by the proposed construction of $\widehat{\Omega}$ and suggests us to perform further analysis on $\sqrt{n}\Omega\nabla L_{n}(\beta^{*})$. To derive the limiting distribution, we shall use Gaussian approximation technique, since the classic central limit theorem fails in high-dimensional setting.

Gaussian approximation was initially invented for high-dimensional independent random variables in Chernozhukov et al., (2013) and further generalized to high-dimensional time series in Zhang and Wu, (2017). Zhang and Cheng, (2017) and Loh, (2018) applied the GA technique in Chernozhukov et al., (2013) to the derivation of asymptotic distribution in linear regression setting. However, data generated from VAR model suffers temporal dependence, which makes the aforementioned techniques unavailable. Although Zhang and Wu, (2017) established such GA results for general time series using dependence adjusted norm, direct application of their theorems does not yield desirable conclusion in ultra-high dimensional setting. This leads us to derive a new GA theorem with better convergence rate, which is achievable thanks to the structure of VAR model.

The next theorem establishes a Gaussian approximation(GA) result for the term $\sqrt{n}\Omega\nabla L_{n}(\beta^{*})$. For a more detailed description of Gaussian approximation procedure, see section XXX.

Since the covariance matrix $D$ of the Gaussian analogue $z_{i}$ is not accessible from the observation $\{X_{i}\}$, we need to give a suitable estimation of $D$ before further performing multiplier bootstrap. The next theorem delivers a consistent estimator for our purpose.

Indeed, under the scaling assumptions, $\|\widehat{D}-D\|_{\max}=o(1)$. With these preparatory results, we are ready to present the main theorem of this paper, which describes a procedure to find the critical value of $\sqrt{n}|\check{\beta}-\beta^{*}|_{\infty}$ using bootstrap.

This result suggests a way to not only find the asymptotic distribution, but also to provide an accurate critical value $c(\alpha)$ using multiplier bootstrap. Under the null hypothesis $H_{0}$, we have $\sqrt{n}|\check{\beta}-\beta^{0}|_{\infty}=\sqrt{n}|\check{\beta}-\beta^{*}|$. This verifies the validity of having (2.11) as a test statistics for simultaneous inference.