Estimation of High-dimensional Nonlinear Vector Autoregressive Models

Estimation of High-dimensional Nonlinear Vector Autoregressive Models

Introduction

The increasing variety of scientific applications has created a growing need for employing a large set of time series (variables) to model complex social and physical systems. This demand arises from various fields, including genomics sharon2013time, neuroscience moller2001instantaneous, pereda2005nonlinear,kato2006statistical, social networks ait2015modeling, economics barigozzi2017network, environmental studies lichstein2002spatial, and communication engineering baddour2005autoregressive. For example, economic policymakers rely on large-scale models of economic indicators sims1980, bernanke2005measuring, banbura2010large, as empirical evidence has shown that they improve forecasts and provide better estimates of how current economic shocks will propagate, which guides policy actions more effectively. Similarly, in genomics and neuroscience, the advent of high-throughput technologies has enabled researchers to collect measurements on hundreds of genes or brain regions shojaie2010discovering, seth2015granger, facilitating comprehensive modeling and deeper insights into biological mechanisms. In social sciences, many key variables are not directly observable but can be inferred through related time series variables, enabling a more nuanced understanding of policy decisions lin2020regularized. Given the wide availability of high-dimensional time series data, understanding their underlying dynamic patterns is crucial for improving practical applications in these domains.

A widely used and informative model for capturing linear temporal dependencies between time series is the vector autoregression (VAR) model. Properties of VAR have been extensively studied in low-dimensional settings; see lutkepohl2005new for a comprehensive overview. Over the past decade, a growing body of literature has leveraged structured sparsity and regularized estimation frameworks to achieve consistent estimation of VAR parameters in high-dimensional settings. basu2015 investigated the theoretical properties of Lasso-penalized high-dimensional VAR models for Gaussian processes. Their result was extended to multi-block VAR models by lin2017 and to factor-augmented VAR models by lin2020regularized. guo2016 introduced a class of VAR models with banded coefficient matrices, which was further developed into spatio-temporal VAR models by gao2019. basu2019low explored high-dimensional VAR models involving low-rank and group-sparse components in network structures. hall2018 studied regularized high-dimensional autoregressive generalized linear models, focusing on Bernoulli and Poisson distributions. Additionally, ghosh2019, ghosh2021strong developed Bayesian VAR models and analyzed their posterior and strong selection consistency. For further related work, see zheng2019testing, pandit2020generalized, wang2022high, wang2023rate, chen2023community, among others.

Although many mechanisms, such as regulatory processes in biology (cf. sima2009 for a survey), involve nonlinear dynamics, research on high-dimensional time series models addressing such dynamics remains limited. mazur2009 and aijo2009 employed Bayesian learning to manage the stochasticity of biological data. lim2015 introduced a family of VAR models using operator-valued kernels to identify nonlinear dynamic systems. zhou2018 proposed a framework for non-parametric autoregressive models within generalized linear models by utilizing reproducing kernel Hilbert spaces, analyzing the convex penalized sparse and smooth estimator. shen2019nonlinear investigated nonlinear structural VAR models with application to brain networks. Additional applications can be found in pereda2005nonlinear, balcilar2016testing, yu2021sparse, among others. Among these works, only zhou2018 provided theoretical guarantees, although their concentration inequalities are not sharp. In this paper, we extend the framework of sparse linear VAR models to sparse non-parametric nonlinear VAR models, with rigorous theoretical guarantees.

This paper has two primary objectives: (i) to develop sharp inequalities for tail probabilities for non-sub-Gaussian nonlinear VAR processes; (ii) to propose a new class of methods for high-dimensional non-parametric VAR models and to apply our inequalities to obtain theoretical properties of $\ell_1$ regularized estimators. It is expected that our framework, inequalities and tools will be useful in other high-dimensional linear and nonlinear VAR problems.

In our theoretical framework, we shall consider the following nonlinear VAR models

where $\epsilon_t\in\mathbb{R}^p, t\in\mathbb{Z}$, are i.i.d. random vectors, $X_t=(X_{t}^{(1)},\ldots,X_{t}^{(p)})^\top\in\mathbb{R}^p$, $h^{(j)}=(h_{1}^{(j)},\ldots,h_{p}^{(j)})^\top$ and $h_k^{(j)}:\mathbb{R}^p\rightarrow \mathbb{R}$, $1\le j\le d, 1\le k\le p,$ are real-valued functions. By stacking lagged vectors, we can let $d=1$ in (ref) and consider the nonlinear VAR(1) model. Then (ref) can be rewritten as

