Reduced-Rank Envelope Vector Autoregressive Model

Abstract

The standard vector autoregressive (VAR) models suffer from overparameterization which is a serious issue for high-dimensional time series data as it restricts the number of variables and lags that can be incorporated into the model. Several statistical methods, such as the reduced-rank model for multivariate (multiple) time series (Velu, Reinsel, and Wichern; Reinsel and Velu; Reinsel, Velu, and Chen) and the Envelope VAR model (Wang and Ding), provide solutions for achieving dimension reduction of the parameter space of the VAR model. However, these methods can be inefficient in extracting relevant information from complex data, as they fail to distinguish between relevant and irrelevant information, or they are inefficient in addressing the rank deficiency problem. We put together the idea of envelope models into the reduced-rank VAR model to simultaneously tackle these challenges, and propose a new parsimonious version of the classical VAR model called the reduced-rank envelope VAR (REVAR) model. Our proposed REVAR model incorporates the strengths of both reduced-rank VAR and envelope VAR models and leads to significant gains in efficiency and accuracy. The asymptotic properties of the proposed estimators are established under different error assumptions. Simulation studies and real data analysis are conducted to evaluate and illustrate the proposed method.

Keywords:

• Envelope model
• Reduced-rank autoregression
• Vector autoregressive model

1 Introduction

With the recent rapid development of information technology, high-dimensional time series data are routinely collected in various fields, like finance, economics, digital signal processing, neuroscience, and meteorology. Classical vector autoregressive (VAR) models are widely used for modeling multivariate time series data due to their ability to capture dynamic relationships among variables in a system and their flexibility. These models are discussed in many time series textbooks, including Hamilton (1994), Lütkepohl (2005), Tsay (2014), Box et al. (2015), Wei (2019), and others. However, the VAR model suffers from overparameterization, particularly when the number of lags and time series increases. While there are several statistical methods for achieving dimension reduction in time series (Park and Samadi 2014, 2020; Cubadda and Hecq 2022b; Samadi and DeAlwis 2023a), they can be inefficient in extracting relevant information from complex data. This is because they fail to differentiate between important and unimportant information, which can obscure the material and useful information. Envelope methods (Cook, Li, and Chiaromonte 2010) use reduced subspaces to link the mean function and dispersion matrix through novel parameterizations. By identifying and removing irrelevant information, the envelope model is based only on useful information and is therefore more efficient. The literature presents various methods for multivariate time series analysis, including the canonical transformation (Box and Tiao 1977), reduced-rank VAR models (Velu, Reinsel, and Wichern 1986; Reinsel, Velu, and Chen 2022), scalar component models (Tiao and Tsay 1989), LASSO regularization of VAR models (Shojaie and Michailidis 2010; Song and Bickel 2011), sparse VAR model based on partial spectral coherence (Davis, Zang, and Zheng 2016), factor modeling (Stock and Watson 2005; Forni et al. 2005; Lam and Yao 2012), envelope VAR models (Wang and Ding 2018; Herath and Samadi 2023a, 2023b), nonlinear VAR models (Samadi, Hajebi, and Farnoosh 2019), and tensor-structure modeling for VAR models (Wang et al. 2022a; Wang, Zheng, and Li 2021), among others.

Let ${\mathbf{y}}_t={(y_{1t},\dots ,y_{qt})}^{\prime }$ be a q-dimensional VAR process of order p, VAR(p), given as (1) $${\mathbf{y}}_t=\mathit{\alpha }+{\mathit{\beta }}_1{\mathbf{y}}_{t-1}+{\mathit{\beta }}_2{\mathbf{y}}_{t-2}+\cdots +{\mathit{\beta }}_p{\mathbf{y}}_{t-p}+{\mathit{\epsilon }}_t, t=1,\dots ,T,$$(1) where ${\mathit{\beta }}_i\in {\mathbb{R}}^{q\times q}$ are coefficient matrices, and ${\mathit{\epsilon }}_t\in {\mathbb{R}}^q$ is a vector white noise process with mean $\mathbf{0}$, and covariance matrix $\mathbf{\Sigma }$, that is, ${\mathit{\epsilon }}_t\sim WN(\mathbf{0},\mathbf{\Sigma })$, and T denotes the sample size. Suppose L is the lag operator and let $\mathit{\beta }(L)$ be the characteristic polynomial function of the model given as $\mathit{\beta }(L)={\mathbf{I}}_q-{\mathit{\beta }}_{1}L-\cdots -{\mathit{\beta }}_p{L}^p$. Then, the time series in (1) is stationary if all the roots of $\det (\mathit{\beta }(L))=0$ are greater than one in modulus.

Modeling high-dimensional multivariate time series is always challenging due to the dependent and high-dimensional nature of the data. Even for moderate dimensions q and p, performing the estimation can be difficult (De Mol, Giannone, and Reichlin 2008; Carriero, Kapetanios, and Marcellino 2011; Koop 2013). The number of coefficient parameters in model (1), $q^{2}p$, can dramatically increase with the dimension q and the lag order p. Therefore, to improve estimation and make inferences on high-dimensional VAR models, it is necessary to restrict the parameter space to a reasonable and manageable number of parameters. To this end, note that the VAR(p) model in (1) can be rewritten as a VAR(1) model as (2) $${\mathbf{y}}_t=\mathit{\alpha }+\mathit{\beta }{\mathbf{x}}_t+{\mathit{\epsilon }}_t,$$(2) where ${\mathbf{x}}_t\in {\mathbb{R}}^{qp\times 1}$ is the vector of lagged variables, that is, ${\mathbf{x}}_t={({\mathbf{y}}_{t-1}^{\prime },{\mathbf{y}}_{t-2}^{\prime },\dots ,{\mathbf{y}}_{t-p}^{\prime })}^{\prime }$, and $\mathit{\beta }=({\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_p)\in {\mathbb{R}}^{q\times qp}$ which encompasses the coefficient matrices of lagged variables.

To address the overparameterization issue of the VAR model, we assume that the autoregressive coefficient matrix $\mathit{\beta }$ in model (2) has a reduced-rank structure, similar to the standard reduced-rank regression (RRR) model (Velu, Reinsel, and Wichern 1986; Anderson 1999, 2002; Reinsel, Velu, and Chen 2022). This technique improves the accuracy of the estimation of $\mathit{\beta }$ by reducing the dimensionality of ${\mathbf{y}}_t$ and ${\mathbf{x}}_t$. The envelope model (Cook, Li, and Chiaromonte 2010) is a new dimension reduction technique that has a different perspective to achieve efficient estimation by linking the mean function and covariance matrix and using the minimal reducing subspace of the covariance matrix. By combining these two methods, we propose a novel dimensionality reduction model for the VAR process called the reduced-rank envelope vector autoregressive (REVAR) model, which extends the idea of envelopes to the reduced-rank VAR model (RRVAR). As a result, our proposed method is more efficient and parsimonious than both methods alone and outperforms both.

The reduced-rank problem first arose in the multivariate regression analysis to achieve dimension reduction by restricting the rank of the coefficient matrix. Anderson (1951) considered the RRR problem for fixed predictors. Then, Izenman (1975) and Reinsel and Velu (1998) studied the RRR in detail. Asymptotic properties for RRR are discussed by Stoica and Viberg (1996), and Anderson (1999). The RRR models have also been investigated by several other authors (Reinsel and Velu 1998; Negahban and Wainwright 2011; Chen, Dong, and Chan 2013; Basu, Li, and Michailidis 2019; Raskutti, Yuan, and Chen 2019).

Reduced-rank model for multivariate (multiple) time series (Velu, Reinsel, and Wichern 1986; Reinsel and Velu 1998; Reinsel, Velu, and Chen 2022) often arises in the multivariate statistics literature when coefficient matrices have low-rank structures. The reduced-rank VAR (RRVAR) model is considered by imposing a low-rank structure on the coefficient matrix of model (2), that is, rank $(\mathit{\beta })=d<q$. As a result, the number of parameters is decreased and the estimation efficiency is improved. The analysis of RRVAR models has connections with some known methodologies such as principal component analysis (PCA) (Rao 1964; Billard, Douzal-Chouakria, and Samadi 2023) and canonical analysis (CCA) (Box and Tiao 1977; Samadi et al. 2017) to achieve dimension reduction and improve predictions. The asymptotic properties of the reduced-rank (RR) and the ordinary least squares (OLS) estimators of the VAR model are studied by Anderson (2002). The RRVAR model has been extended and combined with other approaches in economic and financial modeling, including, common features and RRVAR models (Franchi and Paruolo 2011; Centoni and Cubadda 2015; Cubadda, Hecq, and Telg 2019), structural analysis through RRVAR models (Bernardini and Cubadda 2015; Carriero, Kapetanios, and Marcellino 2016), Bayesian RRVAR models (Carriero, Kapetanios, and Marcellino 2011; Cubadda, Hecq, and Telg 2019), partial least squares approach (Cubadda and Hecq 2011), the vector error correction model (VECM) under cointegration (Lütkepohl 2005, Part II; Hecq, Palm, and Urbain 2006), multivariate autoregressive index models (MAI) (Cubadda Guardabascio and Hecq 2017; Cubadda, Hecq, and Telg 2019), eteroscedastic VAR models (Hetland, Pedersen, and Rahbek in press), time-varying parameter RRVAR models (Brune, Scherrer, and Bura 2022). A detailed review of RRVAR models is given in Cubadda and Hecq (2022a). These models are distinctly different from envelope models because there is no link between the mean function and covariance matrix (see (6) in Section 2.2).

While the reduced-rank VAR model achieves an effective dimensionality reduction, in many applications there are additional structures that can be exploited to achieve even higher dimensionality reduction with lower statistical error. The response envelope proposed by Cook, Li, and Chiaromonte (2010) is another parsimonious approach to achieve dimension reduction and improve the estimation efficiency and prediction accuracy of standard multivariate regression models by parsimoniously decreasing the number of parameters. This method is useful in eliminating immaterial information present in the responses and predictors. The envelope method is effective even when the coefficient matrix is full rank, whereas the reduced-rank approach offers no reduction in this case. This is because the envelope uses the smallest reducing subspace of the covariance matrix that contains the mean function. There are several extensions of the basic envelope methodology to other contexts (Su and Cook 2011; Cook, Helland, and Su 2013; Cook, Forzani, and Zhang 2015; Cook and Zhang 2015a, 2015b; Su et al. 2016; Li and Zhang 2017; Ding and Cook 2018; Lee, M. and Su 2020; Forzani and Su 2021).

Wang and Ding (2018) extended the envelope regression model proposed by Cook, Li, and Chiaromonte (2010) to the VAR model called the envelope VAR (EVAR) model. The EVAR model provides better performance and is more efficient by removing immaterial information from estimation. Rekabdarkolaee et al. (2020) proposed a spatial envelope model for spatially correlated data in multivariate spatial regression. Samadi and DeAlwis (2023b) introduced the envelope matrix autoregressive (MAR) model (for the MAR model, refer to Samadi 2014). Cook, Forzani, and Zhang (2015) proposed a new parsimonious multivariate regression model by combining Anderson’s (1999) RRR model with Cook, Li, and Chiaromonte’s (2010) envelope model, called the reduced-rank envelope model. We incorporate the idea of envelopes (Cook, Li, and Chiaromonte 2010; Wang and Ding 2018) into Velu, Reinsel, and Wichern’s (1986) RRVAR model and propose a novel efficient parsimonious VAR model for high-dimensional time series data. The proposed reduced-rank envelope VAR model combines the advantages and strengths of both the RRVAR and EVAR models which leads to more accurate estimation and higher efficiency.

We use the following notations and definitions throughout this article. All real k × s matrices are denoted as ${\mathbb{R}}^{k\times s}$. The collection of u-dimensional subspaces in a q-dimensional vector space is called the Grassmannian manifold, indicated by ${\mathcal{G}}_{q,u}$. If $\mathbf{A}\in {\mathbb{R}}^{k\times s}$, then $\text{span}(\mathbf{A})\subseteq {\mathbb{R}}^k$ is defined to be the subspace spanned by the columns of $\mathbf{A}$. If $\sqrt{T}(\widehat{\mathit{\theta }}-\mathit{\theta })\underset{}{\overset{\mathcal{D}}{\to }}\mathcal{N}(\mathbf{0},\mathbf{V})$, then the asymptotic covariance matrix of $\sqrt{T}\widehat{\mathit{\theta }}$ is denoted as $\text{avar}\left(\sqrt{T}\widehat{\mathit{\theta }}\right)=\mathbf{V}$. The projection matrix onto $\text{span}(\mathbf{X})$ in the $\mathbf{V}$ inner product is denoted by ${\mathbf{P}}_{\mathbf{X}(\mathbf{V})}=\mathbf{X}{(\mathbf{X}\prime \mathbf{VX})}^{-1}\mathbf{X}\prime \mathbf{V}$, and the projection matrix onto $\text{span}(\mathbf{X})$ in the identity inner product is denoted as ${\mathbf{P}}_{\mathbf{X}}=\mathbf{X}{(\mathbf{X}\prime \mathbf{X})}^{-1}\mathbf{X}\prime$. Let ${\mathbf{Q}}_{\mathbf{X}(\mathbf{V})}=\mathbf{I}-{\mathbf{P}}_{\mathbf{X}(\mathbf{V})}$, and ${\mathbf{Q}}_{\mathbf{X}}=\mathbf{I}-{\mathbf{P}}_{\mathbf{X}}$. The operator $\text{vec}:{\mathbb{R}}^{k\times s}\to {\mathbb{R}}^{ks}$ stacks the columns of a matrix into a column vector, and the operator $\text{vech}:{\mathbb{R}}^{k\times k}\to {\mathbb{R}}^{k(k+1)/2}$ stacks the lower triangular part, including the diagonal, into a vector. Moreover, for any q × q symmetric matrix U, the expansion matrix ${\mathbf{E}}_q\mathrm{}\in \mathrm{}{\mathbb{R}}^{q^2\times q(q+1)/2}$ is defined such that $\text{vec}(\mathbf{U})={\mathbf{E}}_q\text{vech}(\mathbf{U})$, and the contraction matrix ${\mathbf{C}}_q\in {\mathbb{R}}^{q(q+1)/2\times q^2}$ is defined such that $\text{vech}(\mathbf{U})={\mathbf{C}}_q\text{vec}(\mathbf{U})$ (Henderson and Searle 1979). The Moore-Penrose generalized inverse of ${\mathbf{E}}_q$ is defined as ${\mathbf{E}}_q^{†}={({\mathbf{E}}_q^{\prime }{\mathbf{E}}_q)}^{-1}{\mathbf{E}}_q^{\prime }$. The Kronecker product of two matrices $\mathbf{X}$ and $\mathbf{Y}$ is denoted by $\mathbf{X}\otimes \mathbf{Y}$, and the symbol “∼” means identically distributed.

The rest of the article is structured as follows. Section 2 provides a comprehensive review of the reduced-rank VAR and the envelope VAR models, along with the introduction of our proposed approach, the reduced-rank envelope VAR (REVAR) model. The derivation of the maximum likelihood estimators (MLEs) for the parameters of the REVAR model is presented in Section 3. In Section 4, we establish the asymptotic properties of the proposed REVAR estimators with and without normality assumption and compare them with those obtained from the standard VAR, RRVAR, and EVAR models. The algorithms for selecting the lag order, rank, and envelope dimensions are outlined in Section 5. To assess the performance of our proposed REVAR model under different error assumptions, Section 6 presents the results of extensive simulation studies, comparing it with the reduced-rank VAR, envelope VAR, and standard VAR models. Real economic datasets are analyzed in Section 7, and our conclusions are summarized in Section 8. The supplementary materials contain the proofs of lemmas and propositions, as well as additional simulations.

2 Reduced-Rank Envelope Vector Autoregression

2.1 Reduced-Rank VAR Model

Suppose the coefficient matrix in model (2) is rank deficient, that is, $\text{rank}(\mathit{\beta })=d<q$. As a result, it can be written as a product of two lower dimensional matrices, that is, $\mathit{\beta }=\mathbf{AB}$. Then, the reduced-rank VAR (RRVAR) model proposed by Velu, Reinsel, and Wichern (1986) is given by (3) $${\mathbf{y}}_t=\mathit{\alpha }+\mathbf{AB}{\mathbf{x}}_t+{\mathit{\epsilon }}_t,$$(3) where $\mathbf{A}\in {\mathbb{R}}^{q\times d},\mathbf{B}=[{\mathbf{B}}_1,{\mathbf{B}}_2,\dots ,{\mathbf{B}}_p]\in {\mathbb{R}}^{d\times qp}$, and $\text{rank}(\mathbf{A})=\text{rank}(\mathbf{B})=d$.

The conditional log-likelihood function of model (2) under the assumption of normality of ${\mathit{\epsilon }}_t$ can be written as (the initial values $({\mathbf{y}}_{-p+1},\dots ,{\mathbf{y}}_0)$ are assumed to be given) (4) $$L_T(\mathit{\alpha },\mathit{\beta },\mathbf{\Sigma })\propto \frac{-T}{2}\{\hspace{0.17em}\log \hspace{0.17em}|\mathbf{\Sigma }|+\frac{1}{T}\sum _{t=1}^T{({\mathbf{y}}_t-\mathit{\alpha }-\mathit{\beta }{\mathbf{x}}_t)}^{\prime }{\mathbf{\Sigma }}^{-1}({\mathbf{y}}_t-\mathit{\alpha }-\mathit{\beta }{\mathbf{x}}_t)\}.$$(4)

Then, the maximization of $L_T(\mathit{\alpha },\mathit{\beta },\mathbf{\Sigma })$ is performed under the condition that rank $(\mathit{\beta })=d$, or correspondingly under the RR autoregression parameterization $\mathit{\beta }=\mathbf{AB}$ in model (3). The symbol $\propto$ is used to exclude unnecessary additive constants from the likelihood function. Before we present the MLEs of the RRVAR model (3), we introduce some notations.

Let $\mathbf{\Sigma }=\text{cov}({\mathit{\epsilon }}_t)$, and let ${\mathbf{\Gamma }}_k=\text{cov}({\mathbf{y}}_{t-k},{\mathbf{y}}_t)$, $k\ge 0$, denotes the autocovariance matrix function of ${\mathbf{y}}_t$. Also, let ${\mathbf{\Gamma }}_{*}=\text{cov}({\mathbf{x}}_t,{\mathbf{y}}_t)\in {\mathbb{R}}^{qp\times q}$ and ${\mathbf{\Gamma }}_{(p)}=\text{cov}({\mathbf{x}}_t)\in {\mathbb{R}}^{qp\times qp}$ given as $${\mathbf{\Gamma }}_{*}=\left(\begin{array}{c}{\mathbf{\Gamma }}_1\\ {\mathbf{\Gamma }}_2\\ ⋮\\ {\mathbf{\Gamma }}_p\end{array}\right), {\mathbf{\Gamma }}_{(p)}=\left(\begin{array}{cccc}{\mathbf{\Gamma }}_0& {\mathbf{\Gamma }}_1^{\prime }& \dots & {\mathbf{\Gamma }}_{p-1}^{\prime }\\ {\mathbf{\Gamma }}_1& {\mathbf{\Gamma }}_0& \dots & {\mathbf{\Gamma }}_{p-2}^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{\Gamma }}_{p-1}& {\mathbf{\Gamma }}_{p-2}& \dots & {\mathbf{\Gamma }}_0\end{array}\right).$$

