Probabilistic Reduced-Dimensional Vector Autoregressive Modeling with Oblique Projections

Consider a time series $\{\bm{y}_{k}\in\Re^{p}\}_{k=1}^{N+s}$ with $E(\bm{y}_{k})=\bm{0}$ for all $k\in{\llbracket N+s\rrbracket}$. Qin [29] defines that the time series is serially correlated or dynamic if, for at least one $j>0$,

Otherwise, it is serially uncorrelated. Further, as per Qin [29], $\{\bm{y}_{k}\}$ is a reduced-dimensional dynamic (RDD) series if it is serially correlated but $\{\mathbf{a}^{\intercal}\boldsymbol{y}_{k}\}$ is serially uncorrelated for some $\mathbf{a\neq 0}\in\Re^{p}$. If no $\mathbf{a\neq 0}$ exists to render $\{\mathbf{a}^{\intercal}\bm{y}_{k}\}$ serially uncorrelated, $\{\bm{y}_{k}\}$ is full-dimensional dynamic (FDD). Qin [29] develops a latent VAR model to estimate the RDD dynamics of the data series. However, no statistical analysis is given since the latent VAR model does not give the distributions of noise terms.

In this paper, we define a probabilistic reduced-dimensional vector autoregressive series as

where $\{\bm{\varepsilon}_{k}\}$ and $\{\bm{\bar{\varepsilon}}_{k}\}$ are i.i.d. Gaussian noise terms and $\bm{v}_{k}\in\Re^{\ell}$ denotes the vector of $\ell$ DLVs. Both $\bm{\varepsilon}_{j}$ and $\bm{\bar{\varepsilon}}_{j}$ form the innovations and therefore are uncorrelated with the past $\bm{v}_{i}$ for all $i<j$. With $\ell<p$, the loadings matrices $\mathbf{P}\in\Re^{p\times\ell}$, $\mathbf{\bar{P}}\in\Re^{p\times(p-\ell)}$, and $\left[\mathbf{P}\leavevmode\nobreak\ \leavevmode\nobreak\ \mathbf{\bar{P}}\right]\in\Re^{p\times p}$ have full rank. Using the backward-shift operator $q^{-1}$, the latent dynamic system (2) can be converted to

Then, the transfer matrix form of (1)–(2) is given by

which clearly shows the reduced-dimensional dynamics that are commonly observed in routine operational data [29] and economics [34, 14]. Further, for weakly FDD time series, an RDD approximation is often useful to enable rapid downstream decision-making [13]. It should be noted that the latent VAR model (2) can be replaced by a state-space, an ARIMA, or a kernel-based model [8, 38, 27, 17].

The PredVAR model (1)–(2) is uniquely characterized by the tuple $\{\mathbf{P},\{\mathbf{B}_{j}\},\mathbf{\bar{P}}\}$, but it is implicit in the DLVs $\bm{v}_{k}$. To yield an explicit expression of $\bm{v}_{k}$, define $\mathbf{R}\in\Re^{p\times\ell}$ and $\mathbf{\bar{R}}\in\Re^{p\times(p-\ell)}$ such that

which makes $\left[\mathbf{R}\leavevmode\nobreak\ \leavevmode\nobreak\ \mathbf{\bar{R}}\right]$ unique given $\mathbf{P}$ and $\mathbf{\bar{P}}$. The relation (3) also implies

where $\mathbf{P}\mathbf{R}^{\intercal}$ and $\mathbf{\bar{P}}\mathbf{\bar{R}}^{\intercal}$ are two oblique projection matrices because $(\mathbf{P}\mathbf{R}^{\intercal})^{2}=\mathbf{P}\mathbf{R}^{\intercal}$ and $(\mathbf{\bar{P}}\mathbf{\bar{R}}^{\intercal})^{2}=\mathbf{\bar{P}}\mathbf{\bar{R}}^{\intercal}$.

The pair $(\mathbf{R},\mathbf{\bar{R}})$ gives explicitly the DLVs $\bm{v}_{k}$ and the static noise $\bm{\bar{\varepsilon}_{k}}$ from (1),

Equations (4) and (5) jointly explain the dynamic-static decomposition of $\bm{y}_{k}$ and reversely (1) gives the composition of $\bm{y}_{k}$. The reversible relations and the PredVAR model are depicted in Fig. 1 with a block diagram.