Based on model (ref), we shall develop sharp Bernstein-type inequalities. Establishing exponential-type tail probability inequalities for temporal dependent processes is a challenging problem. There has been some effort to derive concentration inequalities for non-i.i.d. processes. For example, generalizations of Bernstein's inequality to $\alpha$-mixing and $\phi$-mixing random variables have been studied in bosq1993, modha1996, samson2000 and merlevede2009, merlevede2011, among others. zhang2021robust provided Bernstein-type inequality for dependent random variables under geometric moment contraction. Exponential-type inequalities were also derived for sums of Markov chains in douc2008, adamczak2008, lemanczyk2021general. Unfortunately, all these inequalities involve extra non-constant factors to account for weak dependence, and are not as sharp as the original Bernstein's inequality for independent random variables. Recently, fan2021hoeffding and jiang2018 established sharp Hoeffding-type inequality and Bernstein-type inequality for stationary Markov dependent random variables. chen2017 derived exponential inequalities and Nagaev-type inequalities for one dimensional linear (or moving average) processes under both short- and long-range dependence. Due to the interactions between temporal and cross-sectional dependence, tail probabilities of high-dimensional time series is much more complicated than the one-dimensional ones. In this work, we establish Bernstein-type inequalities for nonlinear VAR processes. Our inequalities, up to some constants, are as sharp as the classical Bernstein inequality for i.i.d. random variables. To the best of our knowledge, we are among the first to develop such sharp Bernstein-type inequalities for time series. Notably, we do not use the commonly employed “blocking” technique for sequences of dependent random variables hall2018, which allows us to avoid logarithmic factors. Our technical approach can be used to improve existing studies on high-dimensional VAR models, such as in kock2015oracle, jiang2023autoregressive, dahlhaus2023adaptation, wang2023rate.

To study nonlinear dynamical systems from high-dimensional time series data, in this paper, we introduce sparse additive non-parametric VAR models. Our method combines ideas from sparse linear modelling, additive non-parametric regression and VAR models. Each nonlinear function $h_j$, $1\le j\le p$, in model (ref) can be expressed as: $$ h_{j}(x)=\sum_{k=1}^p h_{jk}(x_k),$$ where $x=(x_1,\ldots,x_p)^\top \in \mathbb{R}^p$ and $h_{jk}(\cdot)$ are functions of one dimensional variables. The underlying VAR model is similar to sparse linear regression, but we impose a sparsity constraint on the index set $\{(j,k): h_{jk}(\cdot)\neq 0\}$ of functions $h_{jk}$ that are not identically zero. Then we estimate each nonlinear function $h_{jk}$ in terms of a truncated set of basis functions. ravikumar2009 proposed sparse additive linear models using a basis expansion and LASSO type penalty under i.i.d. data. meier2009 considered a sparsity-smoothness penalty for high-dimensional generalized additive models. koltchinskii2010, raskutti2012minimax and yuan2016 studied a different framework, sparse additive kernel regression, for the cases where the component functions belong to a reproducing kernel Hilbert spaces (RKHS). They penalized the sum of the reproducing kernel Hilbert space norms of the component functions. Their sparse additive linear models are extended to autoregressive generalized linear models in zhou2018. lim2015 introduced operator-valued kernel-based VAR models, and developed proximal gradient descent algorithms. However, their paper does not provide any theoretical guarantees. Recently, duker2025kernel developed an RKHS-based framework for nonlinear VAR processes and derived non-asymptotic probabilistic bounds.

In this work, our method has the nice feature that it decouples smoothness and sparsity. This leads to a simple block coordinate descent algorithm (cf. ravikumar2009) that can be carried out with any non-parametric smoother and scales easily to high-dimensions. Besides, with our new probability inequalities as primary tools, we can analyze the properties of $\ell_1$ regularized estimators under non-Gaussian errors in the context where $p$ is much larger than $n$. Roughly speaking, $p$ can be as large as $e^{n^c}$ for some constant $0<c<1$ if $\epsilon_t$ has finite exponential moments, and the power constant $c$ is related to the truncated number of basis expansion. We shall give a detailed description on how the dependence measures of the processes, the moment condition of the errors, the sparsity of functions and basis expansion affect the rate of convergence and the model selection consistency of the estimator.

