Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,796 characters · 24 sections · 36 citation commands
Quantile Vector Autoregression without Crossing
Vector autoregressive (VAR) models are key tools in macroeconometrics for analyzing multivariate time series. Traditional VAR models focus on conditional means, which limits their ability to capture asymmetric and heterogeneous responses across the distribution of time-series outcomes. Quantile VAR (QVAR) models (e.g., ando2022quantile, ando2024scenario, chavleishvili2024forecasting) overcome this limitation by extending the idea of quantile regression (koenker1978regression) to multivariate time series and provide a more flexible framework by modeling conditional quantiles rather than conditional means.
However, QVAR models encounter two major technical challenges in applications. The primary challenge is quantile crossing, which violates the fundamental requirement that quantile functions must be non-decreasing in the quantile level. When quantile crossing occurs, the model suffers from poor interpretability, incoherent forecasts, and invalid econometric inference, thereby undermining both its practical usefulness and theoretical validity. The second challenge is the quadratic increase in the number of coefficient parameters with the number of time series, often resulting in a loss of estimation efficiency. Although this issue also arises in standard VAR models, it is more severe in QVAR because the parameters to be estimated are functions of quantiles rather than scalars.
In the literature, traditional remedies to the first challenge include post-processing rearrangement (chernozhukov2010quantile) and directly imposing monotonicity via isotonic regression (bondell2010noncrossing). However, these approaches often involve ad hoc corrections, limiting the scalability and interpretability of the resulting models. Recently, ando2025simplex proposed an alternative approach by embedding the domain of the quantile function into a simplex, which they call simplex quantile regression (SQR). The SQR method reformulates the quantile regression problem as an optimization over the simplex space, ensuring non-crossing, smoothness of the quantile function, and computational efficiency. In this paper, we incorporate the approach of ando2025simplex into the QVAR framework and propose a simplex QVAR (SQVAR) approach, which retains the same advantages of the SQR method in the context of QVAR.
Our SQVAR model serves as a fundamental tool for capturing dynamic and heterogeneous interdependence among multiple time series. However, similar to traditional VAR models, it may suffer from the "curse of dimensionality": as the number of time series and the order of AR lags increase, the variance of the estimator can become quickly large. This occurs because the number of parameters grows quadratically with the number of series under study, making statistical inference unreliable without dimension reduction. Thus, as in conventional VAR models, a key challenge in the SQVAR framework is determining which lagged variables are truly relevant for the system's dynamics.
To address this dimensionality issue, we adopt the smoothly clipped absolute deviation (SCAD) regularization method (fan2001variable), given its ability to consistently select relevant lags while allowing us to estimate active parameters without asymptotic bias. Note that, in order to recover the correct AR structure of the original QVAR model, the SCAD penalty must be properly translated to ensure an equivalent penalization structure in the context of SQVAR. Under a suitably constructed simplex embedding and the corresponding SCAD penalty, we show that our proposed estimator based on series approximation can consistently identify the correct AR structure, and enjoys consistency and asymptotic normality for the coefficients of active lags. In addition, we propose a Bayesian information criterion (BIC)-based method for selecting the penalty parameter, and establish its theoretical validity.
Similar to standard VAR models, once a QVAR model is estimated, we can perform numerical simulations using impulse response analysis to study the dynamic effects of macroeconomic or financial shocks on economic variables of interest over time. In this paper, we propose two types of impulse response analysis. The first type is the generalized impulse response analysis, which can be viewed as a QVAR version of pesaran1998generalized. In this approach, the response to a shock is measured in terms of expected outcomes. As in conventional VAR models, this approach averages out the dependence across different quantile levels and may overlook important heterogeneous aspects of economic dynamics, especially in the presence of economic asymmetries, regime changes, or tail events such as financial crises. Alternatively, ando2024scenario introduced scenario-based forecasting error variance decomposition analysis, which allows the effect of shocks to vary across different quantiles of the outcome variable's distribution. Thus, as the second approach, we extend their idea to the QVAR framework. This approach provides a more granular and informative depiction of dynamic relationships among the variables under particular economic conditions specified by the researcher.
To illustrate applications of the proposed method, we analyze daily returns of six U.S. portfolios constructed by the intersections of firm size and past returns over the period from January 1995 to June 2025. The estimation results reveal a heterogeneous dynamic structure across portfolios, providing a detailed view of their complex interactions. We then conduct a scenario-based impulse response analysis for two major events, the 2008 financial crisis and the COVID-19 pandemic, using a low-volatility market period as the baseline. This analysis shows how such shocks generate market instability through cross-portfolio dependence, demonstrating the usefulness of the proposed framework for simulating market dynamics under hypothetical scenarios of interest.
Finally, as an extension, we consider the situation where the dimension of the time series is high. As the number of time series increases, the simplex embedding collapses in the limit, and thus we need to screen relevant variables prior to model transformation. In the literature, FanLv2008 proposed the sure independence screening (SIS) methodology for linear regression, which effectively reduces dimensionality by pre-screening variables based on their marginal utilities. Building on this idea, HeWangHong2013 developed a quantile-adaptive framework that extends SIS to nonlinear variable screening in high-dimensional heterogeneous data settings. While these approaches are useful for our purpose, the response variables in our model are serially dependent, which violates the independence assumptions underlying the existing screening methods. We therefore extend their screening procedure to accommodate high dimensional time series data.
To summarize, this paper makes multiple contributions. First, we develop a simplex quantile VAR (SQVAR) framework that ensures non-crossing quantile curves in multivariate time series. Second, a SCAD-regularized estimation procedure is introduced to mitigate efficiency loss due to the over-parameterization by selecting relevant lagged variables. Third, we derive the asymptotic properties of the estimator, including convergence rates, asymptotic normality, and the consistency of lag selection. Fourth, a BIC-based method for selecting the penalty parameter is provided, and its theoretical validity is established. Fifth, the paper extends impulse response analysis to the QVAR setting, offering both generalized impulse response and scenario-based impulse response approaches. Sixth, we develop a variable screening procedure for analyzing high dimensional time series. Finally, an empirical analysis of U.S. portfolios demonstrates how the proposed framework reveals heterogeneous dynamic dependence across portfolios and enables simulation study of market dynamics.
\paragraph{Paper organization.} The rest of the paper is organized as follows. Section (ref) introduces our working model, the QVAR. In this section, we discuss how to transform a QVAR model into an SQVAR form and provide the stationarity condition of the model. In Section (ref), we describe our SCAD-penalized estimation procedure and a BIC-type model selection method. Asymptotic properties of the proposed method, including the rate of convergence, asymptotic normality, and the consistency of model selection, are presented in Section (ref). In Section (ref), we introduce two new impulse response analysis approaches. Section (ref) presents numerical studies, including a set of Monte Carlo simulations and the empirical application to daily returns of U.S. portfolios. Section (ref) discusses a variable pre-screening procedure for high-dimensional time series. Section (ref) concludes. All technical proofs are relegated to the appendix.
\paragraph{Notation.} For a positive integer $a$, we write $[a] \coloneqq \{1, 2, \ldots, a\}$. $\sigma(X)$ denotes the $\sigma$-field generated by the random variable $X$. For a $p \times q$ matrix $M \coloneqq (m_{ij})$, define the norms $||M|| \coloneqq \left(\sum_{i=1}^p\sum_{j=1}^q m_{ij}^2\right)^{1/2}$, $||M||_2 \coloneqq \left(\overline{\text{eig}}(M^\top M)\right)^{1/2}$, and $||M||_\infty \coloneqq \max_{1 \leq i \leq p}\sum_{j=1}^q |m_{ij}|$, where $M^\top$ is the transpose of $M$, and, for a nonnegative definite matrix $A$, $\overline{\text{eig}}(A)$ and $\underbar{\text{eig}}(A)$ denote its maximum and minimum eigenvalues, respectively. For a function $f$ on $[0,1]$, we denote $||f||_2 \coloneqq \left(\int_0^1 |f(u)|^2 \mathrm{d}u \right)^{1/2}$, and $||f||_\infty \coloneqq \sup_{u \in [0,1]}|f(u)|$. For a random variable $X$, define $||X||_q \coloneqq \left( \mathbb{E}|X|^q \right)^{1/q}$. Finally, $c$, $\overline c$, and $\underline{c}$ possibly with subscripts denote generic positive constants.
Consider an $n$-dimensional vector time series $Y_t = (y_{1t}, \ldots ,y_{nt})^\top$ with $T$ equally-spaced observations $\{Y_1, \ldots , Y_T\}$. For these $n$ time series, we assume that they follow a quantile vector autoregressive (QVAR) model of order $p$, such that their conditional quantile functions are given by
for $\tau \in (0,1)$, where $\mathcal F_t = \sigma\{Y_s : s \le t\}$ is the information set available at time $t$, and $\theta_{0i}$ and $\theta_{li}^{(j)}$ denote unknown functions to be estimated, for $l,i \in [n]$ and $j \in [p]$.
Similar to koenker2006quantile, given the previous outcomes $Y_{t-1}, Y_{t-2}, \ldots$, we can characterize the data generating process (DGP) of $Y_t$ at each $t$ via a random-coefficient specification by introducing an $n$-dimensional uniform rank variable $\bm{U}_t = (U_{1t}, \ldots, U_{nt})$ in place of $\tau$ in (ref). That is, the DGP for the $i$-th time series $y_{it}$ can be expressed as
where $\mu_i \coloneqq \mathbb E [\theta_{0i}(\text{Uniform}[0,1])]$ and $\varepsilon_{it} \coloneqq \theta_{0i}(U_{it}) - \mu_i$. Here, $\varepsilon_{it}$ is interpreted as the innovation error term. The randomness of each $y_{it}$ is due to the randomness from all of $\bm U_t$, $\bm U_{t-1}$, $\ldots$. For the stochastic process $\bm U_t$, we impose the following assumption.
Assumption (ref) allows for arbitrary cross-sectional dependence across the elements of $\bm U_t$. Thus, as in the traditional VAR model, correlation across the innovation terms $(\varepsilon_{1t}, \ldots, \varepsilon_{nt})$ is permitted through the correlation of $U_{it}$'s.
As another important feature of our model, we allow the model (ref) to be potentially over-parameterized in that the following two sets may be nonempty:
where $\mathcal S \coloneqq [n] \times [p]$. This sparse structure is important in practice because researchers usually have no prior information about the lag order $p$, and hence may wish to adopt a relatively large $p$ to avoid misspecification. However, as $p$ increases, the number of coefficients to be estimated grows rapidly, causing a severe finite-sample efficiency issue. Note that this issue is more serious in QVAR than in the standard VAR, in which the parameters to be estimated are scalars, while they are functional parameters in QVAR. To mitigate this problem, we advocate the use of a SCAD-penalized estimator that can automatically detect $\mathcal{S}_{0i}$ and $\mathcal{S}_{1i}$. We also discuss in Section (ref) how to roughly trim irrelevant lags in advance, particularly when the number of time series $n$ is large, as an extension of the present case.
The quantile functions in (ref) are monotonically increasing by definition. To effectively incorporate the monotonicity into estimation procedure, we introduce a notion similar to the maximum effective region (MER) proposed in ando2025simplex: the set of past outcome values where the model satisfies the quantile monotonicity.
Define $N \coloneqq np$, $\Theta_i^{(j)}(\tau) = (\theta_{1i}^{(j)}(\tau), \ldots, \theta_{ni}^{(j)}(\tau))^\top$, $\underbracket{\bm{\Theta}_i(\tau)}_{(N + 1) \times 1} = (\theta_{0i}(\tau), \Theta_i^{(1)}(\tau)^\top, \ldots, \Theta_i^{(p)}(\tau)^\top)^\top$, and $\underbracket{W_t}_{(N + 1) \times 1} = (1, Y_{t - 1}^\top, \ldots, Y_{t - p}^\top)^\top$. Then, for each $i \in [n]$, we can succinctly write (ref) as
Next, we define the strict MER (SMER) as follows.
The SMER is a stronger version of the MER in ando2025simplex by requiring that the conditional quantile function of $y_{it}$ given $W_t = w$ is strictly increasing in $\tau$ for any $w \in \text{SMER}_i$. When $\text{SMER}_i$ is required to coincide with the entire support of $W_t$, this holds if the conditional cumulative distribution function (CDF) of $y_{it}$ given $W_t = w$ is strictly increasing for any $w$ in the support of $W_t$, which is a common requirement in the quantile regression literature. However, as shown below, our estimation procedure requires further that certain points potentially outside the support of $W_t$ also be included in $\text{SMER}_i$. Thus, to guarantee this property, we impose some structural assumptions directly on the form of the true coefficient functions (see Assumption (ref)(ii) below). This assumption is a cost of our SQVAR method, incurred in exchange for its nice properties, and is usually not required in other monotone quantile regression approaches.
Following ando2025simplex, we transform the QVAR model into a simplex QVAR (SQVAR) model with monotonic coefficient functions. Let $\bm v_0, \bm v_1, \ldots, \bm v_N \in \mathbb{R}^{N+1}$ be affine independent reference vertices such that, for any $w$ in the support of $W_t$, there exist nonnegative weights $\{c_j(w)\}_{j=0}^N$ with $\sum_{j=0}^N c_j(w)=1$ satisfying
In other words, $\{c_0(w), c_1(w), \ldots, c_N(w)\}$ represent the barycentric coordinates of $w$ with respect to the vertices $\bm v_0, \bm v_1, \ldots, \bm v_N$. Once such a coordinate system is specified (an example will be given later), the model can be rewritten as
where $C_t = (c_0(W_t), \ldots, c_N(W_t))^\top$, $\bm \Phi_i(U_{it}) = (\phi_{0i}(U_{it}), \ldots, \phi_{Ni}(U_{it}))^\top$, and $\phi_{ji}(U_{it}) \coloneqq \bm v_j^\top \bm \Theta_i(U_{it})$. By construction, if $\bm v_j \in \text{SMER}_i$, then the SQVAR coefficient function $u \mapsto \phi_{ji}(u)$ is strictly increasing. Therefore, if this property holds for all $j$, the nonnegativity of $C_t$ ensures that estimating the $\phi_{ji}$'s under monotonicity constraints automatically yields a monotone conditional quantile curve estimate, as desired. Note also that, for each $w$ to be embedded in the barycentric coordinate system, the support of $W_t$ must be bounded. We formally state these requirements in the following assumption.
Assumption (ref)(i) requires that the observations of $y_{it}$ are bounded within a fixed range. As discussed in bondell2010noncrossing, in the context of linear quantile regression, estimating non-crossing quantile curves on an unbounded domain automatically requires that the conditional quantile function has constant slopes (i.e., a location-shift model). Thus, the boundedness condition in Assumption (ref)(i) allows the coefficient functions to vary with $\tau$. The bound $[\text{lb}_i, \text{ub}_i]$ does not need to be the tightest possible range for $y_{it}$, but it is assumed to be known. In practice, one can simply use the empirical minimum and maximum of $y_{it}$ for $\text{lb}_i$ and $\text{ub}_i$, respectively. Under Assumption (ref), we obtain the following result, which parallels Theorem 3.1 of ando2025simplex.
The last part of Lemma (ref) implies that
must hold, where $\widetilde V_l^{(j)}$ denotes the corresponding $(l,j)$-th row of $\bm V^{-1}$. Note also that since $W_t ^\top \bm \Theta_i(U_{it}) = W_t ^\top \bm V^{-1}\bm \Phi_i(U_{it}) = C_t^\top \bm \Phi_i(U_{it})$, implying $(\bm V^{-1})^{\top} W_t = C_t$; thus, each element of $C_t$ is represented as a linear combination of $(1, y_{1,t-1}, \ldots, y_{n,t-p})$.
Hereinafter, to be consistent with the indexing in the original QVAR model and to improve readability, we shall also use the notations $\bm{\Phi}_i = (\phi_{0i},\phi_{1i}^{(1)}, \ldots, \phi_{ni}^{(1)}, \ldots, \phi_{1i}^{(p)}, \ldots, \phi_{ni}^{(p)})^\top$, $\bm V = (\bm v_0, \bm v_1^{(1)}, \ldots, \bm v_n^{(1)}, \ldots, \bm v_1^{(p)}, \ldots, \bm v_n^{(p)})^\top$, and $C_t = (c_{0t}, c_{1,t-1}, \ldots, c_{n,t-1}, \ldots, c_{1,t-p}, \ldots, c_{n,t-p})^\top$.
Stationarity is a fundamental concept in the analysis of VAR models. Whether the time-series process is stationary or not determines how we estimate and interpret the QVAR model. In the following, we derive the condition for our QVAR model to generate stationary processes.
Let $A_j(\bm U_t) \coloneqq (\theta^{(j)}_{uv}(U_{vt}))_{u,v=1}^n$ be an $n\times n$ matrix and $E(\bm{U}_t) \coloneqq (\theta_{01}(U_{1t}), \ldots, \theta_{0n}(U_{nt}))^\top$. The QVAR model becomes
This can be further written as an extended QVAR model of order one:
where
Moreover, let $\Gamma_{t,0} \coloneqq I_N$ and, for $k \geq 1$, define
Then, we obtain the moving-average expression
if the following assumption is satisfied.
Assumption (ref) leads to $\left\|\Gamma_{t,k}\right\|_2 \to 0$ as $k \to \infty$, which gives the above expression under Assumption (ref)(i). In addition, under Assumption (ref), we have
Note that $\mathbb{E}[\Gamma_{t,k}]$ is independent of $t$ under Assumption (ref), implying the mean stationarity.
Assumption (ref) restricts the relative magnitudes of autoregressive parameters. However, its implications for individual coefficients are difficult to comprehend. The next lemma provides an alternative, easier-to-interpret sufficient condition.
The condition $\varrho \in [0,1/p)$ implies that the dependence on past outcomes needs to be sufficiently weak to achieve the stationarity.
As shown in Lemma (ref), the SQVAR coefficient functions are monotonically increasing by construction. To preserve this monotonicity, we estimate the coefficient functions globally using a series approximation, rather than performing pointwise quantile regression repeatedly.
Let $\{b_h: h = 1,2, \ldots\}$ be a sequence of monotone-spline basis functions, such as I-splines (ramsay1988monotone) or monotone B-splines (he1998monotone), and denote $\bm{b}_H(\tau) = (b_1(\tau), \ldots, b_H(\tau))^\top$, where $H$ denotes the order of the basis functions chosen by the researcher. The first element of $\bm{b}_H(\tau)$ is typically a constant. Suppose that there exist $H \times 1$ vectors $\gamma_{0i}$ and $\gamma^{(j)}_{li}$, for $j=1,\ldots,p$, such that
where the approximation errors vanish as $H \to \infty$. Define $\underbracket{\bm B_H(\tau)}_{(N + 1) \times H(N + 1)} = \text{diag}(\bm{b}_H(\tau)^\top, \ldots, \bm{b}_H(\tau)^\top)$, and $\underbracket{\bm{\gamma}_i}_{H(N + 1) \times 1} = (\gamma_{0i}^\top, \gamma_{1i}^{(1)\top}, \ldots, \gamma_{ni}^{(1)\top}, \ldots, \gamma_{1i}^{(p)\top}, \ldots, \gamma_{ni}^{(p)\top})^\top$. To be consistent with the sparsity structure considered here, $\bm{\gamma}_i$ satisfies
Then, the SQVAR model is approximated as follows:
where $\xi_t(U_{it}) \coloneqq C_t \otimes \bm{b}_H(U_{it})$.
To estimate $\bm{\gamma}_i$, we perform a penalized monotone quantile regression of $y_{it}$ on $\xi_t(\cdot)$. As for the choice of the penalty function, we adopt the SCAD penalty (fan2001variable) due to its ability to automatically select relevant variables while producing asymptotically unbiased estimates. Based on (ref), we consider the following penalty structure:
where $s_\lambda$ is the SCAD penalty function defined as
Here, $\lambda$ is a penalty parameter satisfying $\lambda \to 0$ as $T \to \infty$, and $a$ is a pre-specified constant, which is typically set as $a = 3.7$.
When one uses the coordinate system (ref) in Example (ref), the restriction (ref) is greatly simplified as $\gamma^{(j)}_{li} = \gamma_{0i}$ for all $(l, j) \in \mathcal S_{0i}$. Accordingly, the penalty term takes the following form in this case:
which is conceptually similar to the group-SCAD penalty in wang2007group, wang2008variable.
The choice of $\lambda$ governs both the sparsity pattern and the estimation quality. If $\lambda$ is too small, irrelevant lag effects may remain in the model, leading to higher variance and potential overfitting. If $\lambda$ is too large, relevant lag effects may be excluded, resulting in larger bias and reduced explanatory power. Thus, $\lambda$ must be carefully chosen in a data-adaptive manner to deal with the bias-variance tradeoff. Below, we propose a BIC-type method for selecting the optimal value of $\lambda$.
Let $\{\tau_1, \ldots, \tau_L\}$ be a set of equally spaced grid points in $[0,1]$:
where $L \to \infty$ as $T \to \infty$. Then, we define our SCAD-penalized SQVAR estimator as
where $\mathcal{G}_H = \prod_{j = 1}^{N + 1}\mathcal{G}_{j,H}$ with $\mathcal{G}_{j,H}$ denoting the $j$-th parameter space. The choice of parameter space varies with the basis functions employed. For any basis functions, we assume $|\bm{b}_H(u)^\top \gamma| < \infty$ uniformly in $u \in [0,1]$ for any $\gamma \in \mathcal{G}_{j,H}$, uniformly in $H$. Additionally, for I-splines, each $\mathcal{G}_{j,H}$ must be a subset of $\{\gamma = (\gamma_1, \ldots, \gamma_H) \in \mathbb{R}^H : \gamma_h \ge 0, \; h \in \{2,3, \ldots, H\}\}$, where the first component of $\gamma$ corresponds to the "location" of the function. Similarly, for monotone B-splines, each $\mathcal{G}_{j,H}$ is a subset of $\{\gamma = (\gamma_1, \ldots, \gamma_H) \in \mathbb{R}^H : \gamma_2 \le \cdots \le \gamma_H\}$. The dependence of $\widehat{\bm{\gamma}}_i(\lambda)$ on $\lambda$ is omitted when there is no confusion, and we simply write $\widehat{\bm{\gamma}}_i$.
Once $\widehat{\bm \gamma}_i$ is obtained, we can estimate $\bm \Phi_i(\tau)$ by $\widehat{\bm \Phi}_i(\tau) \coloneqq \bm B_H(\tau) \widehat{\bm \gamma}_i$. Then, the original QVAR coefficients can be recovered by $\widehat{\bm \Theta}_i(\tau) \coloneqq \bm V^{-1} \widehat{\bm \Phi}_i(\tau)$. Moreover, the set of inactive coefficients $\mathcal{S}_{0i}$ is estimated as
In particular, in the case of coordinate system (ref), the estimator of $\mathcal{S}_{0i}$ is given by
and that of the original QVAR coefficients are
where $\widehat{\mathcal{S}}_{1i} \coloneqq \left\{ (l, j) \in \mathcal S : \left\| \widehat \theta_{li}^{(j)} \right\|_2 \neq 0 \right\}$.
In the analysis of VAR models, determining the lagged variables to include in the system is always a central concern. Including too many lags results in inefficient estimates, while omitting relevant lags leads to biased parameter estimates. In our SCAD penalization procedure, the choice of lags is controlled by $\lambda$. Considering the consistency of BIC in selecting the true model in the context of quantile regression (e.g., lian2012note, lee2014model), we employ the following BIC-type criterion:
where $\widehat s_{1i}(\lambda)$ denotes the cardinality of $\widehat{\mathcal S}_{1i}$ under $\lambda$. Then, the optimal $\lambda$ can be obtained by minimizing ${\rm BIC}(\lambda)$.
We summarize our SQVAR estimation procedure as follows:
When conducting prediction or impulse response analysis based on the estimated QVAR model, it is necessary to recover the joint distribution of $\bm U_t$. If the dimension $n$ is small, a nonparametric approach may be feasible. However, for a large $n$, which could occur in practice, any nonparametric method is likely to suffer from the curse of dimensionality. Thus, we suggest a parametric approach and assume that the joint distribution is characterized by a known copula function $\Pi_\kappa$ with correlation parameter(s) $\kappa \in \mathcal{K}$. The parameter space $\mathcal{K}$ depends on the choice of copula. For example, in the case of a multivariate Gaussian copula with uniform correlation, $\kappa$ corresponds to the correlation coefficient and $\mathcal{K} = [-1,1]$.
To estimate the copula parameter, we first need to recover the rank of each observation, $U_{it}$, by inverting the estimated quantile function:
where $\widehat Q_{y_{it}}(\tau \mid \mathcal{F}_{t-1}) \coloneqq \xi_t(\tau)^\top \widehat{\bm{\gamma}}_i$. Note that the inverse function is well-defined because of the strict monotonicity of $\widehat Q_{y_{it}}(\tau \mid \mathcal{F}_{t-1})$ under the proposed approach. Once $\{(\widehat{U}_{1t}, \ldots, \widehat{U}_{nt})\}_{t=1}^T$ are obtained, the copula parameter $\kappa$ is estimated by the maximum likelihood:
where $f_\kappa$ denotes the copula density. Studying the asymptotic properties of $\widehat \kappa$ is beyond the scope of this paper.
In this section, we first study the convergence rate and the limiting distribution of our SQVAR estimator. For expositional simplicity, the basic model assumptions introduced in Section (ref) (Assumptions (ref)-(ref)) are maintained implicitly throughout the following discussion. In addition, we introduce the following assumptions.
Assumption (ref)(i) ensures the boundedness of the coefficient functions, and hence the conditional quantile function as well, over $[0,1]$, in line with Assumption (ref)(i). In view of Assumption (ref)(ii), we eventually require a specific smoothness condition, not mere continuity (e.g., a H\"{o}lder class of appropriate order). However, for simplicity and generality, we do not impose such a condition here. Assumption (ref)(ii) is a weak technical requirement.
Assumption (ref) should be standard in the literature on quantile regression. Assumption (ref)(i) implies the boundedness of the basis functions, which yields $\sup_{u \in [0,1]} ||\bm{b}_H(u)|| = O(\sqrt{H})$. The constant $\pi$ in condition (ref)(ii) represents the smoothness of the functions $\theta_{0i}$ and $\theta_{li}^{(j)}$. Since there is a rich discussion and guidance on this topic in the literature, we omit the details here (see, for example, CHEN20075549). Condition (ref)(iii) is standard.
Assumption (ref) is a collection of conditions requiring the nonsingularity of some matrices. Note that the dimensions of these matrices increase with $H$. Hence, Assumption (ref) implicitly requires $L$ to increase to infinity along with $H$. Assumption (ref) controls the rates of the tuning parameters $\lambda$ and $H$.
Under the above assumptions, we obtain the following convergence results.
A notable implication of Theorem (ref) is that the coefficient functions can be estimated at the parametric rate in the $L_2$-norm. This result is intuitively understandable when viewing our estimator as essentially connecting pointwise QR estimates (which themselves converge at the parametric rate) across different quantile levels using splines. A similar finding is reported in ando2025functional.
Next, we study the limiting distribution of $(\widehat \theta_{0i}(\tau), \widehat \theta_{li}^{(j)}(\tau))$. However, deriving the limiting distribution of the SQVAR estimator under a general coordinate system complicates the analysis and is less convenient for applications. Therefore, in the next two theorems, we restrict our attention to the coordinate system given in (ref).
The next theorem states that, if the penalty parameter is chosen appropriately as in Assumption (ref)(i), the inactive set of coefficient functions can be consistently estimated.
Furthermore, in view of Theorem (ref)(ii), the result in Theorem (ref) also implies the consistency of $\widehat{\mathcal S}_{1i}$. Given these results, we can derive the asymptotic distribution of the active coefficient functions by restricting our attention to the event $\{\widehat{\mathcal S}_{1i} = \mathcal S_{1i}\}$, since this event occurs with probability approaching one.
In this subsection, we establish the consistency of the BIC criterion (ref). We continue to assume that the coordinate system (ref) is employed. Define the pseudo true coefficient parameter under a potentially misspecified active coefficient set $\mathcal S_1 \subseteq \mathcal S$:
and $\epsilon_{it}(\tau, \mathcal S_1) \coloneqq y_{it} - \xi_t(\tau)^\top \bm{\gamma}_i(\mathcal S_1)$. Furthermore, let
and denote the estimator of the true active coefficient set $\mathcal S_{1i}$ under $\widehat \lambda$ by $\widehat{\mathcal S}_{1i}(\widehat \lambda)$.
Note that overfitting (i.e., $\mathcal S_1 \supset \mathcal S_{1i}$) only results in a loss of estimation efficiency and is therefore less problematic than underfitting (i.e., $\mathcal S_1 \subset \mathcal S_{1i}$). To study the behavior of the estimator under underfitted misspecified models, we introduce the following additional assumptions.
Assumptions (ref)(i) and (ii) are parallel with Assumptions (ref)(i) and (ref)(iii), respectively. They are used to derive the convergence rate of the underfitted QR estimator to its pseudo-true value.
Theorem (ref) shows that the consistency of lag selection is achieved by our SCAD-penalized estimator, with the penalty parameter chosen according to (ref). Note that directly solving the minimization problem in (ref) is computationally intensive. Thus, in practice, we suggest using a grid search.
In this section, we study impulse response analysis based on the QVAR model. Impulse responses measure how the system of time series reacts to exogenous shocks given to a series. In the context of QVAR, we consider the shock as a shift in the rank variable $U$.
We propose two types of impulse response analysis. The first type is the generalized impulse response analysis, which can be viewed as a QVAR counterpart to pesaran1998generalized. In this approach, the response to a shock is measured in terms of expected outcomes. While this type of analysis is standard in conventional VAR models, it averages out the dependence at different quantile levels. On the other hand, there are situations on which a particular quantile is of interest. For instance, when assessing Value-at-Risk (VaR) dynamics, it is essential to focus on the tail behavior. To address this need, we propose a second type of impulse response analysis: a scenario-based quantile impulse response approach.
In the literature, Lee2021 proposed a framework for examining the effects of shocks on time series through a quantile impulse response function. However, their approach is subject to the quantile crossing problem, which undermines the reliability of the results. In contrast, our proposed methods are based on non-crossing quantile functions by construction.
We define the generalized impulse response function as follows:
That is, $\text{IRF}_{ij}(h,\tau^\star)$ represents the change in the $h$-step ahead expected outcome $y_{i,t+h}$ when an exogenous shock that shifts the rank variable $U_{jt}$ to a specific quantile level $\tau^\star$ occurs, holding all other information at time $t-1$ fixed. This definition generalizes the conventional impulse response function by characterizing the magnitude of the shock through a shift in the rank variable, rather than through the size of the innovation term.
The computation of $\text{IRF}_{ij}(h,\tau^\star)$ involves Monte Carlo simulation. Specifically, $\mathbb{E}[y_{i,t + h} \mid U_{jt} = \tau^\star, \mathcal{F}_{t-1}]$ can be estimated through the following procedure:
The estimation of $\mathbb{E}[y_{i,t + h} \mid \mathcal{F}_{t-1}]$ can be done similarly by drawing $\{U^*_{it}\}_{i \in [n]}$ from the unconditional copula $\Pi_{\widehat \kappa}$ in Step 1.
ando2024scenario recently proposed a scenario-based quantile network connectedness framework to study how liquidity shocks propagate across U.S. large bank holding companies, particularly under stress scenarios like the COVID-19 pandemic. Based on their idea, we can conduct a scenario-based impulse response analysis in our context.
Consider a rank trajectory $\bm \tau_{ih} \coloneqq \{\tau_{i0}, \tau_{i1}, \ldots, \tau_{ih}\}$ for time series $i$ of length $h + 1$. We refer to this $\bm \tau_{ih}$ as the “scenario" of interest for the $i$-th series. The scenario-based impulse response analysis then investigates how the forecasts of future outcomes would become if the rank variable evolves according to the specified scenario.
The choice of scenario can be arbitrary. For example, if we are interested in particular historical events that occurred at time $t$ for time series $i$, we first estimate the actual quantile levels $\{U_{it}, U_{i,t+1}, \ldots, U_{i,t+h}\}$ from our estimated model. Then, we can use the estimated sequence $\{\widehat U_{it}, \widehat U_{i,t+1}, \ldots, \widehat U_{i,t+h}\}$ as $\bm \tau_{ih}$. As an alternative to focusing on real historical scenarios, one may consider a hypothetical scenario $\bm \tau_{ih}$ based on empirical interest. In this case, it is possible to directly incorporate the assumptions or beliefs of researchers into the scenario.
The scenario-based impulse response analysis can be implemented in the following procedure:
When there is a reference scenario $\overline{\bm \tau}_{ih}$ to measure the impact of a scenario of interest $\bm \tau_{ih}$, we define the impulse response function in this context as
That is, $\text{IRF}_i(h,\bm \tau_{ih}, \overline{\bm \tau}_{ih})$ captures the deviation of the forecasted outcome under the scenario of interest from the reference scenario. Note that, in this scenario-based analysis, it is not necessary to recover the joint distribution of $\bm U_t$.
In this subsection, we numerically compare the performance of the proposed penalized SQVAR method and the standard QR method. Since asymptotically the standard QR also produces quantile curves without crossing, both estimators should perform similarly when the sample size is sufficiently large. However, in small or moderate samples, the monotonicity constraint explicitly imposed on our SQVAR method may affect its finite-sample performance. In particular, it is generally known that imposing a shape restriction helps reduce the estimation variance under weak identification. The following simulation results suggest that this is indeed the case for our estimator.
For the DGP, we consider a trivariate QVAR model with lag order $p=2$, where the first-order lags are all active and the second-order lags are all inactive. The active and inactive lags are selected using the BIC criterion with grid search. The strength of model identification is controlled by scaling the QVAR coefficients $\theta_{li}^{(1)}$ by $1/b$, where $b$ is chosen from $b \in \{1,2, \ldots, 6\}$. That is, the larger $b$ is, the weaker the identification. In Figure (ref), panel (a) reports the boxplots of the root mean squared error (RMSE) averaged over the coefficients for the two estimators when $T = 200$, and panel (b) reports those for $T = 600$, where the number of Monte Carlo replications is set to 500. For more detailed information about the DGP and estimation procedure, please see Appendix (ref).
From Figure (ref), we find that the standard QR estimator tends to perform slightly better than our penalized SQVAR estimator when $b$ is small. It is not surprising in general that an unconstrained estimator achieves smaller RMSE than a constrained estimator. However, as $b$ gets larger, this relationship reverses. In particular, when the sample size is smaller, the robustness of our shape-restricted method to weak identification becomes more apparent. Considering the fact that the standard QR estimator suffers from the quantile crossing problem especially when the sample size is small, as reported in Table (ref), we may conclude that the SQVAR method should be a promising alternative to the standard QR method.
In this subsection, we assess the accuracy of the penalized SQVAR estimator and the BIC lag selection. As in the previous subsection, the experimental setup and detailed results are relegated to Appendix (ref); here we report only the main findings.
Overall, we find that the estimator performs satisfactorily well in terms of RMSE. Estimation errors decrease as the sample size $T$ increases. The number of inner knots for the spline basis appears to have minimum impact on the accuracy. Increasing the number of quantile grid points improves the estimation accuracy slightly.
For the performance of the BIC lag selection, we find that the frequency of correctly identifying the set of active lags increases with the sample size $T$, which is consistent with Theorem (ref). In particular, the frequency that the selected active set contains the true $\mathcal S_1$ reaches 100% when $T$ is large. For the number of inner knots, a more parsimonious estimator tends to perform better possibly due to its smaller variance. For the choice of penalty parameter $\lambda$, although it is highly dependent on the DGP, the results suggest that exploring values in the neighborhood of $\lambda = c_\lambda \ln T / \sqrt{T}$ with $c_\lambda \in [0.5, 1]$ may be a reasonable default.
We apply the proposed model and method to the daily returns of six U.S. portfolios constructed as intersections of two portfolios formed by size (market equity) and three portfolios formed by momentum based on prior (2-12) returns. The data are publicly available from Professor French's website: \url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}.
Specifically, the six portfolios used are: Small-Up, Small-Md, Small-Lo, Big-Up, Big-Md, and Big-Lo. Here, Small/Big refers to the size of market equity, and Up/Md/Lo indicates the level of prior returns, which is referred to as the momentum factor in the literature. For precise definitions of these classifications, see the website cited above. Understanding the dynamic relationships among these portfolios is crucial for investors, policymakers, and risk managers aiming to diversify portfolios, assess systemic risk, and implement effective macroprudential policies.
The data set consists of daily returns for the six portfolios over the time period from January 11, 1995 to June 30, 2025 for $T = 7668$ observations. Figure (ref) plots the time series of these six portfolios. Similar to other return data, the series show no trend, but clear volatility clustering. The volatility patterns of the portfolios seem to depend on the momentum factor.
We apply the proposed SQVAR estimator with the coordinate system given in (ref). The estimation procedure follows the same steps as those in the previous subsections. For the choice of tuning parameters, we set $\sharp k = 1$, $L = 30$, and $p = 6$. The penalty parameter $\lambda$ is selected by minimizing $\mathrm{BIC}(\lambda)$ via a grid search. The lag selection result indicated that all coefficient functions, except for $\theta_{\text{Big-Md}, \text{Small-Up}}^{(6)}$, $\theta_{\text{Big-Md}, \text{Big-Up}}^{(6)}$, and $\theta_{\text{Big-Up}, \text{Big-Up}}^{(5)}$, contribute to our model. For the estimation results of the coefficient functions, for illustrative purposes, we show two pairs of the estimated coefficients in Figure (ref). In the figure, the shaded areas represent pointwise 95% confidence intervals.
In Figure (ref), panel (a) shows the estimated coefficient functions across quantiles $\tau \in (0,1)$ for $\theta_{\text{Small-Md}, \text{Small-Md}}^{(1)}$ and $\theta_{\text{Small-Up}, \text{Small-Up}}^{(1)}$, which correspond to the first-order own-lag effects of these portfolios. Both estimated functions are almost constant over $\tau$, suggesting that the dependence structure they represent is fairly stable and does not vary strongly with market conditions. The coefficient $\theta_{\text{Small-Md}, \text{Small-Md}}^{(1)}$ is consistently negative, around $-0.15$ across most quantiles, with a relatively narrow confidence interval in the central range $(0.25,0.75)$. In contrast, $\theta_{\text{Small-Up}, \text{Small-Up}}^{(1)}$ is consistently positive, hovering around $0.2$, again with tight confidence intervals near the center. A plausible interpretation is that Small-Md firms, characterized by low market capitalization and intermediate past returns, do not generate strong momentum signals and tend to be more exposed to short-term price reversals, leading to negative first-order autocorrelation, whereas Small-Up firms are more likely to attract momentum trading and herding, which can induce positive autocorrelation. The contrast between these two coefficients highlights that different portfolios can exhibit opposite forms of dynamic dependence, and the near constancy across $\tau$ suggests that these effects are structural features of the system under both normal and stressed market conditions.
In panel (b), we show the estimated coefficient functions $\theta_{\text{Small-Lo}, \text{Big-Md}}^{(1)}$ and $\theta_{\text{Small-Up}, \text{Big-Md}}^{(1)}$. In contrast to those in panel (a), both functions depart from constancy in $\tau$, indicating that the strength of the cross-dynamic dependence between the portfolios varies across quantiles. Both estimates are moderately negative over the central quantiles $(0.25,0.75)$, suggesting a robust negative cross dependence under normal market conditions. However, the intervals widen substantially at the tails, reflecting greater uncertainty in these regions.
VaR is a widely used risk measure that captures the potential loss in value of a portfolio under normal market conditions. However, when multiple portfolios exhibit interdependence and spillovers, it is crucial to understand the dependence structure at the tails of the return distribution, as shocks in one sector can propagate to others and alter tail risks dynamically. The proposed scenario-based impulse response analysis can be used to quantify the dependence risk.
To demonstrate our scenario-based analysis, we construct two scenarios corresponding to the following historical events:
To quantify the impact of these events, we define a baseline scenario that reflects stable market conditions. Considering the market to be stable when the CBOE Volatility Index (VIX) is low (\url{https://fred.stlouisfed.org/series/VIXCLS}), we set December 13, 2023 as the starting date of the baseline period. Below, we report the scenario-based average responses $6^{-1}\sum_{i = 1}^6 y_{i, t+h}^*(\bm \tau_{ih})$ and $6^{-1}\sum_{i = 1}^6 y_{i, t+h}^*(\overline{\bm \tau}_{ih})$, for $h = -5, \ldots, -1, 0 \text{ (event date) } , 1, \ldots, 22$, such that the difference between them corresponds to the average of $\text{IRF}_i(h,\bm \tau_{ih}, \overline{\bm \tau}_{ih})$. For negative $h$, we simply report the average of the estimated conditional quantiles on that date. The choice of $h$ = 22 corresponds roughly to the number of trading days in a month.
Figure (ref) presents scenario-based impulse responses under the 2008 financial crisis and the onset of the COVID-19 pandemic, each evaluated relative to the baseline scenario. In the top panel, under the financial crisis scenario, the portfolio returns exhibit large and persistent fluctuations over the forecast horizon. The responses alternate between negative and positive values, with conspicuous downward movements around $h=3$, $10$, and $22$. Relative to the comparatively flat baseline path, the crisis scenario shows substantially higher variability, indicating elevated risk dynamics during this period.
In the bottom panel, under the COVID-19 scenario, the responses display even larger and more irregular movements. Negative responses are especially strong immediately after the event date, with several troughs below $-4$, followed by sharp reversals into positive return. As in the financial crisis scenario, volatility remains elevated throughout the horizon, and there is no clear convergence toward the baseline path.
Overall, Figure (ref) illustrates that both historical events generate substantial departures from baseline tail risk dynamics. These results underscore the usefulness of the scenario-based impulse response framework for quantifying how systemic shocks can alter the dependence structure of portfolio returns.
In some applications, it is natural to consider a high-dimensional QVAR model, where the number of cross-sectional units $n$ is large. This section provides an additional result to address this challenging setting. Note that the result presented below applies to a broad class of high-dimensional autoregressive models, including our QVAR as a special case.
We use a single index $m=(\ell,j)$ with $\ell\in [n]$ (series) and $j\in [p]$ (lag order), so that $X_{m,t} \coloneqq y_{\ell,t-j}$. Let $\tau\in(0,1)$ be a fixed quantile level. For each time series $i \in [n]$ and predictor index $m=(\ell,j)$, we write the marginal conditional $\tau$-th quantile as $Q_{y_{it}}(\tau \mid X_{m,t}\big)$, and denote the unconditional $\tau$-th quantile of $y_{it}$ by $Q_{y_{it}}(\tau) \coloneqq \inf\{y:\Pr(y_{it}\le y)\ge\tau\}$. Based on the following observation \[ Y \text{ and } X_{m,t} \text{ are independent } \iff Q_{y_{it}}(\tau \mid X_{m,t}\big) - Q_{y_{it}}(\tau) = 0 \quad \text{for all } \tau \in (0, 1), \] HeWangHong2013 developed a variable screening procedure by introducing a quantile-adaptive framework for high-dimensional heterogeneous data. Similar to HeWangHong2013, we define the active set for the $i$-th series by \[ M_{i}=\{ m:\ Q_{y_{it}}(\tau \mid \mathcal F_{t-1})\ \text{functionally depends on}\ X_{m,t}\, ~ {\rm for~some}~\tau \in (0,1)\}. \] To quantify the effect of $X_{m,t}$ on $y_{it}$ at a specific $\tau$ and similar to HeWangHong2013, we consider the following measure: \[ d_{i,m,\tau}(X_{m,t}) \coloneqq Q_{y_{it}}(\tau \mid X_{m,t}\big) - Q_{y_{it}}(\tau\big). \] To estimate $Q_{y_{it}}(\tau \mid X_{m,t}\big)$ and $Q_{y_{it}}(\tau)$, we can use $\widehat Q_{y_{it}}(\tau \mid X_{m,t}\big) = \widehat b_{i,m,0}(\tau) + \widehat b_{i,m,1}(\tau)X_{m,t}$ with $\{\widehat b_{i,m,0}(\tau), \widehat b_{i,m,1}(\tau)\}= \operatorname*{argmin}_{b_{0},b_{1}} \sum_{t=1}^T \rho_\tau(y_{it} - b_{0}-X_{m,t}b_{1})$, and $\widehat Q_{y_{it}}(\tau)=\widehat F^{-1}_{y_i}(\tau)$, which is the $\tau$-th sample quantile function based on $\{y_{i1}, \dots, y_{iT}\}$.
Define \[ \widehat{d}_{i,m,\tau}(X_{m,t}) \coloneqq \widehat Q_{y_{it}}(\tau \mid X_{m,t}\big)-\widehat Q_{y_{it}}(\tau) = \widehat b_{i,m,0}(\tau) + \widehat b_{i,m,1}(\tau)X_{m,t}-\widehat F^{-1}_{y_i}(\tau). \] We expect $\widehat{d}_{i,m,\tau}(X_{m,t})$ to be close to zero if $X_{m,t}$ is independent of $y_{it}$. The independence screening is based on the magnitude of the estimated marginal components \[ \|\widehat{d}_{i,m,\tau}\|_2^2 \coloneqq \frac{1}{T} \sum_{t=1}^T \widehat{d}_{i,m,\tau}(X_{m,t})^2. \] More specifically, we select the subset of variables $\|\widehat{d}_{i,m,\tau}\|_2^2 \ge \nu_T $ for some $\tau \in (0,1)$ where $\nu_T$ is a predefined threshold. In practical implementation, we employ discrete points of $\tau \in (0,1)$ as $\mathcal T :=\{ \tau_1,\ldots,\tau_A\}$, and check the condition $\|\widehat{d}_{i,m,\tau}\|_2^2 \ge \nu_T$ at each grid point: \[ \widehat M_{i} = \left\{ m : \|\widehat{d}_{i,m,\tau_a}\|_2^2 \ge \nu_T \quad {\rm for~some}~ \tau_a \in \mathcal T \right\}. \] Note that unlike HeWangHong2013, we are trying to detect the active set by considering the entire $\tau$ under dependent time series data. We impose the following assumptions.
Theorem (ref) indicates that, with probability tending to one, the screened model contains all true active predictors for all quantile levels and all equations. Therefore, when the dimension $n$ of QVAR is very high, we can apply the screening prior to performing the simplex transformation.
Similar to HeWangHong2013, we can study the number of selected variables after the screening. Define $Z_{m,it} \coloneqq\widehat d_{i,m,\tau}(X_{m,t})^2 - \mathbb E[\widehat d_{i,m,\tau}(X_{m,t})^2]$. For each candidate index $m$ define the screening score \[ S_{i,m} = \mu_{i,m} + \frac{1}{T}\sum_{t=1}^T Z_{m,it}, \qquad \mu_{i,m} = \mathbb E[S_{i,m}]. \]
The above theorem suggests that if a signal is well separated, then the model obtained after screening is $O(1)$. We can relax this result by allowing the size of active signals $|M_i|$ gradually increases with $T$. However, this is out of the scope of this paper.
In this paper, we introduced a simplex quantile VAR (SQVAR) framework that enforces non-crossing quantile curves in QVAR models. For estimation, we developed a SCAD-penalized monotone series method. We established asymptotic properties of the estimator, including the rate of convergence, asymptotic normality, and the consistency of model selection under some regularity conditions. In addition, we derived a BIC-type criterion for choosing the SCAD penalty parameter and proved its theoretical validity. Furthermore, we extended impulse response analysis to the proposed framework by developing two novel approaches. An empirical application to U.S. portfolio data demonstrated that the proposed method captures heterogeneous dynamics across return distributions and provides useful insights into market risk, especially during financial crises.
There are several directions for future research. An important open question is to consider the joint estimation of multiple time series, which could improve efficiency by explicitly accounting for cross-sectional dependence (e.g., jun2009efficient, petrella2019joint). In addition, several inferential problems remain open, such as statistical inference for the impulse responses, uniform inference for the estimated coefficient functions and quantile functions, and error variance decomposition in the QVAR context. Addressing these topics would further enhance the scope and applicability of the proposed framework.