Given a series of observations ${\mathbf{y}}_{-p+1},\dots ,{\mathbf{y}}_0,{\mathbf{y}}_1,\dots ,{\mathbf{y}}_T$, their sample covariance matrices are denoted by ${\widehat{\mathbf{\Gamma }}}_{*}$ and ${\widehat{\mathbf{\Gamma }}}_{(.)}$ with the divisor T. Without loss of generality, suppose the sample lagged vector ${\mathbf{x}}_t$ is centered, therefore, $\overline{\mathbf{x}}=\mathbf{0}$. For ease of notations, we shall drop the subsubscript “t” in ${\mathbf{y}}_t$ and ${\mathbf{x}}_t$ when they are subscripted. Let ${\widehat{\mathbf{\Gamma }}}_k=\frac{1}{T}{\sum }_{t=1}^T({\mathbf{y}}_{t-k}-\overline{\mathbf{y}})({\mathbf{y}}_t-\overline{\mathbf{y}})\prime$ ($k=0,1,\dots$), and ${\widehat{\mathbf{\Gamma }}}_{(p)}=\frac{1}{T}{\sum }_{t=1}^T{\mathbf{x}}_t{\mathbf{x}}_t^{\prime }$ be the sample autocovariance matrices of ${\mathbf{y}}_t$ and ${\mathbf{x}}_t$ at lag k, respectively; and ${\widehat{\mathbf{\Gamma }}}_{*}=\frac{1}{T}{\sum }_{t=1}^T{\mathbf{x}}_t({\mathbf{y}}_t-\overline{\mathbf{y}})\prime$ be the sample cross-covariance matrix between ${\mathbf{x}}_t$ and ${\mathbf{y}}_t$. Moreover, let ${\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}={\widehat{\mathbf{\Gamma }}}_0-{\widehat{\mathbf{\Gamma }}}_{*}^{\prime }{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1}{\widehat{\mathbf{\Gamma }}}_{*}$ be the sample covariance matrix of residuals, and ${\widehat{\mathbf{\Gamma }}}_{\mathbf{y}°\mathbf{x}}={\widehat{\mathbf{\Gamma }}}_{*}^{\prime }{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1}{\widehat{\mathbf{\Gamma }}}_{*}$ be the sample cross-covariance matrix of the fitted vectors, resulting from a vector autoregression of ${\mathbf{y}}_t$ on ${\mathbf{x}}_t$. The sample canonical correlation matrix between ${\mathbf{y}}_t$ and ${\mathbf{x}}_t$ is defined as ${\mathbf{C}}_{\mathbf{y},\mathbf{x}}=\mathrm{}{\widehat{\mathbf{\Gamma }}}_0^{-1/2}{\widehat{\mathbf{\Gamma }}}_{*}^{\prime }{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2}$, where ${\mathbf{C}}_{\mathbf{x},\mathbf{y}}={\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{\prime }$. Truncated matrices are denoted by superscripts, for example, ${\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{(d)}$ and ${\widehat{\mathbf{\Gamma }}}_{*}^{(d)}$ are formed by the d eigenvectors associated with the d largest singular values of ${\mathbf{C}}_{\mathbf{y},\mathbf{x}}$ and ${\widehat{\mathbf{\Gamma }}}_{*}$, respectively.

Note that the decomposition $\mathit{\beta }=\mathbf{AB}$ is not unique, as it can be expressed as $\mathit{\beta }=\mathbf{A}{\mathbf{C}}^{-1}\mathbf{C}\mathbf{B}$ for any nonsingular matrix $\mathbf{C}\in {\mathbb{R}}^{q\times q}$. The MLEs of the RRVAR model parameters, obtained by Reinsel and Velu (1998), involve normalization constraints for identifiability, such as $\mathbf{B}{\mathbf{\Gamma }}_{(p)}{\mathbf{B}}^{\prime }={\mathbf{\Lambda }}^2,{\mathbf{A}}^{\prime }\mathit{\rho }\mathbf{A}={\mathbf{I}}_d$, where ${\mathbf{\Lambda }}^2=\text{diag}({\lambda }_1^2,\dots ,{\lambda }_d^2)$ and ${\lambda }_i^2$’s are the eigenvalues of ${\mathit{\rho }}^{1/2}{\mathbf{\Gamma }}_{*}^{\prime }{\mathbf{\Gamma }}_{(p)}^{-1}{\mathbf{\Gamma }}_{*}{\mathit{\rho }}^{1/2}$, for any positive definite matrix $\mathit{\rho }$. However, the parameters of interest $\mathit{\beta }$ and $\mathbf{\Sigma }$ are identifiable, with $\text{span}(\mathbf{A})=\text{span}(\mathit{\beta })$, $\text{span}(\mathbf{B}\prime )=\text{span}(\mathit{\beta }\prime )$. We propose a new framework for ML estimation of $\mathit{\beta }$ and $\mathbf{\Sigma }$ that does not require constraints on $\mathbf{A}$ or $\mathbf{B}$, and works for any decomposition $\mathit{\beta }=\mathbf{AB}$ with $\text{rank}(\mathbf{A})=\text{rank}(\mathbf{B})=d$.

The MLEs of the parameters of the RRVAR model (3), which maximize (4), are associated with the canonical correlations (Anderson 2002), and given by ${\widehat{\mathit{\alpha }}}_{\text{ RRVAR}}=\overline{\mathbf{y}}$ and (5) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{\text{ RRVAR}}={\widehat{\mathbf{\Gamma }}}_0^{1/2}{\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{(d)}{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2}, {\widehat{\Sigma }}_{\text{ RRVAR}}={\widehat{\mathbf{\Gamma }}}_0-{\widehat{\mathit{\beta }}}_{\text{ RRVAR}}{\widehat{\mathbf{\Gamma }}}_{*}\\ ={\widehat{\mathbf{\Gamma }}}_0^{1/2}({\mathbf{I}}_q-{\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{(d)}{\mathbf{C}}_{\mathbf{x},\mathbf{y}}^{(d)}){\widehat{\mathbf{\Gamma }}}_0^{1/2}.\end{array}$$(5)

There are different estimators for the RR parameters $\mathbf{A}$ and $\mathbf{B}$ in the literature based on different constraints on $\mathbf{A}$ and $\mathbf{B}$. All of them can be reproduced by decomposing the rank-d estimator ${\widehat{\mathit{\beta }}}_{\text{ RRVAR}}$ given in (5). The OLSVAR estimators of $\mathit{\beta }$ and $\mathbf{\Sigma }$ can be obtained by replacing the truncated sample canonical correlation matrices ${\mathbf{C}}_{(.)}^{(d)}$ with their untruncated versions, that is, ${\widehat{\mathit{\beta }}}_{\text{OLSVAR}}={\widehat{\mathbf{\Gamma }}}_0^{1/2}{\mathbf{C}}_{\mathbf{y},\mathbf{x}}{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2}$ and ${\widehat{\Sigma }}_{\text{OLSVAR}}={\widehat{\mathbf{\Gamma }}}_0^{1/2}({\mathbf{I}}_q-{\mathbf{C}}_{\mathbf{y},\mathbf{x}}{\mathbf{C}}_{\mathbf{x},\mathbf{y}}){\widehat{\mathbf{\Gamma }}}_0^{1/2}$.

2.2 Envelope VAR Model

The envelope model is a parsimonious model introduced by Cook, Li, and Chiaromonte (2010) for multivariate regression that achieves efficiency gains in estimation and better prediction performance. An envelope model for time series data has been developed by Wang and Ding (2018) in the context of VAR models. Before proceeding further, we introduce the following definitions which will be used in the following sections (see also Cook, Li, and Chiaromonte 2010; Cook, Forzani, and Zhang 2015).

Definition 2.1.

A subspace $\mathcal{R}\subseteq {\mathbb{R}}^q$ is defined to be a reducing subspace of $\mathbf{M}\in {\mathbb{R}}^{q\times q}$ ($\mathcal{R}$ reduces M), if and only if $\mathbf{M}={\mathbf{P}}_{\mathcal{R}}\mathbf{M}{\mathbf{P}}_{\mathcal{R}}+{\mathbf{Q}}_{\mathcal{R}}\mathbf{M}{\mathbf{Q}}_{\mathcal{R}}$.

The envelope methodology is based on the notion of reducing subspace, which is essential to the development of theory and practice in functional analysis (Conway 1990).

Definition 2.2.

Let $\mathbf{M}\in {\mathbb{R}}^{q\times q}$ and $\mathcal{S}\subseteq \text{span}(\mathbf{M})$. Then the M-envelope of $\mathcal{S}$ is defined as the intersection of all reducing subspaces of M that contains $\mathcal{S}$ and is denoted by ${\mathcal{E}}_{\mathbf{M}}(\mathcal{S})$.

Since the intersection of any two reducing subspaces of M is again a reducing subspace of M, Definition 2.2 ensures the existence and uniqueness of envelopes. For ease in notation, we use ${\mathcal{E}}_{\mathbf{M}}(\mathcal{D})$ instead of ${\mathcal{E}}_{\mathbf{M}}\{\text{span}(\mathbf{D})\}$. It can be characterized that the M-envelope is the span of some subset of the eigenvectors of M (Cook, Li, and Chiaromonte 2010).

Let $\mathcal{S}\subseteq {\mathbb{R}}^q$ such that (i) ${\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t|{\mathbf{x}}_t\sim {\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t$, and (ii) $\text{cov}({\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t,{\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t|{\mathbf{x}}_t)=0$, where ${\mathbf{P}}_{\mathcal{S}}$ is the projection matrix onto $\mathcal{S}$, and ${\mathbf{Q}}_{\mathcal{S}}$ is the orthogonal projection. ${\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t$ and ${\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t$ are called material and immaterial parts of ${\mathbf{y}}_t$ respectively. Conditions (i) and (ii) together imply that any dependence of ${\mathbf{y}}_t$ on its lagged values ${\mathbf{x}}_t$ must be concentrated in the material part of the autoregression, that is, ${\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t$, while ${\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t$ is invariant to changes in ${\mathbf{x}}_t$ and becomes white noise and immaterial to the estimation of $\mathit{\beta }$. Therefore, ${\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t$ carry all material information in ${\mathbf{y}}_t$, and a change in ${\mathbf{x}}_t$ can affect the distribution of ${\mathbf{y}}_t$ only via ${\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t$. If ${\mathit{\epsilon }}_t$ is normally distributed, then condition (ii) implies that ${\mathbf{P}}_{\mathcal{S}}{\mathbf{y}}_t$ is independent of ${\mathbf{Q}}_{\mathcal{S}}{\mathbf{y}}_t$ given ${\mathbf{x}}_t$. Under the envelope model, it can be shown that (i) and (ii) hold if and only if (a) $\text{span}(\mathit{\beta })=\mathcal{B}\subseteq \mathcal{S}$, and (b) $\mathcal{S}$ is a reducing subspace of $\mathbf{\Sigma }$, that is, $\mathbf{\Sigma }={\mathbf{P}}_{\mathcal{S}}\mathbf{\Sigma }{\mathbf{P}}_{\mathcal{S}}+{\mathbf{Q}}_{\mathcal{S}}\mathbf{\Sigma }{\mathbf{Q}}_{\mathcal{S}}$ (Cook, Li, and Chiaromonte 2010). That is, the $\mathcal{S}$ is a reducing subspace of $\mathbf{\Sigma }$ that contains $\mathcal{B}$, and the intersection of all such subspaces is called $\mathbf{\Sigma }$-envelope of $\mathcal{B}$, and is denoted as ${\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$. Let $u=\dim ({\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B}))$, then under the RRVAR model (3), ${\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})={\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{A})$, where $\mathcal{A}=\text{span}(\mathbf{A})$ and $u\ge d$. This is because $\dim ({\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B}))\ge \dim (\mathcal{B})=\text{rank}(\mathit{\beta })=d$.

Let $\mathbf{\Phi }\in {\mathbb{R}}^{q\times u}$ be an orthogonal basis of ${\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$ with an orthogonal complement basis ${\mathbf{\Phi }}_0\mathrm{}\mathrm{}\in \mathrm{}{\mathbb{R}}^{q\times (q-u)}$ such that $(\mathbf{\Phi },{\mathbf{\Phi }}_0)$ is an orthogonal matrix. Since $\text{span}(\mathit{\beta })\subseteq {\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})=\text{span}(\mathbf{\Phi })$, there exists a $\mathit{\xi }\mathrm{}\in \mathrm{}{\mathbb{R}}^{u\times qp}$ such that $\mathit{\beta }=\mathbf{\Phi }\mathit{\xi }$. Since ${\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$ is spanned by a subset of $\mathbf{\Sigma }$’s eigenvectors, there exist symmetric positive definite matrices $\mathbf{\Omega }\in {\mathbb{R}}^{u\times u}$ and ${\mathbf{\Omega }}_0\in {\mathbb{R}}^{(q-u)\times (q-u)}$ such that $\mathbf{\Sigma }=\mathbf{\Phi }\mathbf{\Omega }{\mathbf{\Phi }}^{\prime }+{\mathbf{\Phi }}_0{\mathbf{\Omega }}_0{\mathbf{\Phi }}_0^{\prime }$. That is, $\mathit{\xi }\in {\mathbb{R}}^{u\times qp}$ contains the coordinates of $\mathit{\beta }$ in terms of the basis matrix $\mathbf{\Phi }$; and $\mathbf{\Omega }$ and ${\mathbf{\Omega }}_0$ carry the coordinates of $\Sigma$ relative to $\mathbf{\Phi }$ and ${\mathbf{\Phi }}_0$, respectively. Thus, the u-dimensional envelope model is summarized as (6) $$\mathit{\beta }=\mathbf{\Phi }\mathit{\xi }, \mathbf{\Sigma }=\mathbf{\Phi }\mathbf{\Omega }{\mathbf{\Phi }}^{\prime }+{\mathbf{\Phi }}_0{\mathbf{\Omega }}_0{\mathbf{\Phi }}_0^{\prime }.$$(6)

This provides a link between $\mathit{\beta }$ and $\mathbf{\Sigma }$. That is, the white noise variation is decomposed into the variation related to the material part of ${\mathbf{y}}_t$, that is, $\mathbf{\Phi }\mathbf{\Omega }{\mathbf{\Phi }}^{\prime }$ = var $({\mathbf{P}}_{\mathbf{\Phi }}{\mathbf{y}}_t)$, and the variation related to the immaterial part of ${\mathbf{y}}_t$, that is, ${\mathbf{\Phi }}_0{\mathbf{\Omega }}_0{\mathbf{\Phi }}_0^{\prime }$ = var $({\mathbf{P}}_{{\mathbf{\Phi }}_0}{\mathbf{y}}_t)$. These two are respectively material and immaterial to the estimation of $\mathit{\beta }$.

Lam and Tam (2011) and Lam and Yao (2012) developed an approach that uses information from lagged autocovariance matrices via eigendecomposition to estimate the factor loading space. This approach has recently been extended to VAR models by Cubadda and Hecq (2022a, 2022b), and Wang et al. (2022b). Although it follows that the orthogonal linear combinations of the factors convey all the relevant information for the estimation (like ${\mathbf{\Phi }}^{\prime }{\mathbf{y}}_t$ in the EVAR model), and projection onto its orthogonal complement (like ${\mathbf{\Phi }}_0^{\prime }{\mathbf{y}}_t$ in the EVAR model) gives white noise, it is distinctly different from the EVAR model as there is no required connection between $\mathit{\beta }$ and the error covariance matrix $\mathbf{\Sigma }$. The main difference between the two approaches is that the envelope methodology exploits connections between $\mathit{\beta }$ and $\mathbf{\Sigma }$ using a minimal reducing subspace, and the information of the basis matrix $\mathbf{\Phi }$ comes from both the conditional mean and the structure of covariance matrix.

2.3 Reduced-Rank Envelope VAR Model

We propose a new parsimonious model, called the reduced-rank envelope vector autoregressive model (REVAR) model, by incorporating the idea of envelope VAR models into the reduced-rank VAR model. This model improves the accuracy and efficiency of VAR estimation. Given that $\text{span}(\mathbf{\Phi })={\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$ with $\dim ({\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B}))=u\le q$, the RRVAR model (3) can be reparameterized into an envelope structure as (7) $$\mathit{\beta }=\mathbf{AB}=\mathbf{\Phi }\mathit{\xi }=\mathbf{\Phi }\mathit{\nu }\mathbf{B}, \mathbf{\Sigma }=\mathbf{\Phi }\mathbf{\Omega }{\mathbf{\Phi }}^{\prime }+{\mathbf{\Phi }}_0{\mathbf{\Omega }}_0{\mathbf{\Phi }}_0^{\prime },$$(7) where the parameterization $\mathit{\beta }=\mathbf{\Phi }\mathit{\xi }$ with $\mathit{\xi }\in {\mathbb{R}}^{u\times qp}$ forms the EVAR model (6) proposed by Wang and Ding (2018), and the columns of $\mathit{\nu }\in {\mathbb{R}}^{u\times d},u\ge d$, are the coordinates of A with respect to the basis $\mathbf{\Phi }$. Then the REVAR model is summarized as follows (8) $${\mathbf{y}}_t=\mathit{\alpha }+\mathbf{\Phi }\mathit{\nu }\mathbf{B}{\mathbf{x}}_t+{\mathit{\epsilon }}_t, \mathbf{\Sigma }=\mathbf{\Phi }\mathbf{\Omega }{\mathbf{\Phi }}^{\prime }+{\mathbf{\Phi }}_0{\mathbf{\Omega }}_0{\mathbf{\Phi }}_0^{\prime },$$(8) where $\mathbf{\Omega },{\mathbf{\Omega }}_0$ are the same as those in (6). There are no further constraints on $\mathbf{A},\mathbf{B},\text{or} \mathit{\nu }$ in this model except that all three matrices have rank d. Notice that the constrained parameters $\mathbf{A},\mathbf{B},\mathbf{\Phi },{\mathbf{\Phi }}_0,\mathit{\xi },\mathit{\nu },\mathbf{\Omega },\text{and} {\mathbf{\Omega }}_0$ in RRVAR, EVAR, and REVAR models are not unique, however, $\mathit{\beta }$ and $\mathbf{\Sigma }$ in (7) are unique. In order to compare the properties of the REVAR model with the RRVAR model, we assume $\mathbf{\Phi }$ is known.

Lemma 1.

Suppose $\mathbf{\Phi }$ is known, then the maximum likelihood estimators (MLEs) of the parameters of the REVAR model (8), which maximize (4) are as follows, $$\begin{array}{cc}{\widehat{\mathit{\alpha }}}_{\mathbf{\Phi }}\hfill & =\overline{\mathbf{y}}, {\widehat{\mathit{\beta }}}_{\mathbf{\Phi }}=\mathbf{\Phi }{\widehat{\mathit{\nu }}}_{\mathbf{\Phi }}{\stackrel{\mathbf{⁁}}{\mathbf{B}}}_{\mathbf{\Phi }}=\mathbf{\Phi }{\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2}{\mathbf{C}}_{\mathbf{\Phi }\prime \mathbf{y},\mathbf{x}}^{(d)}{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2}\hfill \\ {\widehat{\Sigma }}_{\mathbf{\Phi }}\hfill & =\mathbf{\Phi }{\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2}\left({\mathbf{I}}_u-{\mathbf{C}}_{\mathbf{\Phi }\prime \mathbf{y},\mathbf{x}}^{(d)}{\mathbf{C}}_{\mathbf{x},\mathbf{\Phi }\prime \mathbf{y}}^{(d)}\right){\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2}\mathbf{\Phi }\prime +{\mathbf{Q}}_{\mathbf{\Phi }}{\widehat{\mathbf{\Gamma }}}_0{\mathbf{Q}}_{\mathbf{\Phi }}.\hfill \end{array}$$

Note that the MLEs obtained in Lemma 1 are based on $\mathbf{\Phi }$, however, the final estimators ${\widehat{\mathit{\beta }}}_{\mathbf{\Phi }}$ and ${\widehat{\Sigma }}_{\mathbf{\Phi }}$ depend on $\mathbf{\Phi }$ only through $\text{span}(\mathbf{\Phi })={\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$. That is, for any orthogonal matrix $\mathbf{O}\in {\mathbb{R}}^{u\times u}$ we have ${\widehat{\mathit{\beta }}}_{\mathbf{\Phi }}={\widehat{\mathit{\beta }}}_{\mathbf{\Phi }\mathbf{O}}$ and ${\widehat{\Sigma }}_{\mathbf{\Phi }}={\widehat{\Sigma }}_{\mathbf{\Phi }\mathbf{O}}$.

As a result, when the envelope is known, we should concentrate on the reduced response time series $\mathbf{\Phi }\prime {\mathbf{y}}_t$ and obtain the rank-d RRVAR estimator ${\widehat{\mathit{\nu }}}_{\mathbf{\Phi }}{\stackrel{\mathbf{⁁}}{\mathbf{B}}}_{\mathbf{\Phi }}$ of the vector autoregression of $\mathbf{\Phi }\prime {\mathbf{y}}_t$ on ${\mathbf{x}}_t$. Since ${\widehat{\Sigma }}_{\mathbf{\Phi }}={\mathbf{P}}_{\mathbf{\Phi }}{\widehat{\Sigma }}_{\mathbf{\Phi }}{\mathbf{P}}_{\mathbf{\Phi }}+{\mathbf{Q}}_{\mathbf{\Phi }}{\widehat{\Sigma }}_{\mathbf{\Phi }}{\mathbf{Q}}_{\mathbf{\Phi }}$, therefore, according to Definition 2.1, the estimator ${\widehat{\Sigma }}_{\mathbf{\Phi }}$ of ${\Sigma }_{\mathbf{\Phi }}$ is now reduced by $\text{span}(\mathbf{\Phi })$. Thus, $\text{span}(\mathbf{\Phi })$ is a reducing subspace of ${\widehat{\Sigma }}_{\mathbf{\Phi }}$ that also includes $\text{span}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi }})$, that is, the envelope subspace is maintained by the structure of the given ML estimators of the coefficient matrix. We obtain the ML estimator $\widehat{\Phi }$ and show that the REVAR estimators for the coefficient matrix $\mathit{\beta }$ and covariance matrix $\mathbf{\Sigma }$ are identical to the estimators in Lemma 1 by substituting $\mathbf{\Phi }$ with $\widehat{\mathbf{\Phi }}$.

In comparing the models, notice that if the envelope dimension u is equal to the response dimension, that is, $d<u=q$, then there is no immaterial information to be removed by the envelope VAR model and the REVAR model collapses to the RRVAR model (3), that is, $\mathbf{\Phi }={\mathbf{I}}_q$ and $\mathbf{A}=\mathit{\nu }$. If d = u, then the REVAR model is equivalent to the EVAR model. When d = u = q, the REVAR model collapses to the standard VAR model. There is an extreme situation when $p>q=1$, then both the RRVAR and EVAR models collapse to the standard univariate time series model. In this situation, there would not be any reduction. If $q>p=1$, then both the RRVAR and REVAR models gain efficiency. Extreme cases in the REVAR model are numerically explored in Supplement S8.4 to support these assertions.

3 Estimation

3.1 Parameterization for Each Model

We define the following estimable functions for each model. Let h represent the parameter vector of the standard VAR model (1), and let $\mathit{\psi }$, $\mathit{\delta }$, and $\mathit{\theta }$ denote the parameter vectors of the RRVAR, EVAR, and REVAR models, respectively. Since the estimator of the common parameter $\mathit{\alpha }$ is $\widehat{\mathit{\alpha }}=\overline{\mathbf{y}}$ for all methods, and $\overline{\mathbf{y}}$ is asymptotically independent of the estimators of $\mathit{\beta }$ and $\mathbf{\Sigma }$; therefore, $\mathit{\alpha }$ is omitted from all models. Hence, we have (9) $$\begin{array}{c}\mathbf{h}=\left(\begin{array}{c}\text{vec}(\mathit{\beta })\\ \text{vech}(\mathbf{\Sigma })\end{array}\right), \mathit{\psi }=\left(\begin{array}{c}\text{vec}(\mathbf{A})\\ \text{vec}(\mathbf{B})\\ \text{vech}(\mathbf{\Sigma })\end{array}\right),\\ \mathit{\delta }=\left(\begin{array}{c}\text{vec}(\mathbf{\Phi })\\ \text{vec}(\mathit{\xi })\\ \text{vech}(\mathbf{\Omega })\\ \text{vech}({\mathbf{\Omega }}_0)\end{array}\right),\mathit{\theta }=\left(\begin{array}{c}\text{vec}(\mathbf{\Phi })\\ \text{vec}(\mathit{\nu })\\ \text{vec}(\mathbf{B})\\ \text{vech}(\mathbf{\Omega })\\ \text{vech}({\mathbf{\Omega }}_0)\end{array}\right),\end{array}$$(9) where $\mathbf{h}=({\mathbf{h}}_1^{\prime },{\mathbf{h}}_2^{\prime })\prime$, $\mathit{\psi }=({\mathit{\psi }}_1^{\prime },{\mathit{\psi }}_2^{\prime },{\mathit{\psi }}_3^{\prime })\prime ,\mathit{\delta }=({\mathit{\delta }}_1^{\prime },\dots ,{\mathit{\delta }}_4^{\prime })\prime$, and $\mathit{\theta }=({\mathit{\theta }}_1^{\prime },\dots ,{\mathit{\theta }}_5^{\prime })\prime$, respectively. Note that $\mathbf{h}=\mathbf{h}(\mathit{\psi })$ for the RRVAR model, $\mathbf{h}=\mathbf{h}(\mathit{\delta })$ for the EVAR model, and $\mathbf{h}=\mathbf{h}(\mathit{\theta })$ for the REVAR model. To compare the models, we use $\mathcal{T}(·)$ to represent the total number of parameters (NOP) in $\mathbf{h},\mathit{\psi },\mathit{\delta }$, and $\mathit{\theta }$. Then, the parameter count for each model is

1. Standard VAR(p) model, ${\mathcal{T}}_{\text{OLSVAR}}=\mathcal{T}(\mathbf{h})=q^{2}p+q(q+1)/2$;

2. Reduced-rank VAR(p) model, ${\mathcal{T}}_{\text{ RRVAR}}=\mathcal{T}(\mathit{\psi })=d\{q(p+1)-d\}+q(q+1)/2$;

3. Envelope VAR(p) model, ${\mathcal{T}}_{\text{EVAR}}=\mathcal{T}(\mathit{\delta })=\mathit{uqp}+q(q+1)/2$;

4. Reduced-rank envelope VAR(p) model, ${\mathcal{T}}_{\text{REVAR}}=\mathcal{T}(\mathit{\theta })=d\{qp+u-d\}+q(q+1)/2$.

The reduced NOP from the standard VAR model to the RRVAR model is $\mathcal{T}(\mathbf{h})-\mathcal{T}(\mathit{\psi })=(qp-d)(q-d)\ge 0$, and that is even further reduced from the RRVAR model to the REVAR model by the number of $\mathcal{T}(\mathit{\psi })-\mathcal{T}(\mathit{\theta })=d(q-u)\ge 0$. Moreover, compared to the standard VAR model, the reduced NOP by the EVAR model is $\mathcal{T}(\mathbf{h})-\mathcal{T}(\mathit{\delta })=qp(q-u)\ge 0$, and it is further reduced by $\mathcal{T}(\mathit{\delta })-\mathcal{T}(\mathit{\theta })=(qp-d)(u-d)\ge 0$ from the EVAR model to the REVAR model. It is important to note that envelope models achieve efficiency gains not only from parsimony but mainly from the structure of the covariance matrix. To evaluate the impact of covariance matrix structure and variations in the immaterial and material matrices on the performance of envelope models, we conducted simulations by adjusting the ratios of immaterial variation to material variation, as presented in Supplement S8.5.

3.2 Maximum Likelihood Estimation for the REVAR Model

In this section, we derive the ML estimators of the REVAR model (8) for a given lag order p, rank d, and envelope dimension u. How to select p, d, and u are discussed in Section 5. The MLEs for the REVAR model can be obtained by replacing h with $\mathbf{h}(\mathit{\theta })$ in the log-likelihood function (4) as $L_T(\mathit{\alpha },\mathit{\beta }(\mathit{\theta }),\mathbf{\Sigma }(\mathit{\theta }))\equiv L_T(\mathit{\alpha },\mathbf{B},\mathit{\nu },\mathbf{\Omega },{\mathbf{\Omega }}_0,\mathbf{\Phi }|p,d,u)$ and maximizing it with respect to all parameters other than $\mathbf{\Phi }$. Since $\text{span}(\mathbf{\Phi })$ lies in a Grassmann manifold (Cook, Li, and Chiaromonte 2010), therefore, we cannot analytically find a unique optimal value of $\mathbf{\Phi }$. As it is discussed in Proposition 1, we can obtain $\widehat{\mathbf{\Phi }}$ from minimization over a Grassmannian.

To obtain the MLE of $\mathbf{\Phi }$, let $\mathbf{D}\in {\mathbb{R}}^{q\times u}$ be a semiorthogonal matrix and define the standardized version of $\mathbf{D}\prime {\mathbf{y}}_t\in {\mathbb{R}}^u$ as ${\mathbf{z}}_{\mathbf{D}}={(\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0\mathbf{D})}^{-1/2}\mathbf{D}\prime {\mathbf{y}}_t$ with sample covariance matrix ${\mathbf{I}}_u$. Suppose ${\widehat{\omega }}_i(\mathbf{D}),i=1,\dots ,u$, is the ith eigenvalue of ${\widehat{\mathbf{\Gamma }}}_{{\mathbf{z}}_{\mathbf{D}}|\mathbf{x}}^{-1}$ which is also the ith eigenvalue of ${(\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}\mathbf{D})}^{-1/2}(\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0\mathbf{D})$ ${(\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}\mathbf{D})}^{-1/2}$.

Proposition 1.

Let $L_T(\mathit{\alpha },\mathbf{B},\mathit{\nu },\mathbf{\Omega },{\mathbf{\Omega }}_0,\mathbf{\Phi }|p,d,u)$ denotes the log-likelihood function of the REVAR model (8). Then, it achieves its maximum at $\widehat{\mathbf{\Phi }}=\arg {\min }_{\mathbf{D}}{\mathbf{F}}_T(\mathbf{D}|p,d,u)$, where (10) $${\mathbf{F}}_T(\mathbf{D}|p,d,u)=\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0^{-1}\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|{\mathbf{I}}_u-{\widehat{\mathbf{\Gamma }}}_{{\mathbf{z}}_{\mathbf{D}}°\mathbf{x}}^{(d)}|$$(10) (11) $$=\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0^{-1}\mathbf{D}|+\sum _{i=d+1}^u\hspace{0.17em}\log \hspace{0.17em}\{{\widehat{\omega }}_i(\mathbf{D})\},$$(11) where the minimization is over the Grassmannian of dimension $u\in {\mathbb{R}}^q$, that is, ${\mathcal{G}}_{q,u}$.

The expansion of the objective function in (10) offers a method of interpreting the log-likelihood functions for envelope-based models. For instance, the objective function of the EVAR model (6) (Wang and Ding 2018) can be similarly expanded as follows (12) $$\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\widehat{\mathbf{\Gamma }}}_0^{-1}\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|{\mathbf{I}}_u-{\widehat{\mathbf{\Gamma }}}_{{\mathbf{z}}_{\mathbf{D}}°\mathbf{x}}|,$$(12) which is similar to the expression (10) except that the last term in (12) is based on the ordinary least squares method, whereas the last term in (10) is based on the reduced-rank VAR model. In other words, the two expressions in (10) and (12) are the same when d = u.

Proposition 2.

The sample objective function ${\mathbf{F}}_T(\mathbf{D}|p,d,u)$ in (11) converges in probability to its population counterpart $\mathbf{F}(\mathbf{D}|p,u)=\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime \mathbf{\Sigma }\mathbf{D}|+\hspace{0.17em}\log \hspace{0.17em}|\mathbf{D}\prime {\mathbf{\Gamma }}_0^{-1}\mathbf{D}|$ uniformly in D, as $T\to \infty$. The estimator $\widehat{\mathbf{\Phi }}=\arg {\min }_{\mathbf{D}}{\mathbf{F}}_T(\mathbf{D}|p,d,u)$ is consistent, and ${\widehat{\mathcal{E}}}_{\mathbf{\Sigma }}(\mathcal{B})=\text{span}\{\arg {\min }_{\mathbf{D}}{\mathbf{F}}_T(\mathbf{D}|p,u)\}$.

The following Proposition provides the MLEs of the REVAR model parameters in (8).

Proposition 3.

The ML estimators for the REVAR model (8) that maximize the conditional log-likelihood function in (4) are ${\widehat{\mathit{\alpha }}}_{\text{REVAR}}=\overline{\mathbf{y}}$, $\widehat{\mathbf{\Phi }}=\arg {\min }_{\mathbf{D}}{\mathbf{F}}_T(\mathbf{D}|p,d,u)$, and $$\begin{array}{c}{\widehat{\mathbf{\Omega }}}_0={\widehat{\mathbf{\Phi }}}_0^{\prime }{\widehat{\mathbf{\Gamma }}}_0{\widehat{\mathbf{\Phi }}}_0, {\widehat{\Omega }}_{}={\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2}\left({\mathbf{I}}_u-{\mathbf{C}}_{\mathbf{\Phi }\prime \mathbf{y},\mathbf{x}}^{(d)}{\mathbf{C}}_{\mathbf{x},\mathbf{\Phi }\prime \mathbf{y}}^{(d)}\right){\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2},\\ {\widehat{\mathit{\beta }}}_{\text{REVAR}}=\widehat{\mathbf{\Phi }}\widehat{\mathit{\nu }}{\widehat{\mathbf{B}}}_{\text{REVAR}}=\widehat{\mathbf{\Phi }}{\widehat{\mathbf{\Gamma }}}_{\mathbf{\Phi }\prime \mathbf{y}}^{1/2}{\mathbf{C}}_{\mathbf{\Phi }\prime \mathbf{y},\mathbf{x}}^{(d)}{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2},\\ {\widehat{\mathbf{\Sigma }}}_{\text{REVAR}}=\widehat{\mathbf{\Phi }}\widehat{\mathbf{\Omega }}{\widehat{\mathbf{\Phi }}}^{\prime }+{\widehat{\mathbf{\Phi }}}_0{\widehat{\mathbf{\Omega }}}_0{\widehat{\mathbf{\Phi }}}_0^{\prime }.\end{array}$$

Note that rank $({\widehat{\mathit{\beta }}}_{\text{REVAR}})=d (\le u)$, and $\text{span}({\widehat{\mathit{\beta }}}_{\text{REVAR}})\subseteq {\mathcal{E}}_{\mathbf{\Sigma }}(\mathcal{B})$. In comparison to the RRVAR model, ${\widehat{\Sigma }}_{\text{REVAR}}={\mathbf{P}}_{\widehat{\mathbf{\Phi }}}{\widehat{\Sigma }}_{\text{REVAR}}{\mathbf{P}}_{\widehat{\mathbf{\Phi }}}+{\mathbf{Q}}_{\widehat{\mathbf{\Phi }}}{\widehat{\Sigma }}_{\text{REVAR}}{\mathbf{Q}}_{\widehat{\mathbf{\Phi }}}$ now has an envelope structure. If u = q, then $\mathbf{\Phi }={\mathbf{I}}_q,\mathbf{A}=\mathit{\nu }$, and the REVAR model is reduced to the RRVAR model.

The full Grassmann (FG) optimization method can be computationally slow and expensive, particularly in high-dimensional problems. Therefore, we used the FG algorithm for the small VARs, but for higher dimensional cases, we adapt and employ the one-dimensional (1D) algorithm proposed by Cook and Zhang (2016). The 1D algorithm is computationally more efficient and robust compared to the FG optimization. This is due to the fact that it decomposes the u-dimensional FG optimization into a series of u one-dimensional optimization problems and does not require an initial guess (see Cook and Zhang 2018).

4 Asymptotic Properties

In this section, we establish the asymptotic properties of the proposed ML estimators assuming both Gaussian and non-Gaussian white noise processes. We derive their asymptotic distributions under the reduced-rank envelope VAR model and compare their asymptotic efficiencies to those of the OLSVAR, RRVAR, and REVAR models. Related asymptotic results for estimators of the parameters of the RRVAR, EVAR, and OLSVAR models can be found in Anderson (2002), Wang and Ding (2018), and Lütkepohl (2005), respectively. The asymptotic comparison between ${\widehat{\mathit{\beta }}}_{\text{ RRVAR}}$ and ${\widehat{\mathit{\beta }}}_{\text{OLSVAR}}$ can be found in Anderson (2002). Therefore, our focus here is on comparing ${\widehat{\mathit{\beta }}}_{\text{REVAR}}$ and ${\widehat{\mathit{\beta }}}_{\text{ RRVAR}}$. The asymptotic results of ${\widehat{\mathit{\beta }}}_{\text{REVAR}}$ over ${\widehat{\mathit{\beta }}}_{\text{EVAR}}$ are comparable to those of ${\widehat{\mathit{\beta }}}_{\text{ RRVAR}}$ over ${\widehat{\mathit{\beta }}}_{\text{OLSVAR}}$. This is due to the rank reduction constraint present in the material response time series $\mathbf{\Phi }\prime {\mathbf{y}}_t$.

Let ${\mathbf{J}}_{\mathbf{h}}$ be the Fisher information of $\mathbf{h}=({\mathbf{h}}_1^{\prime },{\mathbf{h}}_2^{\prime })\prime =(\text{vec}\prime (\mathit{\beta }),\text{vech}\prime (\mathbf{\Sigma }))\prime$. Then, the asymptotic covariance matrix of the OLSVAR estimator of h is given by Lütkepohl (2005) (13) $$\begin{array}{c}\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{OLSVAR}}\right)={\mathbf{J}}_{\mathbf{h}}^{-1}=\left(\begin{array}{cc}{\mathbf{J}}_{\mathit{\beta }}^{-1}& \mathbf{0}\\ \mathbf{0}& {\mathbf{J}}_{\mathbf{\Sigma }}^{-1}\end{array}\right)\\ =\left(\begin{array}{cc}{\mathbf{\Gamma }}_{(p)}^{-1}\otimes \mathbf{\Sigma }& \mathbf{0}\\ \mathbf{0}& 2{\mathbf{E}}_q^{†}(\mathbf{\Sigma }\otimes \mathbf{\Sigma }){\mathbf{E}}_q^{†}\prime \end{array}\right).\end{array}$$(13)

This is also the asymptotic covariance matrix of the unrestricted ML estimator of h. Let us define the gradient matrices of $\mathbf{h}(\mathit{\psi })$ and $\mathbf{h}(\mathit{\theta })$ as (14) $$\begin{array}{c}\mathbf{H}=\frac{\partial \mathbf{h}(\mathit{\psi })}{\partial \mathit{\psi }}=\left(\begin{array}{ccc}\frac{\partial {\mathbf{h}}_1}{\partial {\mathit{\psi }}_1^{\prime }}& \dots & \frac{\partial {\mathbf{h}}_1}{\partial {\mathit{\psi }}_3^{\prime }}\\ \frac{\partial {\mathbf{h}}_2}{\partial {\mathit{\psi }}_1^{\prime }}& \dots & \frac{\partial {\mathbf{h}}_2}{\partial {\mathit{\psi }}_3^{\prime }}\end{array}\right),\\ \mathbf{R}=\frac{\partial \mathbf{h}(\mathit{\theta })}{\partial \mathit{\theta }}=\left(\begin{array}{ccc}\frac{\partial {\mathbf{h}}_1}{\partial {\mathit{\theta }}_1^{\prime }}& \dots & \frac{\partial {\mathbf{h}}_1}{\partial {\mathit{\theta }}_5^{\prime }}\\ \frac{\partial {\mathbf{h}}_2}{\partial {\mathit{\theta }}_1^{\prime }}& \dots & \frac{\partial {\mathbf{h}}_2}{\partial {\mathit{\theta }}_5^{\prime }}\end{array}\right).\end{array}$$(14)

The asymptotic covariance matrices of the RRVAR and REVAR estimators, that is, ${\widehat{\mathbf{h}}}_{\text{ RRVAR}}=\mathbf{h}(\widehat{\mathit{\psi }})$ and ${\widehat{\mathbf{h}}}_{\text{REVAR}}=\mathbf{h}(\widehat{\mathit{\theta }})$, are provided in Proposition 4, which can be obtained by using the asymptotic theory of overparameterized structural models proposed by Shapiro (1986).

Proposition 4.

Suppose ${\mathit{\epsilon }}_t\sim \mathcal{N}(\mathbf{0},\mathbf{\Sigma })$, then it can be shown that $\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{OLSVAR}}\right)={\mathbf{J}}_{\mathbf{h}}^{-1}$, $\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{ RRVAR}}\right)=\mathbf{H}{(\mathbf{H}\prime {\mathbf{J}}_{\mathbf{h}}\mathbf{H})}^{†}\mathbf{H}\prime$, and $\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{REVAR}}\right)=\mathbf{R}{(\mathbf{R}\prime {\mathbf{J}}_{\mathbf{h}}\mathbf{R})}^{†}\mathbf{R}\prime$. The differences between the asymptotic covariances are as follows $$\begin{array}{c}\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{OLSVAR}}\right)-\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{ RRVAR}}\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}={\mathbf{J}}_{\mathbf{h}}^{-1/2}{\mathbf{Q}}_{{\mathbf{J}}_{\mathbf{h}}^{1/2}\mathbf{H}}{\mathbf{J}}_{\mathbf{h}}^{-1/2}\ge 0,\\ \text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{ RRVAR}}\right)-\text{avar}\left(\sqrt{T}{\widehat{\mathbf{h}}}_{\text{REVAR}}\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}={\mathbf{J}}_{\mathbf{h}}^{-1/2}\left({\mathbf{P}}_{{\mathbf{J}}_{\mathbf{h}}^{1/2}\mathbf{H}}-{\mathbf{P}}_{{\mathbf{J}}_{\mathbf{h}}^{1/2}\mathbf{R}}\right){\mathbf{J}}_{\mathbf{h}}^{-1/2}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}={\mathbf{J}}_{\mathbf{h}}^{-1/2}{\mathbf{P}}_{{\mathbf{J}}_{\mathbf{h}}^{1/2}\mathbf{H}}{\mathbf{Q}}_{{\mathbf{J}}_{\mathbf{h}}^{1/2}\mathbf{R}}{\mathbf{J}}_{\mathbf{h}}^{-1/2}\ge 0,\end{array}$$