The rest of the paper is structured as follows. Section (ref) presents Bernstein-type inequalities for nonlinear VAR processes in (ref) under Lipschitz condition and different types of moment conditions for the error processes. In Section (ref), we first formulate an $\ell_1$ regularized optimization problem for nonlinear VAR models on the population level that induces sparsity. Then we derive a sample version of the problem using basis expansion. Theoretical properties that analyze the effectiveness of the estimators in the high-dimensional setting are also presented. Simulation studies and real data analysis are carried out in Sections (ref) and (ref), respectively. Proofs of theorems and technical lemmas are contained in Section (ref).

We now introduce some notation. For a vector $x=(x_1,\ldots,x_p)^\top$, define $\|x\|_q = (|x_1|^q+\ldots+|x_p|^q)^{1/q}$, $q\ge 1$, $\|x\|=\|x\|_2$, $\|x\|_{\infty}=\max_{1\le j\le p} |x_j|$, and $\mathrm {abs}(x):=(|x_1|,\ldots,|x_p|)^\top.$ For a matrix $A = (a_{ij})$, write $|A|_{\infty}=\max_{i,j}|a_{ij}|$, the Frobenius norm $\|A\|_F = (\sum_{ij} a_{ij}^2)^{1/2}$, the spectral norm $\|A\|_2=\max_{\|x\|_2\le 1} \|Ax\|_2$ and the matrix infinity norm $\|A\|_{\infty}=\max_{i}\sum_j |a_{ij}|$. Let $\lambda_{\min}(A)$ (resp. $\lambda_{\max}(A)$) be the minimum (resp. maximum) eigenvalue of $A$. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$, write $a_n=O(b_n)$ (resp. $a_n\asymp b_n$) if there exists a constant $C$ such that $|a_n|\leq C |b_n|$ (resp. $1/C \leq a_n/b_n\leq C$) holds for all sufficiently large $n$, and write $a_n=o(b_n)$ if $\lim_{n\to\infty} a_n/b_n =0$.

Let $\epsilon_t, t\in\mathbb Z$, be i.i.d. random vectors and $\mathcal F_k=(\ldots,\epsilon_{k-1},\epsilon_k)$. Define projection operator $P_k$, $k\in\mathbb{Z}$, by $P_k(\cdot)=\mathbb E(\cdot|\mathcal F_k)-\mathbb E(\cdot|\mathcal F_{k-1}).$ Let $(\epsilon_k')$ be an i.i.d. copy of $(\epsilon_k)$. For $X_t=g(\ldots,\epsilon_{t-1},\epsilon_t)$, where $g$ is a measurable function, we define the coupled version $X_{t,\{k\}}=g(\ldots,\epsilon_{k-1},\epsilon_{k}',\epsilon_{k+1},\ldots,\epsilon_t)$, which has the same distribution as $X_t$ with $\epsilon_{k}$ in the latter replaced by an i.i.d. copy $\epsilon_{k}'$.

Bernstein Inequalities for Nonlinear VAR Processes

Exponential inequalities play a fundamental role in high-dimensional inference. Differently from i.i.d. random variables, directly applying concentration inequalities for dependent random variables to high-dimensional time series problems may lead to suboptimal results in many cases, due to the interrelationship between temporal and cross-sectional dependencies. zhang2017,zhang2021convergence,han2023high introduced new dependence measures to describe temporal and cross-sectional dependence of high-dimensional time series, then derived Fuk-Nagaev type inequalities for heavy tailed random vectors to study statistical properties of sample mean vector, spectral density matrix estimation and robust $M$-estimation, respectively. In this section, we shall present new and powerful inequalities for tail probabilities of nonlinear vector autoregressive (VAR) processes. The processes can be non-Gaussian. In Theorem (ref), we provide Bernstein-type inequalities for nonlinear VAR process under finite moment condition and exponential moment condition, respectively. In contrast, exponential inequalities provided in basu2015 are only applicable to Gaussian processes and linear VAR models with Gaussian innovation vectors (cf. Proposition 2.4 therein).

To establish exponential inequalities, we introduce the following assumptions on the function $h$ and the errors $\epsilon_t$ in model (ref). Recall that $\| \cdot\|_\infty$ is the matrix infinity norm.

The above assumption requires componentwise Lipschitz condition for nonlinear VAR processes. This assumption can be easily extended to nonlinear VAR($d$) processes. See also chen1993, diaconis1999, jarner2001, shao2007, fan2008 and chen2016 for nonlinear autoregressive processes. Intuitively, $\rho$ quantifies the strength of dependence. For example, in one dimensional AR(1) model, $X_t=\rho X_{t-1}+\epsilon_t$. Larger $\rho$ suggests stronger dependence.