Given any pair of nonsingular matrices $\mathbf{\bar{M}}\in\Re^{(p-\ell)\times(p-\ell)}$ and $\mathbf{M}\in\Re^{\ell\times\ell}$, (1) can be rewritten as

which indicates that the tuple $(\mathbf{P},\bm{v}_{k},\mathbf{\bar{P}},\bm{\bar{\varepsilon}}_{k})$ cannot be uniquely identified since it can be exactly represented by $(\mathbf{P}\mathbf{M}^{-1},\mathbf{M}\bm{v}_{k},\mathbf{\bar{P}}\mathbf{\bar{M}}^{-1},\mathbf{\bar{M}}\bm{\bar{\varepsilon}}_{k})$. However, all equivalent realizations differ by similarity transformation only.

The non-uniqueness of the matrices offers the flexibility to extract desired coordinates or features of the DLVs. For example, we can enforce the following conditions to achieve identifiability.

1. 1.

   $\mathbf{P}$ orthogonal or the covariance $\mathbf{\Sigma}_{\bm{v}}=\mathbf{I}$ or diagonal, and independently,

2. 2.

   $\mathbf{\bar{P}}$ orthogonal or the covariance $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}=\mathbf{I}$ or diagonal if $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}$ is non-singular.

3. 3.

   If $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}$ is singular, we can represent $\bm{\bar{\varepsilon}}_{k}=\mathbf{Q}\bm{\bar{\varepsilon}}_{k}^{*}$ with $\mathbf{Q}$ being full-column rank and $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}^{*}}=\mathbf{I}$ having a lower dimension than $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}$.

In this paper, however, the DLV and static loadings matrices $\mathbf{P}$ and $\mathbf{\bar{P}}$ need not be mutually orthogonal, requiring an oblique projection for the DLV modeling.

These options have been used in the literature to enforce diagonal $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}$ and $\mathbf{\Sigma}_{\bm{\varepsilon}}$ as in [2], which does not consider rank-deficient innovations, or enforce $\mathbf{\Sigma}_{\bm{v}}=\mathbf{I}$ as in [29] to achieve the CCA objective. Further, the DLVs in [29] are arranged in a non-increasing order of predictability. A rank-deficient $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}}$ corresponds to a rank-deficient $\mathbf{\Sigma}_{\bm{y}}$, which is studied in this paper. These freedoms allow us to obtain desirable realizations of the PredVAR model.

A reduced-rank VAR (RRVAR) model is specified as follows: instead of estimating a full-rank coefficient matrix for each lag in a regular VAR model, we impose a factor structure on the coefficient matrices, such that they can be written as a product of lower-rank matrices, namely,

where $\mathbf{\acute{R}}^{\intercal}\mathbf{\acute{P}}\in\Re^{\ell\times\ell}$ is invertible with $\ell<p$ and $\bm{e}_{k}\in\Re^{p}$ is the full-dimensional innovations vector driving $\{\bm{y}_{k}\}$.

The reduced-rank VAR is analogous to the reduced-rank regression (RRR) for linear regression [16, 32]. Just as RRR requires specific solutions that differ from ordinary least squares, RRVAR also requires a specific solution differing from that of regular VAR. Theorem 1 shows that the PredVAR and RRVAR models are equivalent.

Theorem 1 is proven in Appendix A. Since PredVAR and RRVAR represent the same reduced-dimensional dynamics equivalently, one can choose either form for the convenience of the specific analysis and applications. It is worth noting that $\bm{e}_{k}$ is independent of the similarity transforms permitted in (6).

Given the measured data $\{\bm{y}_{k}\}_{k=1}^{N+s}$, the following likelihood function should be maximized:

Considering the RRVAR model (8), maximizing the likelihood amounts to minimizing the following function constrained by $\mathbf{R}^{\intercal}\mathbf{P}=\mathbf{I}$ and $\mathbf{R}^{\intercal}\mathbf{\Sigma}_{\bm{e}}\mathbf{\bar{R}}=\bm{0}$ from (9),

Then, it follows from (3), (4), and (10) that