where $†$ denotes the Moore-Penrose inverse. Moreover, we have $\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}})\}\mathrm{}\ge \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{ RRVAR}})\}\ge \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}})\}$. Similarly in compared to the envelope VAR estimator, $\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}})\}\ge \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{EVAR}})\}\ge \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}})\}$.

Proposition 5.

Suppose ${\mathit{\epsilon }}_t\sim \mathcal{N}(\mathbf{0},\mathbf{\Sigma })$, and the rank of the coefficient matrix $\mathit{\beta }$ is d. Then, $\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}}-\mathit{\beta })$ and $\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{ RRVAR}}-\mathit{\beta })$ asymptotically follow normal distributions with mean zero and the following covariance matrices (15) $$\begin{array}{c}\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}})\}\\ ={\mathbf{\Gamma }}_{(p)}^{-1}\otimes \mathbf{\Sigma },\\ \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{ RRVAR}})\}\\ =({\mathbf{I}}_{q^{2}p}-{\mathbf{Q}}_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})}\otimes {\mathbf{Q}}_{\mathbf{A}({\mathbf{\Sigma }}^{-1})})\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}})\}\end{array}$$(15) (16) $$=\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{A}}\mathbf{Q}{\prime }_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})})\}+\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{B}})\},$$(16) where $\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{A}}\mathbf{Q}{\prime }_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})})\}=({\mathbf{Q}}_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})}{\mathbf{\Gamma }}_{(p)}^{-1})\otimes ({\mathbf{P}}_{\mathbf{A}({\mathbf{\Sigma }}^{-1})}$ $\mathbf{\Sigma })$, and $\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{B}})\}=({\mathbf{P}}_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})}{\mathbf{\Gamma }}_{(p)}^{-1})\otimes \mathbf{\Sigma }$.