The following theorem presents a Bernstein-type inequality for bounded Lipschitz continuous functions, under both the finite moment condition and the exponential moment condition of the error vectors $\epsilon_t$, respectively.

Theorem (ref)(i) addresses the finite moment case for the error vectors $\epsilon_t$ (cf. Assumption (ref) (i)). If the error vectors $\epsilon_t, t\in\mathbb{Z},$ satisfy stronger moment condition than merely having a finite $q$-th moment, we can expect a stronger inequality than (ref). Indeed, when $\epsilon_t$ has subexponenial tail (Assumption (ref) (ii)), we obtain an improved Bernstein-type inequality in (ref). Different from Theorem (ref) (i), in Theorem (ref)(ii), function $g$ can be unbounded.

\paragraph{Proof Sketch.} The proof of Theorem (ref) is quite involved. The key steps involve employing a martingale decomposition and deriving a sharp bound for the martingale differences. To be more specific, without loss of generality, assume $\|G\|_1=1$ with $G$ defined in Assumption (ref). Recall that $\mathcal F_k=(\ldots,\epsilon_{k-1},\epsilon_k)$ and the projection operator $P_k(\cdot)=\mathbb E(\cdot|\mathcal F_k)-\mathbb E(\cdot|\mathcal F_{k-1})$, for $k\in\mathbb{Z}$. The summation can then be decomposed into a sum of martingale differences:

For $X_t=\mathcal G(\ldots,\epsilon_{t-1},\epsilon_t)$, where $\mathcal G$ is some measurable function, following wu2005nonlinear, we define the coupled version $$X_{t,\{k\}}=\mathcal G(\ldots,\epsilon_{k-1},\epsilon_{k}',\epsilon_{k+1},\ldots,\epsilon_t).$$ For $x=(x_1,\ldots,x_p)^\top$, write $\mathrm {abs}(x)=(|x_1|,\ldots,|x_p|)^\top$. Since the mapping $h$ is componentwise Lipschitz continuous, by induction, we have $\mathrm {abs}(X_t-X_{t,\{k\}})\leq H^{t-k} \mathrm {abs}(\epsilon_k-\epsilon_k')$. Hence

Since the function \( g(\cdot) \) is bounded by \( M \), it follows that \( |P_k(g(X_t))| \leq 2M \). Therefore, combining this with (ref), we obtain