For each $j\in{\llbracket s\rrbracket}$, differentiating $L^{\bm{y}}(\{\mathbf{B}_{j}\},\mathbf{P},\mathbf{\Sigma}_{\bm{e}},\mathbf{R})$ with respect to $\mathbf{B}_{j}$ and setting it to zero lead to

Rearranging the above equation gives

On the other hand, denoting the one-step prediction

it follows that (11) can be rearranged as

Differentiating (13) with respect to $\mathbf{P}$ and $\mathbf{\Sigma}_{\bm{e}}^{-1}$ and setting them to zero lead to

Substituting (14) into (15) leads to

It is clear that the necessary condition of optimality involves nonlinear equations of the parameters $(\{\mathbf{B}_{j}\},\mathbf{P},\mathbf{\Sigma}_{\bm{e}},\mathbf{R})$. To solve them iteratively, we separate the solution into two EM steps. First, given $\mathbf{\hat{R}}$ and $\{\mathbf{\hat{B}}_{j}\}$ respectively, the E-step of the EM procedure gives

for $k=1,2,\cdots,N+s$.

Second, (17) and (18) can be used to estimate the model parameters in the M-step. Formulate matrices

which include the formation of $\mathbf{V}_{s}$ and $\mathbf{Y}_{s}$. Substituting the estimated DLV vectors (17) into (12) leads to

Further, denoting

(18) and (19) lead to

where $\mathbf{\Pi}_{\mathbb{V}}\triangleq{\mathbb{V}}({\mathbb{V}}^{\intercal}{\mathbb{V}})^{-1}{\mathbb{V}}^{\intercal}$ is a projection matrix, which shows that $\mathbf{\hat{V}}_{s}$ is the orthogonal projection of $\mathbf{V}_{s}$ onto the range of $\mathbb{V}$. It follows from (14) and (16) that

where $\mathbf{\Pi}_{\mathbf{\hat{V}}_{s}}\triangleq\mathbf{\hat{V}}_{s}(\mathbf{\hat{V}}_{s}^{\intercal}\mathbf{\hat{V}}_{s})^{-1}\mathbf{\hat{V}}_{s}^{\intercal}$ is a projection matrix.

Finally, we need to update $\mathbf{\hat{R}}$ to complete the iteration loop. Based on (9), pre-multiplying $\mathbf{\hat{R}}^{\intercal}$ and post-multiplying $\mathbf{\hat{R}}$ to (22) lead to

From the orthogonal projection in (20), it is straightforward to show that $\mathbf{\hat{V}}_{s}^{\intercal}(\mathbf{V}_{s}-\mathbf{\hat{V}}_{s})=\mathbf{0}$, which gives $\mathbf{\hat{V}}_{s}^{\intercal}\mathbf{V}_{s}=\mathbf{\hat{V}}_{s}^{\intercal}\mathbf{\hat{V}}_{s}$. The last term in (23) becomes

Further, (23) becomes

The best predictive DLV model is one that has the smallest prediction error covariance $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$ given $\mathbf{\hat{R}}$ in the sense that $\mathbf{\hat{\Sigma}}^{\prime}_{\bm{\varepsilon}}-\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$ is positive semi-definite for any other $\mathbf{\hat{\Sigma}}^{\prime}_{\bm{\varepsilon}}$. However, since $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$ can be made arbitrarily small by scaling $\mathbf{\hat{R}}$, we need to fix the scale of one of the terms in (24). In [29], the scale is fixed by enforcing $\mathbf{V}_{s}^{\intercal}\mathbf{V}_{s}=\mathbf{I}$. In this paper, an equivalent scaling is implemented in the next subsection with

where $\mathbf{\hat{\Sigma}}_{\bm{y}}=\mathbf{Y}_{s}^{\intercal}\mathbf{Y}_{s}/N$. Then (24) becomes

Therefore, maximizing $\mathbf{\hat{\Sigma}}_{\hat{\bm{v}}}=\mathbf{\hat{R}}^{\intercal}\mathbf{Y}_{s}^{\intercal}\mathbf{\Pi}_{\mathbf{\hat{V}}_{s}}\mathbf{Y}_{s}\mathbf{\hat{R}}/N$ subject to (25) leads to the smallest $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$.