The asymptotic variance in (15) follows from Anderson (2002, eq. 5.22(5) $$\begin{array}{c}{\widehat{\mathit{\beta }}}_{\text{ RRVAR}}={\widehat{\mathbf{\Gamma }}}_0^{1/2}{\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{(d)}{\widehat{\mathbf{\Gamma }}}_{(p)}^{-1/2}, {\widehat{\Sigma }}_{\text{ RRVAR}}={\widehat{\mathbf{\Gamma }}}_0-{\widehat{\mathit{\beta }}}_{\text{ RRVAR}}{\widehat{\mathbf{\Gamma }}}_{*}\\ ={\widehat{\mathbf{\Gamma }}}_0^{1/2}({\mathbf{I}}_q-{\mathbf{C}}_{\mathbf{y},\mathbf{x}}^{(d)}{\mathbf{C}}_{\mathbf{x},\mathbf{y}}^{(d)}){\widehat{\mathbf{\Gamma }}}_0^{1/2}.\end{array}$$(5) ). Moreover, the asymptotic results in Proposition 5 depend on the decomposed components $\mathbf{A}$ and $\mathbf{B}$ only through their orthogonal projections ${\mathbf{Q}}_{\mathbf{A}({\mathbf{\Sigma }}^{-1})}$ and ${\mathbf{Q}}_{\mathbf{B}({\mathbf{\Gamma }}_{(p)})}$, respectively. Therefore, these results are valid for any decomposition $\mathit{\beta }=\mathbf{AB}$ that satisfies $\mathbf{A}\in {\mathbb{R}}^{q\times d}$ and $\mathbf{B}\in {\mathbb{R}}^{d\times qp}$.

Proposition 6.

Suppose ${\mathit{\epsilon }}_t\sim \mathcal{N}(\mathbf{0},\mathbf{\Sigma })$, then under the REVAR model $\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}}-\mathit{\beta })$ follows asymptotic normal distribution with zero mean and covariance (17) $$\begin{array}{c}\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}})\}\\ =\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi }})\}+\text{avar}\{\sqrt{T}\text{vec}({\mathbf{Q}}_{\mathbf{\Phi }}{\widehat{\mathit{\beta }}}_{\mathit{\nu },\mathbf{B}})\}\end{array}$$(17) (18) $$\begin{array}{c}=\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi },\mathit{\nu }}\mathbf{Q}{\prime }_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})})\}+\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi },\mathbf{B}})\}\\ +\text{avar}\{\sqrt{T}\text{vec}({\mathbf{Q}}_{\mathbf{\Phi }}{\widehat{\mathit{\beta }}}_{\mathit{\nu },\mathbf{B}})\},\end{array}$$(18) where from (16) we have $\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi },\mathit{\nu }}\mathbf{Q}{\prime }_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})})\}=\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{A}}\mathbf{Q}{\prime }_{\mathbf{B}\prime ({\mathbf{\Gamma }}_{(p)})})\}$.