Since $\|v_t\|_1$ decays exponentially fast, for all sufficiently large $t$, one shall expect the first term $v_{t-k}^\top\mathbb E(\mathrm {abs}(\epsilon_k-\epsilon_k') \big |\mathcal F_k)$ to be small. Then by carefully leveraging between the two terms as detailed in Lemma (ref), we obtain that $$\mathbb E(e^{|\xi_k|h})<\infty$$ for any $h\leq h^*$ some constant $h^*>0$. Since $\xi_k$'s are martingale differences,

where the conditional expectation in the last line can be shown to be bounded for any $h\leq h^*$ with the bound denoted by $c$. Hence

The above applies for $-n\leq k\leq n.$ For $k<-n$, we can show that those terms are negligible. The desired result then follows by Markov's inequality and recursively using (ref) for $-n\leq k\leq n.$

It should be emphasized that our Bernstein-type concentration inequalities are sharp, and does not contain any annoying extra logarithmic terms. These inequalities are useful for handling non-Gaussian VAR problems.

Sparse additive nonlinear VAR models

In this section, we study sparse additive nonlinear VAR models. We first introduce a population-level optimization problem and then derive a sample-based algorithm through basis expansion. Our theoretical analysis builds on the technical tools developed in previous section.

The model

Assume that we are provided with observed time series data $X_1,\ldots,X_n \in \mathbb{R}^p$, which are sampled from a dynamical system involving $p$ variables. Our primary goal is to infer the direct influence that each variable $j$ exerts on every other variable $k$ (with $k\neq j$, $1\le k \le p$). For instance, in the case of linear VAR models, the evolution of the system is often characterized by $X_t=GX_{t-1}+\epsilon_t$, where $G$ is a $p\times p$ coefficient matrix, and $\epsilon_t$ represents noise. In our study, we assume that a first-order stationary model provides a sufficient approximation of the temporal dependencies within the system. Accordingly, we recall the nonlinear VAR framework in (ref),

where the function $h$ can capture potentially complex, nonlinear dynamics.

In this section, we propose a new class of high-dimensional, sparse, additive non-parametric VAR models. Here, each component $h_j$ of the function $h$ is assumed to decompose additively in terms of the individual components of the state vector $x\in\mathbb{R}^p$. Specifically, we posit that for each variable $j$

where each function $h_{jk}:\mathbb{R}\rightarrow \mathbb{R}$ captures the individual contribution of the $k$-th variable to the dynamics of the $j$-th variable.

Let $\Pi$ denote the joint distribution of the vector $X_t$, and let $\Pi_k$ denote the marginal distribution of the $k$-th component $X_{t}^{(k)}$ for each $1\le k\le p$. For practical purposes, we define the $L_2(\Pi_k)$-norm of the function $h_{jk}$ as

This definition is particularly relevant because it allows us to accommodate functions $h_{jk}$ that might not be Lebesgue integrable over the entire real line; instead, the integrability is considered relative to the distribution $\Pi_k$.

The classical nonlinear ridge regression is defined as

where the norms measure the overall discrepancy and the smoothness penalty on each component. To encourage sparsity in high-dimensional settings, we replace the squared norm $\| h_{jk}\|_{\Pi_k,2}^2$ with the norm $\| h_{jk}\|_{\Pi_k,2}$ itself. This substitution leads to a population-level penalized least squares estimator defined by the optimization problem

Here, $h$ is decomposed as in (ref) into a sum of univariate functions, and $\mathcal I_{k}$ is an appropriate function class for the $k$-th component. In practice, the norm $\| h_{jk}\|_{\Pi_k,2}$ can be estimated empirically by $(n^{-1}\sum_{t=1}^nh_{jk}^2(X_{t-1}^{(k)}))^{1/2}$.

By decomposing $h_{j}$ into additive components, our framework enhances interpretability and computational efficiency in high-dimensional settings. The imposed sparsity helps to pinpoint which variables have a direct influence on the dynamics. This model is an extension of the sparse additive models developed for the i.i.d. case ravikumar2009 and is especially relevant when the system exhibits nonlinear structure that traditional linear models fail to capture, while still preserving a structure that is amenable to rigorous analysis and estimation.

For each $k\in\{1,\ldots,p\}$, let $\mathcal H_k$ denote the Hilbert subspace $L_2(\Pi_k)$ consisting of measurable functions $f(\cdot)$ satisfying $\mathbb E f(X_t^{(k)})=0$ and the norm $\|f\|= (\mathbb E f^2(X_t^{(k)}) )^{1/2} < \infty$. The inner product on $\mathcal H_k$ is defined as $$\langle f,g \rangle=\mathbb E\big( f(X_t^{(k)}) g(X_t^{(k)})\big).$$ We denote by $\mathcal H=\mathcal H_1\oplus \mathcal H_2\oplus \cdots \oplus H_p$ the Hilbert space of functions of $(x_1,\ldots,x_p)$ that admit an additive representation $m(x)=\sum_k f_k(x_k)$ with each $f_k\in\mathcal H_k, k=1,\ldots,p$.

We now impose the following assumption on our basis expansion.

For example, we can choose the Fourier basis functions to satisfy Assumption (ref). In standard nonparametric regression such as ravikumar2009, covariates are often assumed to be bounded (i.e., to have compact support). Similarly, in our nonlinear VAR framework we assume $h_{jk}(x)$ in (ref) have compact support for mathematical convenience and tractability; see also raskutti2012minimax,zhou2018. Many of our results can be extended to the case of unbounded $h_{jk}(x)$ via truncation arguments with proper tail decay conditions. We omit such arguments for the sake of presentation simplicity. For example, the Fourier basis satisfies this assumption. This assumption implies that the tail of the expansion satisfies $\sum_{l=L+1}^\infty (b_{jl}^{(l)*})^2 \le C^2L^{-2\beta}$, which corresponds to the functional class condition of ravikumar2009 and is a standard requirement in basis expansion methods. The parameter $\beta$ captures the level of smoothness, effectively linking our function class to a function space. Although one could allow $\beta$ to vary adaptively with $k$, we confine ourselves to a common smoothness level in this work.

Let $L=L_n$ be a truncation parameter, and let $h_{jk}^{(L)}$ be the approximation of $h_{jk}$ defined by

In this formulation, $h_{jk}^{(L)}$ is interpreted as the projection of $h_{jk}$ onto the truncated set of basis functions $\{\psi_{k,1}, \ldots, \psi_{k,L}\}$. Then, for $1\le j,k\le p,$ the model can be written as

is the reminder term and captures the bias introduced by truncating the basis expansion.

We now define the oracle coefficients on the population level for the basis expansion and the design matrix. For any $x=(x_1,\ldots,x_p)^\top \in \mathbb{R}^p$, set vectors

Let $r_t=(r_{t}^{(1)},\ldots,r_{t}^{(p)})^\top$. With these definitions, the model can be rewritten in a compact form as

where

Consequently, the solution to our optimization problem (ref) can be approximately estimated by solving

This formulation can be viewed as a functional version of the group lasso, and the standard convexity arguments guarantee the existence of a minimizer.

Compared with the approach in lim2015, which employs operator-valued reproducing kernels for VAR models, our formulation offers a key advantage: it decouples the smoothness and sparsity components. This separation allows us to employ a block coordinate descent algorithm (cf. ravikumar2009) to efficiently construct the estimator. In the following section, we leverage the technical tools developed in Section (ref) to establish the theoretical properties of our $\ell_1$-regularized estimator, under the assumption that the particular smoother in (ref) is used.

Asymptotic properties

To facilitate the theoretical analysis, we impose the following assumptions on the functions $h_{jk}$ ($1\le j,k\le p$) and the basis expansions. For a function $f:\mathbb{R}^d\rightarrow\mathbb{R},$ denote $\|f\|_2:=(\int_{\mathbb{R}^d} f^2(x)d x)^{1/2}$.

Condition (ref) in Assumption (ref) is similar to the smallest population eigenvalue conditions commonly used in high-dimensional statistics raskutti2011minimax,vandegeer2014. In addition, it parallels the population minimum eigenvalue condition in Assumption 4 of chen2015optimal and Assumption S.3 of belloni2019conditional for sieve basis expansion functions. If the marginal density of $X_t^{(k)}$ satisfies $0< f_{\min} \le f_k(x)\le f_{\max} <\infty$ for $1\le k\le p $ and almost all $x\in[-c_0,c_0]$, then

This verifies condition (ref) in Assumption (ref). In the following Proposition (ref), we use concentration inequalities to establish the sample version of Assumption (ref).

Assumption (ref) imposes a sparsity condition on the nonlinear functions. Structural sparsity condition is often used in high-dimensional setting, for example, cai2011 in covariance matrix estimation. To achieve convergence rates without an additional factor of $p$, as is typically desired in high-dimensional settings, global boundedness of the quantities in Assumption (ref) is usually required, as in koltchinskii2010sparsity. However, raskutti2012minimax finds an elaborate way to circumvent this requirement when studying sparse additive models with RKHS components.

The following Proposition (ref) establishes an upper bound on the remainder term $\|r_t\|_\infty$ as a function of the smoothness level $\beta$, the number of basis functions $L$, and the sparsity level $s_0$. Moreover, the quantity $\frac{1}{n}\sum_{t=1}^n [h_{jk}(X_{t-1}^{(k)})- h_{jk}^{(L)}(X_{t-1}^{(k)}) ]^2$ serves as a measure of the $L_2$ bias between $h_{jk}$ and its orthogonal projection onto the finite-dimensional subspace spanned by the chosen basis functions.

Formally, we have the following asymptotic properties for the $\ell_1$ regularized estimators. Theorem (ref) shows how the rate of convergence of $\hat b-b^*$ and the errors of the estimated functions $\hat h_{jk}$ depend on the sparsity of functions, basis expansions, the dependence strength of the processes and the moment condition.

Observe that since $s\le s_0p$, the bounds in (ref), (ref) and (ref) imply that

where $b^*$ and $b_{j}^*$ are defined in (ref) and (ref). The quantity $\rho$ measures the dependence strength of the processes, and the constant $\mu_e$ encodes the moment condition. Theorem (ref) shows that, provided Assumptions (ref) and (ref)(ii) hold with $\rho\le \rho_0<1$ and $\rho_0$ is a constant, neither the dependence strength nor the moment constant $\mu_e$ affects these convergence rates. The second terms in (ref) and (ref) quantify the bias due to truncating the basis expansion. Moreover, Theorem (ref) implies that if the noise $\epsilon_t$ has finite exponential moments, then we may allow the dimension $p$ to grow as fast as $e^{n^c}$ for some constant $0<c<1$; the exponent $c$ depends on the chosen truncation level $L$ of basis expansion.

It is instructive to compare the two terms in the tuning requirement $\lambda$ from (ref). In the case with relative low dimension $\mathrm {log}(p)\lesssim s_0^2n L^{1-2\beta}$ and low basis number $L\lesssim s_0^{2/(2\beta-1)}(n/\mathrm {log} n)^{1/(2\beta-1)}$, the basis-expansion bias term $s_0L^{1-\beta}$ dominates. On the other hand, if the dimension $p$ is large such that $\mathrm {log}(p)\gtrsim s_0^2n L^{1-2\beta}$ or basis number $L$ is large with $L\gtrsim s_0^{2/(2\beta-1)}(n/\mathrm {log} n)^{1/(2\beta-1)}$, the stochastic term $(n^{-1}L\mathrm {log}(pL))^{1/2}$ becomes the leading factor.

Next, we turn to model-selection consistency. In place of Assumptions (ref), we present an alternative condition that directly targets the support of each component. To simplify the notation, let $\Psi_{S_j}(X_{t})=(\psi_{k} (X_{t}^{(k)})^\top, k\in S_j)$ be the truncated feature vector in $\mathbb{R}^{L\cdot \textrm{Card}(S_j)}$, where $\psi_{k} $ is defined in (ref). We then assemble these vectors into the block-diagonal matrix

This assumption corresponds to the condition of ravikumar2009,ravikumar2010. Similar to Assumption (ref), (ref) and (ref) are also standard, and are commonly imposed for high-dimensional regression analysis. Besides, (ref) relates to the incoherence condition, see e.g. wainwright2009,ravikumar2010. In the following proposition, we establish a sample version of Assumption (ref).

In Theorem (ref), we show that, under certain conditions, our method recovers the sparsity pattern asymptotically. Recall $S=\{(j,k): h_{jk}(\cdot)\not\equiv 0,1\le j, k\le p\}$. Then $S=\{(j,k): b_{j,k}^*\neq 0,1\le j, k\le p\}$. Let $\hat S_n:=\{(j,k): \hat b_{j,k}\neq 0,1\le j, k\le p\}$.

In a $p$-dimensional vector time series, the pattern of direct influences among variables can be represented by a binary adjacency matrix $A = (a_{jk}) \in \{0,1\}^{p\times p}$, where \[ a_{jk} =

\] In a linear VAR model $X_t = G\,X_{t-1} + \epsilon_t$, this network structure is typically inferred from the nonzero entries of the transition matrix $G$, which is often assumed to be sparse hall2018. A theory-free principle was advocated in sims1980 for inferring economic relations between variables of linear VARs.

In our nonlinear VAR framework, each component function $h_{jk}$ quantifies the influence of $k$ on $j$. Moreover, the group lasso formulation in (ref) yields a sparse estimate $\hat b$, so that many blocks $\hat b_{j,k}$ are exactly zero. We therefore define the estimated adjacency matrix $\hat A=(\hat a_{jk})$ by \[ \hat a_{jk} =

\] Since $\hat A$ need not be symmetric, it encodes a directed graph. Our Theorem (ref) then guarantees model selection consistency for $\hat A$, ensuring that the true influence network is recovered with high probability. We demonstrate the proposed network estimation method on real data in Section (ref).

Simulation Studies

In this section, we shall evaluate the numerical performance of the proposed estimation procedures of nonlinear VAR models.

We design three different patterns of the binary transition matrix (network matrix, see Section (ref)) $A$: random, band, cluster. Typical realizations of these patterns are illustrated in Figure (ref). The pattern “cluster” has block diagonal structure, where each block is of dimension $10\times 10$ and satisfies the pattern “random”. In each dimension $j$, $1\le j\le p$, we randomly assign 5 nonzero functions, according to the pattern of the transition matrix. The relevant nonzero component functions are given by

where $\Phi(\cdot,0.5,1)$ is the Gaussian probability distribution function with mean 0.5 and standard deviation 1. In other words, for each $j$ with $1\le j\le p$, we randomly select 5 functions $h_{jk} $ ($1\le k\le p$) to be the above nonzero functions. The rest $p-5$ functions of $h_{jk}$ ($1\le k\le p$) are all zeros. Elementary calculation shows that this nonlinear VAR process is stable and satisfies Assumption (ref). In order to ensure reasonable signal to noise ratio, the error processes $\epsilon_t$ are generated from $0.2 N(0,1)$.

In all the conducted experiments, we assess the model selection performance of our procedure using the area under the receiver operating characteristic curve (AUROC) and the area under the Precison-Recall curve (AUPR) ignoring the sign (positive negative influence), where the ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) and the precision-recall curve is a plot of the precision against the recall. Define TPR, FPR, precision and recall as follows

Here TP and TN stand for true positives and true negatives, respectively, and FP and FN stand for false positives/negatives. We choose a set of data dimensions $p=20, 50,100$ while the sample size is $n= 50,100,200,500$, respectively. The empirical values reported in Tables (ref) are averages over 1000 replications.

It can be seen from Table (ref) that the proposed estimation procedure of nonlinear VAR model performs fairly well as reflected in both AUROC and AUPR. In particular, when the sample size is moderate ($n\ge 100$), our method provides pretty good AUROC in all cases. As expected, when the sample size $n$ increases, our method performs better. And both AUROC and AUPR decreases as the dimension $p$ increase. Besides, our proposed method makes no significant differences in terms of 3 patterns of transition matrix.

Real Data Analysis

We now apply our nonlinear VAR model to the analysis of a real biological gene regulatory network time series expression data. The network is an {\it E. coli} SOS DNA repair system, which has been well studied in biology, see e.g, ronen2002. The main function of the SOS signaling pathway is to regulate cellular immunity and repair DNA damage. We consider an eight gene network, part of the SOS DNA repair network in the bacteria {\it E. coli}. The time series gene expression data set of the network was collected by ronen2002. The data are kinetics of 8 genes, that is, lexA, recA, ruvA, polB, umuDC, uvrA, uvrD, uvrY, where lexA and recA are the key genes in the pathway. The 8 genes were measured at 50 instants which are evenly spaced by 6 min intervals.

We compare the performance of our method with the Lasso regularized linear VAR method (basu2015). The tuning parameter $\lambda$ in both methods and the number of basis function $L$ are chosen by time series cross-validation procedure (see han2015). Figure (ref) represents the bacterial SOS DNA repair system. Figure (ref) shows the real SOS DNA repair network, which contains 9 edges. Figures (ref) and (ref) show the inferred gene regulatory networks using our nonlinear VAR model and the $\ell_1$ regularized linear VAR model, respectively. In Figure (ref), one can see that our method finds 6 out of the 9 edges in the target network and identifies lexA as the hub gene for this network. Our method identifies most interactions except lexA$\to$ruvA, lexA$\to$uvrY and recA$\to$lexA. In comparison, in the Figure (ref), the $\ell_1$ regularized linear VAR model recognizes only 4 out of the 9 true edges, and predicts a wrong edge. Furthermore, our proposed method gives the area under ROC curve 0.8116 and the area under Precison-Recall curve 0.6836. While, the $\ell_1$ regularized linear VAR model gives AUROC 0.7222 and AUPR 0.6036. In summary, our proposed method has a better performance than the regularized linear VAR model on the SOS DNA repair network, although none of these two methods can faithfully recover all of the edges. This phenomenon also confirms that there exists nonlinear dynamics in the gene regulatory networks.

Proofs

Write random variable $\bm\xi \in {\mathcal L}^m $, $m\ge 1$, if the $m$-norm $\|\bm \xi \|_m := (\mathbb E |\bm\xi|^m)^{1/m} < \infty$. Denote $\|\bm \xi \|:=\|\bm \xi \|_2$. Let $\mathcal F_k=(\ldots,\epsilon_{k-1},\epsilon_k)$, $\mathcal F_k^n=\{\epsilon_{k},\ldots,\epsilon_n\}$, and $\mathbb E_0(X)=X-\mathbb E X.$ Define projection operator $P_k(\cdot)=\mathbb E(\cdot|\mathcal F_k)-\mathbb E(\cdot|\mathcal F_{k-1})$, $k\in\mathbb{Z}$. Let $(\epsilon_k')_{ k\in\mathbb{Z}}$ be an i.i.d. copy of $(\epsilon_k)_{k\in\mathbb{Z}}$, so that $\epsilon_i,$ $\epsilon_j'$, $i,j\in\mathbb{Z}$ are i.i.d. For any $X_t=\mathcal G(\ldots,\epsilon_{t-1},\epsilon_t)$, where $\mathcal G$ is a measurable function, we define the coupled version $X_{t,\{k\}}=\mathcal G(\ldots,\epsilon_{k-1},\epsilon_{k}',\epsilon_{k+1},\ldots,\epsilon_t)$. If $k>i$, then $X_{t,\{k\}}=X_t$.

Proofs of Theorems in Section (ref)

Proofs of Theorems in Section (ref)

By (ref), for vector $b=(b_{j,k})_{1\leq j,k\leq p}$ and $b_{j,k}\in\mathbb{R}^L$, define the $(2,\alpha)$ group structure norm

where $\alpha\ge 1$. For instance, with the choice $\alpha=1$, this norm corresponds to the regularizer that underlies the group Lasso. For $\alpha=\infty$,