The constraint (25) can be implemented by normalizing the data based on $\mathbf{\hat{\Sigma}}_{\bm{y}}$, a technique used in [29, 14]. The constraint (25) indicates that $\mathbf{\hat{R}}$ can only be found in the eigenspace of $\mathbf{\hat{\Sigma}}_{\bm{y}}$ corresponding to the non-zero eigenvalues. Thus, assuming $\textrm{rank}(\mathbf{\hat{\Sigma}}_{\bm{y}})=r\leq p$, we must have the number of DLVs $\ell\leq r$.

By performing eigenvalue decomposition (EVD) on

where $\mathbf{U}\in\Re^{p\times r}$, $\mathbf{\tilde{U}}\in\Re^{p\times(p-r)}$, and $\mathbf{D}\in\Re^{r\times r}$ contains the non-zero eigenvalues in non-increasing order, we define a normalized data vector $\bm{y}_{k}^{*}\in\Re^{r}$ with

which represents $\bm{y}_{k}^{*}$ with the following normalized PredVAR model

where $\mathbf{P}^{*}\in\Re^{r\times\ell}$, $\mathbf{\bar{P}}^{*}\in\Re^{r\times(r-\ell)}$, and we make $\mathbf{\Sigma}_{\bm{\bar{\varepsilon}}^{*}}=\mathbf{I}$ without any loss of generality. Then, it follows from (26) and (27) that

Forming $\mathbf{Y}_{i}^{*}$ the same way as $\mathbf{Y}_{i}$ for $i=0,1,\cdots,s$, it is straightforward to show that $\mathbf{\hat{\Sigma}}_{\bm{y}^{*}}=\mathbf{Y}_{s}^{*\intercal}\mathbf{Y}_{s}^{*}/N=\mathbf{I}$. By using the uncorrelated condition of the DLVs and static noise, we have

i.e., $\mathbf{P}^{*\intercal}\mathbf{P}^{*}=\mathbf{I}$, $\mathbf{\bar{P}}^{*\intercal}\mathbf{\bar{P}}^{*}=\mathbf{I}$, and $\mathbf{P}^{*\intercal}\mathbf{\bar{P}}^{*}=\mathbf{0}$. Comparing the above relation to (3), it is clear that $\mathbf{P}^{*}$ and $\mathbf{R}^{*}$ coincide with each other, i.e., $\mathbf{R}^{*}=\mathbf{P}^{*}$. Therefore, we generate the DLVs with

giving $\mathbf{R}=\mathbf{U}\mathbf{D}^{-\frac{1}{2}}\mathbf{R}^{*}$. Therefore, the optimal PredVAR solution is converted to finding $\mathbf{R}^{*}$, and thus $\mathbf{R}=\mathbf{U}\mathbf{D}^{-\frac{1}{2}}\mathbf{R}^{*}$, to make $\bm{v}_{k}$ most predictable.

To find the estimate $\mathbf{\hat{R}}^{*}$, the constraint (25) becomes

which again makes $\mathbf{\hat{R}}^{*}$ an orthogonal matrix. Equation (23) becomes

Therefore, to minimize $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$ subject to (30) with $\ell$ DLVs, $\mathbf{\hat{R}}^{*}$ must contain the eigenvectors of the $\ell$ largest eigenvalues of $\mathbf{Y}_{s}^{*\intercal}\mathbf{\Pi}_{\mathbf{\hat{V}}_{s}}\mathbf{Y}_{s}^{*}/N$. Performing EVD on

where $\mathbf{\Lambda}$ contains the eigenvalues in non-increasing order, it is convenient to choose

as the optimal solution, making

contain the leading canonical correlation coefficients (i.e., $R^{2}$ values) of the DLVs and

the residual covariance. It is noted that if some diagonal elements of $\mathbf{\Lambda}(1:\ell,1:\ell)$ are zero, the noise covariance $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$ is rank-deficient, and the corresponding DLVs have perfect predictions.

To summarize, the whole EM iteration solution boils down to updating $\mathbf{\hat{R}}^{*}$ and $\mathbb{\hat{B}}$ until convergence. With the converged $\mathbf{\hat{R}}^{*}$ and $\mathbb{\hat{B}}$, we can calculate $\mathbf{\hat{V}}_{s}$, $\mathbf{\hat{R}}$, $\mathbf{\hat{P}}$ from (28), and $\mathbf{\hat{\Sigma}}_{\bm{e}}$ from (22).