The asymptotic advantages of the REVAR over the RRVAR model can be obtained by using the results in Propositions 5 and 6. By subtracting (18) from (16), we have (19) $$\begin{array}{cc}\hfill & \text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{ RRVAR}})\}-\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}})\}\hfill \\ \hfill & =\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{B}})\}\hfill \\ \hfill & -\text{avar}\{\sqrt{T}\text{vec}({\widehat{\mathit{\beta }}}_{\mathbf{\Phi },\mathbf{B}})\}-\text{avar}\{\sqrt{T}\text{vec}({\mathbf{Q}}_{\mathbf{\Phi }}{\widehat{\mathit{\beta }}}_{\mathit{\nu },\mathbf{B}})\}\ge 0,\hfill \end{array}$$(19) where for a given $\mathbf{B}\in {\mathbb{R}}^{d\times qp}$, the estimators are defined as ${\widehat{\mathit{\beta }}}_{\mathbf{B}}={\widehat{\mathbf{A}}}_{\mathbf{B}}\mathbf{B},{\widehat{\mathit{\beta }}}_{\mathbf{\Phi },\mathbf{B}}=\mathbf{\Phi }{\widehat{\mathit{\nu }}}_{\mathbf{\Phi },\mathbf{B}} \mathbf{B}={\widehat{\mathbf{A}}}_{\mathbf{\Phi },\mathbf{B}}\mathbf{B}$, and ${\widehat{\mathit{\beta }}}_{\mathit{\nu },\mathbf{B}}={\widehat{\mathbf{\Phi }}}_{\mathit{\nu },\mathbf{B}} \mathit{\nu }\mathbf{B}={\widehat{\mathbf{A}}}_{\mathit{\nu },\mathbf{B}}\mathbf{B}$.

4.1 Asymptotic Properties Under Nonnormality

Suppose ${\widehat{\mathbf{h}}}_{\text{OLSVAR}}=(\text{vec}\prime ({\widehat{\mathit{\beta }}}_{\text{OLSVAR}}),\text{vech}\prime ({\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}))\prime$ is the OLSVAR estimator of h under the unrestricted VAR model, and ${\widehat{\mathbf{h}}}_{\text{REVAR}}=\mathbf{h}(\widehat{\mathit{\psi }})$ denote the REVAR estimator. Moreover, we assume that ${\mathbf{h}}_0$ and ${\mathit{\psi }}_0$ are the true values of h and $\mathit{\psi }$, respectively. The following objective function is obtained after maximizing $L_T(\mathit{\alpha },\mathit{\beta },\mathbf{\Sigma })$ in (4) with respect to $\mathit{\alpha }$ $$\begin{array}{c}{L}_T(\mathit{\beta },\mathbf{\Sigma })\propto \frac{-T}{2}\{\hspace{0.17em}\log \hspace{0.17em}\left|\mathbf{\Sigma }\right|+\text{tr}[{\mathbf{\Sigma }}^{-1}({\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}\\ +({\widehat{\mathit{\beta }}}_{\text{OLSVAR}}-\mathit{\beta }){\widehat{\mathbf{\Gamma }}}_{(p)}({\widehat{\mathit{\beta }}}_{\text{OLSVAR}}-\mathit{\beta })\prime )]\}.\end{array}$$

We now consider $L_T(\mathit{\beta },\mathbf{\Sigma })$ as a function of h and ${\widehat{\mathbf{h}}}_{\text{OLSVAR}}$, and define the discrepancy function as $\mathcal{F}(\mathbf{h},{\widehat{\mathbf{h}}}_{\text{OLSVAR}})=2/T\{L_T({\widehat{\mathit{\beta }}}_{\text{OLSVAR}},{\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}})-L_T(\mathit{\beta },\mathbf{\Sigma })\}$, which satisfies the necessary conditions of Shapiro’s (1986). It can be shown that the second derivative of $\mathcal{F}(\mathbf{h},{\widehat{\mathbf{h}}}_{\text{OLSVAR}})$, that is, ${\mathbf{J}}_{\mathbf{h}}\mathrm{}\mathrm{}=1/2({\partial }^2\mathcal{F}(\mathbf{h},{\widehat{\mathbf{h}}}_{\text{OLSVAR}})/\partial \mathbf{h}\partial \mathbf{h}\prime )$ evaluated at ${\widehat{\mathbf{h}}}_{\text{OLSVAR}}\mathrm{}\mathrm{}={\mathbf{h}}_0$ and $\mathbf{h}={\mathbf{h}}_0$ is the Fisher information matrix for h when ${\mathit{\epsilon }}_t$ is normally distributed. The following proposition gives the asymptotic distribution of ${\widehat{\mathbf{h}}}_{\text{REVAR}}$ without the normality assumption on ${\mathit{\epsilon }}_t$.

Proposition 7.

Assume that the error terms of the reduced-rank envelope VAR model (8) are iid and the fourth moments of ${\mathit{\epsilon }}_t$ are finite. Then, $\sqrt{T}({\widehat{\mathbf{h}}}_{\text{OLSVAR}}-{\mathbf{h}}_0)\underset{}{\overset{\mathcal{D}}{\to }}\mathcal{N}(\mathbf{0},\tilde{\mathbf{V}})$ for some positive definite covariance matrix $\tilde{\mathbf{V}}$, and $\sqrt{T}({\widehat{\mathbf{h}}}_{\text{REVAR}}-{\mathbf{h}}_0)\underset{}{\overset{\mathcal{D}}{\to }}\mathcal{N}(\mathbf{0},\tilde{\mathbf{Z}})$, with $\tilde{\mathbf{Z}}=\mathbf{R}{(\mathbf{R}\prime {\mathbf{J}}_{\mathbf{h}}\mathbf{R})}^{†}\mathbf{R}\prime {\mathbf{J}}_{\mathbf{h}}\tilde{\mathbf{V}}{\mathbf{J}}_{\mathbf{h}}\mathbf{R}{(\mathbf{R}\prime {\mathbf{J}}_{\mathbf{h}}\mathbf{R})}^{†}\mathbf{R}\prime$, where ${\mathbf{J}}_{\mathbf{h}}$ and $\mathbf{R}$ are defined in (13) and (14), respectively. Particularly, $\sqrt{T}(\text{vec}({\widehat{\mathit{\beta }}}_{\text{REVAR}})-\text{vec}(\widehat{\mathit{\beta }}))\underset{}{\overset{\mathcal{D}}{\to }}\mathcal{N}(\mathbf{0},{\tilde{\mathbf{Z}}}_{11})$, where ${\tilde{\mathbf{Z}}}_{11}$ is the top-left block matrix of $\tilde{\mathbf{Z}}$ of dimension $q^{2}p\times q^{2}p$.

The $\sqrt{T}$-consistency of ${\widehat{\mathit{\beta }}}_{\text{REVAR}}$, the reduced-rank envelope VAR estimator, relies on the $\sqrt{T}$-consistency of both ${\widehat{\mathit{\beta }}}_{\text{OLSVAR}}$ and ${\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}$, despite the nonnormality of the error terms, and the properties of the discrepancy function $\mathcal{F}(\mathbf{h},{\widehat{\mathbf{h}}}_{\text{OLSVAR}})$. The asymptotic covariance matrix ${\tilde{\mathbf{Z}}}_{11}$ can be estimated conveniently by plugging in the estimated covariance matrix $\tilde{\mathbf{V}}$ into $\tilde{\mathbf{Z}}$, but its accuracy depends on the distribution of ${\mathit{\epsilon }}_t$ for any fixed sample size.

5 Selections of p, d, and u of VAR Models

5.1 Lag Order (p) Selection

The selection of the lag order (p) is an empirical problem and a critical element in the specification of VAR models. Therefore, the first step in VAR analysis is to determine the lag order. To do this, several model selection criteria are employed. The common approach is to fit VAR(p) models with lag orders $p=0,\dots ,p_{\max }$ and select the value of p that minimizes some model selection information criteria (IC). The general form of IC for VAR models is given by $IC(p)=\hspace{0.17em}\ln \hspace{0.17em}|\widehat{\mathbf{\Sigma }}(p)|+c_T\varphi (p)$ (Lütkepohl 2005; Tsay 2014), where $\widehat{\mathbf{\Sigma }}(p)=T^{-1}{\sum }_{t=1}^T{\widehat{\mathit{\epsilon }}}_t{\widehat{\mathit{\epsilon }}}_t^{\prime }$, which is the estimated residual covariance matrix for a model of order p, c_T is a sequence that depends on the sample size T, $\varphi (p)$ is a penalty function that penalizes large VAR(p) models. The most common information criteria are the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), given as $$\begin{array}{c}AIC(p)=\hspace{0.17em}\ln \hspace{0.17em}|\widehat{\mathbf{\Sigma }}(p)|+\frac{2}{T}{\mathcal{T}}_{\text{OLSVAR}},\\ BIC(p)=\hspace{0.17em}\ln \hspace{0.17em}|\widehat{\mathbf{\Sigma }}(p)|+\frac{\hspace{0.17em}\ln \hspace{0.17em}(T)}{T}{\mathcal{T}}_{\text{OLSVAR}},\end{array}$$ where ${\mathcal{T}}_{\text{OLSVAR}}=q^{2}p+q(q+1)/2$, is the number of parameters in the standard VAR model.