The static part of the model can be calculated from the oblique complement of the dynamic counterpart. It is straightforward to calculate $\mathbf{\hat{\bar{P}}^{*}}$ from (29). From (26) and (27), we have

Therefore, to match (1) and satisfy (3), we can choose

which give the static noise

The pseudocode of the PredVAR algorithm is shown in Algorithm 1. The data are first normalized based on $\mathbf{\hat{\Sigma}}_{\bm{y}}$ before the EM iteration. Then, the estimates $(\mathbf{\hat{R}^{*}},\mathbb{\hat{B}})$ are iterated with the EM-based PredVAR solution. After convergence, we obtain $(\{\mathbf{\hat{B}}_{j}\},\mathbf{\hat{P}},\mathbf{\hat{R}},\mathbf{\hat{\Sigma}}_{\bm{e}})$ for the RRVAR model and $(\mathbf{\hat{\bar{P}}},\mathbf{\hat{\bar{R}}},\mathbf{\hat{\Sigma}}_{\bm{\bar{\varepsilon}}})$ for the static noise model.

In real applications, it is possible that the covariances $\mathbf{\hat{\Sigma}}_{\bm{y}}$, $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$, and $\mathbf{\hat{\Sigma}}_{\bm{\bar{\varepsilon}}}$ can be rank-deficient. Algorithm 1 handles the rank deficiency of $\mathbf{\hat{\Sigma}}_{\bm{y}}$ by the initial EVD step. Since $\mathbf{\hat{\Sigma}}_{\bm{v}}=\mathbf{I}$ is enforced in this algorithm, the null space of $\mathbf{\hat{\Sigma}}_{\bm{y}}$ cannot be part of $\bm{v}_{k}$, and therefore, it must be part of $\mathbf{\hat{\Sigma}}_{\bm{\bar{\varepsilon}}}$. On the other hand, with the enforced $\mathbf{\hat{\Sigma}}_{\bm{v}}=\mathbf{I}$, it is possible to have rank-deficient $\mathbf{\hat{\Sigma}}_{\bm{\varepsilon}}$, which corresponds to DLVs perfectly predictable. Therefore, we give the following remarks.

We give the following effective strategy for initializing $\hat{\mathbf{{R}}}^{*}$, which is consistent with the updating of $\hat{\mathbf{{R}}}^{*}$ in Algorithm 1. Specifically, similar to $\mathbb{V}$, form the augmented matrix

Replacing $\mathbf{\hat{V}}_{s}$ with $\mathbb{Y}^{*}$ in (32) and performing EVD give rise to

where $\mathbf{\Pi}_{\mathbb{Y}^{*}}=\mathbb{Y}^{*}(\mathbb{Y}^{*\intercal}\mathbb{Y}^{*})^{-1}\mathbb{Y}^{*\intercal}$ and $\mathbf{\Lambda}_{0}$ contains the eigenvalues in a non-increasing order, we can initialize $\mathbf{\hat{R}}^{*}$ as $\mathbf{W}_{0}(:,1:\ell)$.

The underlying idea is to estimate the error covariance matrix $\mathbf{\hat{\Sigma}}_{\bm{e}^{f}}$ from the following full-rank VAR model

which makes

Interestingly, the canonical analysis of time series in [3] estimates the VAR error covariance matrix first, and then performs EVD to select $\ell$ canonical components for the smallest eigenvalues of the error covariance matrix, which is equivalent to the selection of $\mathbf{\hat{R}}^{*}$ in this initialization step. Therefore, the work[3] serves as an initial step of the PredVAR iteration only.

The PredVAR solution finds $\mathbf{R}$ to minimize the prediction error covariance of DLVs subject to $\mathbf{R}^{\intercal}\mathbf{P}=\mathbf{I}$. Thus, given $\mathbf{P}$ and $\mathbf{\Sigma}_{\bm{e}}$, we solve the following optimization problem under the Loewner order [4]

for the best predictability. The optimization result leads to the following theorem that gives rise to uncorrelated innovations of $\bm{\varepsilon}_{k}$ and $\bm{\bar{\varepsilon}}_{k}$.