5.2 Rank (d) Selection

The implicit assumption of RR models is that the coefficient matrix is not of full rank. Optimal rank selection has an important role in dimension reduction. In the context of RRR, Bura and Cook (2003) proposed a rank selection test that follows a Chi-squared distribution and only requires the finite second moments of the response variables. The rank of $\mathit{\beta }$ can be determined using the test statistic ${\mathcal{M}}_d=T{\sum }_{j=d+1}^q{\lambda }_j^2$, where ${\lambda }_1\ge \cdots \ge {\lambda }_q$ are the eigenvalues of the qp × q matrix ${\widehat{\mathit{\beta }}}_{\text{std}}=\sqrt{(T-qp-1)/T} {\widehat{\mathbf{\Gamma }}}_{(p)}^{1/2}{\widehat{\mathit{\beta }}}_{\text{OLSVAR}}^{\prime }{\widehat{\mathbf{\Gamma }}}_{\mathbf{y}|\mathbf{x}}^{-1/2}$. Under the null hypothesis H₀: d = d₀, the test statistic ${\mathcal{M}}_{d_0}$ has an asymptotic Chi-square distribution, that is, ${\mathcal{M}}_{d_0}\sim {\mathcal{X}}_{(qp-d_0)(q-d_0)}^2$ (Bura and Cook 2003). To find the rank d, one can compute a sequence of test statistics ${\mathcal{M}}_{d_0}$, for $d_0=0,1,\dots ,q-1$, and compare them to the percentiles of their corresponding null distributions ${\mathcal{X}}_{(qp-d_0)(q-d_0)}^2$. The test procedure stops at the first nonsignificant test of H₀: d = d₀. Then, d₀ is an estimate of the rank of $\mathit{\beta }$.

5.3 Envelope Dimension (u) Selection

To select the envelope dimension u, there are two ways. One approach is to simultaneously determine d and u by seeking (d, u) from (0, 0) to (q, q) ($d\le u\le q$), and selecting the (d, u) pair that has the lowest IC value. The IC to determine the optimal (d, u) pair are $$\begin{array}{c}AIC(d,u)=2{\mathcal{T}}_{\text{REVAR}}-2{\widehat{L}}_{d,u},\\ BIC(d,u)=\hspace{0.17em}\log \hspace{0.17em}(T){\mathcal{T}}_{\text{REVAR}}-2{\widehat{L}}_{d,u},\end{array}$$ where ${\mathcal{T}}_{\text{REVAR}}$ is the number of REVAR model parameters, and ${\widehat{L}}_{d,u}$ is the maximized log-likelihood function calculated at the MLEs in Proposition 3. Alternatively, we can first find d using the method in Section 5.2, then search for the value of u from d to q–1 that minimizes AIC or BIC. In this case, under the null hypothesis $H_0:u=u_0$, the test statistic ${\mathcal{M}}_{d,u_0}=2({\widehat{L}}_{d,q}-{\widehat{L}}_{d,u_0})$ has an asymptotic Chi-square distribution with $(q-u_0)d$ degree of freedom, that is, ${\mathcal{M}}_{d,u_0}\sim {\mathcal{X}}_{(q-u_0)d}^2$. The latter method is computationally more efficient.

6 Simulation Studies

In this section, we compare the performance of our proposed REVAR model with the RRVAR, the EVAR model, and the OLSVAR model using simulation studies under different data generating processes (DGP). We simulate data under various parameter settings from model (8) with ${[\Omega ]}_{ij}={(-0.9)}^{|i-j|}$ and ${[{\Omega }_0]}_{ij}=5{(-0.5)}^{|i-j|}$ for $\mathbf{\Omega }$ and ${\mathbf{\Omega }}_0$, respectively. The semiorthogonal matrix $\Phi$ and matrix $\mathbf{B}$ are generated from the uniform distribution on (0, 1). The entries of $\mathit{\nu }$ and p presample observations (${\mathbf{y}}_0,\dots ,{\mathbf{y}}_{-p+1}$) are generated from the standard normal distribution. Then, $\mathbf{\Phi }$ and $\mathit{\beta }=\mathbf{\Phi }\mathit{\nu }\mathbf{B}$ are standardized so that $\mathbf{\Phi }\prime \mathbf{\Phi }={\mathbf{I}}_u$ and $\left|\right|\mathit{\beta }|{|}_F=1$, where $\left|\right|.|{|}_F$ is the Frobenius norm, and the coefficient matrix $\mathit{\beta }$ satisfies the stationary condition. Estimation errors are obtained by comparing the estimated coefficient matrix $\widehat{\mathit{\beta }}$ to the true coefficient matrix $\mathit{\beta }$ by using $\left|\right|\mathit{\beta }-\widehat{\mathit{\beta }}|{|}_F$. We conduct simulation studies under various error distributions with sample sizes of $T=160, 270, 450, 740, 1200, \text{and} 2000$. Each scenario is replicated 100 times, and the minimum and maximum standard error ratios of coefficient estimates ($r_{\min }=\min (S{E}_{\mathrm{M}}/S{E}_{\text{REVAR}})$ and $r_{\max }=\max (S{E}_{\mathrm{M}}/S{E}_{\text{REVAR}})$) are calculated for the M = OLSVAR, EVAR, and RRVAR models compared to our proposed REVAR model. Matlab codes are available upon request.

6.1 Simulation Studies with Normal Errors

In this section, data with normal errors are analyzed. Table 1 shows the total number of parameters for different VAR models in our simulation studies with various combinations of $(d,u,p,q)$. The REVAR model has the fewest parameters, making it more parsimonious than other models. Consequently, it significantly outperforms them, as shown in Figure 1.

Fig. 1 Average estimation errors of moderate-dimensional VAR models with different combinations of $(d,u,p,q)$ against sample size. Refer to Table 1 for details.

Table 1 Total number of parameters (NOP) in each VAR model correspond to Figure 1.

	Large envelop	Nearly full rank	Typical Scenario	p > 1
$(d,u,p,q)$	$(2,6,1,7)$	$(5,6,1,7)$	$(3,4,1,7)$	$(3,4,2,7)$
OLSVAR	77	77	77	126
RRVAR	52	73	61	82
EVAR	70	70	56	84
REVAR	50	68	52	73

Figure 1 displays the impact of envelope dimension and rank on the relative performance of each method. The simulation results in Figure 1 are obtained using the full Grassmannian (FG) algorithm. In the scenario with a large envelope dimension (Figure 1(a)), the RRVAR approach outperforms the OLSVAR approach, while the EVAR model shows a relatively smaller improvement. The REVAR approach has a slight advantage over RRVAR in this scenario. In the second case (Figure 1(b)), where $\mathit{\beta }$ is almost full rank, the RRVAR method outperforms the OLSVAR by a small margin. Both the EVAR and REVAR approaches demonstrate significant improvement over the RRVAR and OLSVAR, where the REVAR is superior to all others. For the third scenario (Figure 1(c)) neither the EVAR nor RRVAR approaches are favored. Both the RRVAR and EVAR approaches exhibit similar behavior and show substantial improvements compared to the OLSVAR approach. The final scenario (Figure 1(d)) is similar to the third but with a higher number of lags (p).

Our proposed REVAR model combines the advantages and strengths of both reduced-rank and envelope VAR models, demonstrating its superiority over alternative methods in simulation studies.

To further investigate the performance of the proposed REVAR model, particularly in higher-dimensional cases, we demonstrate its effectiveness in four additional examples involving medium and large VARs. Table 2 compares the total number of parameters for each model with different combinations of $(d,u,p,q)$. Similar to Figure 1, Figure 2 illustrates the impact of envelope dimension and rank on the comparative performances of each approach. The top row of Figure 2 graphically presents the comparison of the four VAR models for medium-sized VARs with p = 1 and p > 1, while the bottom row shows the comparison for relatively large VARs. All simulation results demonstrate the superiority of our REVAR model with significant improvements over other VAR models. Due to computational costs, the 1D algorithm is used for the simulation studies of medium and large VARs, as it offers faster computation compared to the FG algorithm.

Fig. 2 Average estimation errors of higher-dimensional VAR models with different combinations of $(d,u,p,q)$ against sample size. Refer to Table 2 for details.

Table 2 Total number of parameters (NOP) in each VAR model correspond to Figure 2.

	Medium VAR	Medium VAR with $p\mathrm{}>\mathrm{}1$	Large VAR	Large VAR with $p\mathrm{}>\mathrm{}1$
$(d,u,p,q)$	$(5,10,1,20)$	$(5,10,2,20)$	$(10,20,1,40)$	$(10,20,2,40)$
OLSVAR	610	1010	2420	4020
RRVAR	385	485	1520	1920
EVAR	410	610	1620	2420
REVAR	335	435	1320	1720

Figure 3 summarizes the minimum and maximum asymptotic standard error ratios of the estimated coefficients for each of the OLSVAR, EVAR, and RVAR models in relation to the proposed REVAR model. The top row of Figure 3 displays these ratios versus the sample size when $(d,u,p,q)=(5,10,1,20)$, while the bottom row shows the ratios for $(d,u,p,q)=(10,20,1,40)$. All ratios in Figure 3 are greater than one, indicating that the REVAR model achieves the highest efficiency gains by significantly reducing the standard errors of the estimated coefficient matrix compared to the three comparative VAR models.

Fig. 3 Minimum (left panels) and maximum (right panels) asymptotic standard error ratios of coefficient estimates for each VAR model with respect to the proposed REVAR model.

6.2 Simulation Studies with Nonnormal Errors

In this section, simulations are conducted with nonnormal errors using ${\mathit{\epsilon }}_t={\mathbf{\Sigma }}^{1/2}{{\rm Y}}_t$, where ${{\rm Y}}_t$ is a vector of iid random variables with a mean of $\mathbf{0}$ and a covariance matrix of ${\mathbf{I}}_q$. We considered four different distributions, that is, normal, uniform, t₆, and ${\chi }_6^2$ distributions. Figure 4 compares the relative performance of each method across these distributions, in terms of the impact of the envelope dimension and rank. Our experiments demonstrated that the REVAR estimator consistently has the lowest average estimation error compared to the other models across all four distributions. Similar results were obtained for higher-dimensional VAR models, which are presented in Figures S1 and S2 in Supplement S8.1. All the minimum and maximum asymptotic standard error ratios of each VAR model relative to our proposed REVAR model exceed one, indicating that the REVAR model consistently achieves higher efficiency gains (see Figures S3, S4, and S5 in Supplement S8.1).

Fig. 4 Average estimation errors for $(d,u,p,q)=(3,4,1,7)$ with Normal, Uniform, t-Student, and Chi-square ${\mathit{\epsilon }}_t$ distributions, plotted against sample size. Refer to Table 1 for details.

Fig. 5 Average estimation errors of moderate-dimensional VAR models with different $(d,u,p,q)$ for stochastic volatility martingale difference sequence (SV-MDS) errors against sample size.

6.3 Simulation Studies with Martingale Difference Sequence Errors

In this section, we adopt a martingale difference error structure for the analysis. Let $\left\{{\mathit{\epsilon }}_t\right\}$ be a q-dimensional martingale difference sequence with respect to the increasing sequence of σ-fields $\{{\mathcal{F}}_t$}, where ${\mathit{\epsilon }}_t$ are ${\mathcal{F}}_t$-measurable and $E({\mathit{\epsilon }}_t|{\mathcal{F}}_t)=0$ for all $1\le t\le T$. We also introduce a martingale sequence $\{{\mathit{\zeta }}_t,\sigma ({\mathit{\zeta }}_t)\}$ for each t, and define the error term as ${\mathit{\epsilon }}_t={\mathit{\zeta }}_{t+1}-{\mathit{\zeta }}_t$. To simulate the error structure, we generate ${\mathit{\zeta }}_1$ from $\mathcal{N}(0,\mathbf{\Sigma })$ distribution and subsequently generate the rest of the sequence from the conditional distribution ${\mathit{\zeta }}_{t+1}|{\mathit{\zeta }}_t\sim \mathcal{N}({\mathit{\zeta }}_t,\mathbf{\Sigma })$. This simulation method has been used in previous studies (Zhao, Wang, and Zhang 2011; Zhou and Lin 2013). We have also implemented this simulation approach for nonnormal scenarios.

Furthermore, we adopt a stochastic volatility martingale difference sequence model to capture the time-varying volatility in the time series data. The error terms are modeled as $$\begin{array}{c}{\mathit{\epsilon }}_t={\mathbf{e}}_t\hspace{0.17em}\exp \hspace{0.17em}({\mathit{\sigma }}_t)\hspace{1em}\text{with}\hspace{1em}{\mathit{\sigma }}_t=0.25{\mathit{\sigma }}_{t-1}+0.05{\mathbf{u}}_t,\\ {\mathbf{e}}_t\stackrel{\text{iid}}{\sim }\mathcal{N}(\mathbf{0},\mathbf{\Sigma }),\hspace{1em}{\mathbf{u}}_t\stackrel{\text{iid}}{\sim }\mathcal{N}(\mathbf{0},{\mathbf{V}}_u),\end{array}$$ where ${\mathbf{V}}_u=({\nu }_{u,kl})\in {\mathbb{R}}^{p\times p}$ with ${\nu }_{u,kl}={0.9}^{|k-l|}$. This simulation method has been previously utilized by Chang, Jiang, and Shao (in press) and Escanciano and Velasco (2006) for the univariate and multivariate martingale difference hypothesis testing problems, respectively.

Figure S6 displays the impact of envelope dimension and rank when errors are generated from martingale difference sequences, while Figure 5 presents the results for errors generated from stochastic volatility martingale difference sequences. Across both figures, the simulation results consistently demonstrate the superior performance of our proposed REVAR model, revealing significant improvements over other VAR models. Similar findings were observed for higher-dimensional VAR models, as shown in Figures S7 and S8 in Supplement S8.2. Furthermore, all the minimum and maximum asymptotic standard error ratios of each VAR model relative to our proposed REVAR model were found to be greater than one, indicating the consistent efficiency gains achieved by the REVAR model (refer to Figures S9 and S10 in Supplement S8.2)

We conduct a pseudo-real-time forecasting exercise utilizing an expanding window of data to calculate h-step ahead forecasts (${\widehat{y}}_{t+h|t}$). The accuracy of these forecasts is evaluated using the root mean square forecast error (RMSFE): (20) $${\text{RMSFE}}_h\mathrm{}=\mathrm{}\sqrt{\frac{1}{T-T_0-\mathcal{H}+1}\sum _{t=T_0+\mathcal{H}-h}^{T-h}{({\widehat{y}}_{t+h|t}-y_{t+h})}^2}, h=1,\dots ,\mathcal{H},$$(20) where T₀ and T indicate the start and the end of the evaluation sample, respectively, and $\mathcal{H}$ is the maximum forecast horizon of interest. Our forecasting simulations encompass various error distributions, such as normal, nonnormal, and martingale difference errors with stochastic volatility (SV). The evaluation sample consists of the last 25% of the data.

Table S2 and Table S3 in Supplement S8.3 present a comprehensive summary of the root mean square forecast error (RMSFE_h ) for different forecast horizons ($h=1,2,3,4$), averaged over the evaluation samples across various error distributions and different DGP scenarios, with a sample size of T = 700. The results consistently demonstrate that our REVAR model outperforms other VAR models in terms of forecast accuracy. Furthermore, all average asymptotic standard error ratios relative to the REVAR model exceed one, indicating greater efficiency gains and improved accuracy in coefficient estimation compared to the comparative VAR models.

Table 3 summarizes the asymptotic Chi-squared test results for rank selection, envelope dimension, and lag order selections using the BIC criterion at a 0.05 significance level. We performed the simultaneous selection of d and u. The REVAR model with dimensions $(d,u,p,q)=(3,5,1,7)$ and $(2,4,2,8)$ are utilized, with a total NOP in REVAR model of 55 and 72, respectively. The table reports the percentages of correctly identified rank d and envelope VAR dimension u and lag order p for different sample sizes. Most of the other combinations of $(d,u,p,q)$ employed in the simulation studies yielded similar results, although in some cases, the selection of u and d was less accurate and resulted in higher percentages of overestimation. Consistent with the findings of Forzani and Su (2021), our results indicate that when the information criteria fail to select the true dimensions (u and d), they tend to overestimate these dimensions. Although this overestimation leads to a loss of efficiency, it does not introduce bias into the estimation process. Conversely, underestimation of d and u can indeed lead to biased outcomes. Nevertheless, overestimation of d and u is typically not a major concern (Cook, Forzani, and Zhang 2015).

Table 3 Percentage selection of the true dimensions (d and u) and lag order (p).

	$(d,u,p,q)=(3,5,1,7)$			$(d,u,p,q)=(2,4,2,8)$
	p Selection	d Selection	u Selection	p Selection	d Selection	u Selection
T	BIC	Chi-squared test	BIC	BIC	Chi-squared test	BIC
240	100%	93%	88%	100%	83%	54%
450	100%	95%	92%	100%	94%	90%
700	100%	97%	99%	100%	94%	99%
950	100%	98%	100%	100%	93%	98%
1200	100%	98%	100%	100%	94%	99%

7 Real Data Analysis

We analyzed four quarterly macroeconomic datasets obtained from the Federal Reserve Economic Quarterly Data (FRED-QD) website, with $T=244$ observations from 1959Q1 to 2019Q4. The first dataset, NIPA, consists of q = 8 macroeconomic variables from the National Income and Product Accounts, including real Gross Domestic Product (GDP) and its components. The second dataset consists of q = 11 variables related to Industrial Production, the third dataset includes q = 20 variables related to Price, and the fourth dataset encompasses q = 8 variables related to Money and Credit. Table S6 in Supplement S9 provides a detailed overview of the variables in each dataset. To ensure stationarity, appropriate transformations were applied to each variable group, such as taking the first difference of the logarithmic series for the 8 NIPA variables. Further details regarding variable descriptions and transformations can be found in McCracken and Ng (2020).

The rank ($\widehat{d}$) is determined using the Chi-squared test (Section 5.2) at a significance level of 0.01. The lag order ($\widehat{p}$) and the envelope dimension ($\widehat{u}$) are selected based on the BIC criterion as described in Sections 5.1, and 5.3, respectively. We conducted pseudo-real-time forecasting experiments using an expanding window scheme. The process begins with a specific initial fraction of the full sample (1959Q1-2004Q4) and recursively calculated forecasts over the evaluation sample (2005Q1-2019Q4). These forecasts are computed using the stationary bootstrap scheme proposed by Politis and Romano (1994) with 100 bootstrap samples. We assessed the forecasting performance of our proposed REVAR model using the ${\text{RMSFE}}_h$ for different forecast horizons (h) and compared it with other VAR models.

Tables 4 and 5 present the RMSFE_h values ($h=1,2,3,4$) for the NIPA and Price datasets, respectively. Supplementary Tables S4 and S5 in Supplement S9 provide the same results for the Money and Credit, and the Industrial Production datasets, respectively. The results consistently demonstrate the superior forecasting performance of our proposed REVAR model compared to the three competitive VAR models. Moreover, the average asymptotic standard error ratios ($r_{\text{avg}.}$) of the OLSVAR, EVAR, and RRVAR models relative to our proposed REVAR model are consistently greater than one, which indicates significant efficiency gains achieved by the REVAR model.

Table 4 Pseudo-real-time forecasting performance with bootstrap for the NIPA dataset (1959Q1-2019Q4) using an expanding window scheme, evaluated from 2005Q1 to 2019Q4.

				RMSFEh
$(\widehat{d},\widehat{u},\widehat{p},q)$	Model(M)	NOP	$r_{\text{avg}.}$	h = 1	h = 2	h = 3	h = 4
	OLSVAR	100	1.2912	0.059709	0.059433	0.059647	0.059671
(3, 4, 1, 8)	EVAR	68	1.0854	0.059133	0.05888	0.059074	0.05909
	RRVAR	75	1.0029	0.059278	0.058971	0.059047	0.059086
	REVAR	63	_	0.058950	0.058666	0.058850	0.058788

Table 5 Pseudo-real-time forecasting performance with bootstrap for the Price dataset (1959Q1-2019Q4) using an expanding window scheme, evaluated from 2005Q1 to 2019Q4.

				RMSFEh
$(\widehat{d},\widehat{u},\widehat{p},q)$	Model(M)	NOP	$r_{\text{avg}.}$	h = 1	h = 2	h = 3	h = 4
	OLSVAR	610	1.2236	0.13261	0.13114	0.13214	0.13187
(13, 14, 1, 20)	EVAR	490	1.0203	0.13194	0.13065	0.13156	0.13112
	RRVAR	561	1.1201	0.13228	0.13105	0.13182	0.13157
	REVAR	483	_	0.13184	0.13063	0.13150	0.13101

In order to evaluate the effectiveness of the proposed REVAR model with a larger number of variables, Table 6 summarizes the RMSFE_h values ($h=1,2,3,4$) using the first 43 variables from the FRED-QD dataset. These results not only show the superior performance of the REVAR model over other VAR models but also highlight its significant efficiency gains and forecasting improvements compared to smaller models.

Table 6 Pseudo-real-time forecasting performance with bootstrap for the initial 43 out of 47 macroeconomic variables in Table S6 (1959Q1-2019Q4), evaluated from 2005Q1 to 2019Q4 using an expanding window scheme.

				RMSFEh
$(\widehat{d},\widehat{u},\widehat{p},q)$	Model(M)	NOP	$r_{\text{avg}.}$	h = 1	h = 2	h = 3	h = 4
	OLSVAR	2795	1.4202	0.21493	0.21371	0.21381	0.21396
(18, 19, 1, 43)	EVAR	1763	1.0406	0.20709	0.20664	0.20735	0.20668
	RRVAR	2170	1.0594	0.21003	0.20956	0.20962	0.20925
	REVAR	1738	_	0.20644	0.20615	0.20681	0.20631

8 Discussion

In this article, we proposed the reduced-rank envelope VAR (REVAR) model, a parsimonious vector autoregressive model that incorporates the concept of envelopes into the reduced-rank VAR framework. The proposed REVAR approach offers a powerful approach to detect relevant and irrelevant information in VAR models and significantly improve estimation accuracy and efficiency by removing the immaterial information and variations from the VAR process. We derived maximum likelihood (ML) estimators for the REVAR model and investigated their asymptotic properties under various error assumptions, including normality, nonnormality, and martingale difference sequence error assumptions. Through extensive simulation studies and real data analysis, we compared the performance of the proposed REVAR model with that of the standard VAR model, EVAR model, and RRVAR model. Our results consistently demonstrated the superiority of the REVAR model in terms of estimation accuracy, efficiency gains, and forecasting performance across different scenarios and datasets. The technical details, including proofs of propositions and additional simulations, can be found in the Supplementary Materials, providing further insights into the properties of the REVAR model and its estimators.

The REVAR model exhibits remarkable speed. For instance, fitting the model with $(d,u,p,q)=(3,4,1,7)$ in Figure 1 averaged just 0.07123 sec, while fitting the $(5,10,1,20)$ model in Figure 2 took approximately 0.61587 sec, both with a sample size of $T=450$. These timings were achieved using an 8-core 3.2 GHz Apple M1 processor (Table S1).

For VAR models with a moderate number of variables and short lag lengths, OLS or ML estimations typically involve hundreds of parameters in coefficient and covariance matrices. More parameters can lead to less precise estimates and larger standard errors, especially with a fixed sample size. Applying the proposed REVAR model can parsimoniously reduce the number of model parameters, thereby improving the estimation efficiency and prediction accuracy. Evidence suggests that increasing the VAR dimension beyond 20 in economic applications does not yield substantial gains (see Banbura, Giannone, and Reichlin 2010; Koop 2013). However, our simulations and empirical results indicate that the REVAR model is capable of handling large datasets with up to 50 variables while maintaining robust performance.

Furthermore, our proposed REVAR model can be extended to a time-varying parameter model with stochastic volatility. In recent literature, there has been a growing recognition of the importance of stochastic volatility models in forecasting, as opposed to purely time-varying coefficients (Chan and Eisenstat 2018). It is worth noting that any subset of components within the decomposition of ${\Sigma }_t$ in (7) can vary over time. However, by employing a common reducing subspace for all ${\Sigma }_t$, it becomes feasible to implement a REVAR model with stochastic volatility. While these extensions are not straightforward, they are feasible and currently under investigation.

Acknowledgments

Thanks to the Editor, Professor Atsushi Inoue, an Associate Editor, and anonymous referees for their insightful and constructive feedback, which has greatly enhanced the article. Grateful to Professors Raja Velu and Dennis Cook for their helpful suggestions.

Supplementary Materials

The online supplementary materials contain all technical proofs of lemmas and propositions, along with a detailed dataset description, and additional simulation and empirical results.

Disclosure Statement

No potential conflict of interest was reported by the author(s).