Stable Reduced-Rank VAR Identification

Abstract

The vector autoregression (VAR) has been widely used in system identification, econometrics, natural science, and many other areas. However when the state dimension becomes large the parameter dimension explodes. So rank reduced modelling is attractive and is well developed. But a fundamental requirement in almost all applications is stability of the fitted model. And this has not been addressed in the rank reduced case. Here, we develop, for the first time, a closed-form formula for an estimator of a rank reduced transition matrix which is guaranteed to be stable. We show that our estimator is consistent and asymptotically statistically efficient and illustrate it in comparative simulations.

1 Introduction

The vector autoregression (VAR) is widely applied in control, signal processing, time series and econometrics. The lag-one VAR (VAR(1)) is more general than it seems since any higher order VAR can be written as a state space model whose state equation is a VAR(1) [12],[27]. More recently, the VAR(1) has been used in reinforcement learning [23, 1, 17, 19].

High dimensional time series have now become very common across many disciplines, e.g. epidemiology [28], financial economics [20]. So rank reduced modelling has gained increasing attention [20],[4],[5].

In physical system, stability is crucial. And this has led to a continuing development of stability guaranteed estimators, e.g. for state space modelling. A common approach employs a perturbation minimization procedure, where a preliminary unstable least-squares estimator is stabilized by the smallest additive perturbation. For example, Mari et al. [15] consider a semidefinite programming (SP) problem with a Lyapunov stability constraint. The problem can be solved by standard linear matrix inequality (LMI) methods. However, the computation becomes prohibitive with higher-dimensional states. Miller and de Callafon [16] extend the SP method to include more eigenvalue constraints, besides stability. Boots et al. [2] and the first algorithm in [14] use singular value constraints instead of eigenvalue constraints, which can be too conservative. The second algorithm in [14] uses a line search method based on the gradient sampling and the reliability is sensitive to the user-defined step size parameter. Tanaka and Katayama [24] and Jongeneel et al. [10] project the unstable estimator onto a stable region by solving linear quadratic regulator (LQR) problems. Jongeneel et al. [10] show that the computational efficiency is greatly improved in higher-dimension applications, and provide error analysis and statistical guarantees.

Other approaches include, Chui and Maciejowski [6] who iteratively augment the unstable estimator until its largest eigenvalues have the modulus of the user-defined parameter. However, this method distorts the estimator and introduces additional bias. Van Gestel et al. [26] introduce a regularization term to the least squares problem so that the variance of the estimators is reduced and in the meantime, the stability is guaranteed. However, the regularization is too conservative. Umenberger et al. [25] consider the state space modelling and take the maximum likelihood approach to guarantee stability. However, their method is computationally expensive.

In recent work, we [21] solved a state space estimation problem by using a Burg-type forwards-backwards (FB) optimization [3, 22, 18] on the error residuals and remarkably find a closed-form stable solution that is computationally cheap and involves no tuning parameters.

Among the above works, [10] deal with pure VAR(1) and the others consider the state space setup where subspace identification methods are needed before fitting a VAR(1) model [27]. However, none of the above stability-enforced methods deals with the reduced rank modeling.

In this paper, we extend the work in [21], using the FB approach to obtain a reduced rank (RR)-VAR(1) with guaranteed stability. The FB approach generates an estimator that is no more computationally expensive than the least squares RR method of [20]. Further, no tuning parameters are required.

The rest of the paper is organized as follows. We first introduce the forwards and backwards VAR(1) models in Section 2, and then develop for the first time:

1. in Section 3 a closed-form stable estimator for the full-rank VAR(1),

2. In Section 4, based on (i), a closed-form stable estimator for the RR-VAR(1),

3. in Section 5, statistical consistency, and asymptotic efficiency for the new estimator.

The results are illustrated in comparative simulations in Section 6. Section 7 contains conclusions.

We use the following notations.
$\rho(A)$ is the spectral radius of $A$.
$\|A\|=\sqrt{\operatorname{tr}(AA^{\prime})}$ is the Frobenius norm of $A$.
$A>0$ means $A$ is positive definite.
$A\geq 0$ denotes positive semi-definiteness.
w.p.1 means with probability 1.

2 VAR(1) and the Backwards Model

In this section, we review the (forwards) VAR(1) and its associated backwards model. The VAR(1) is a Markov process generated by

$$y_{t}=Fy_{t-1}+w_{t},\quad t=1,\dotsm,T,$$ (2.1)

where $y_{t}$ is the observed $n-$vector time series, $w_{t}$ is a zero-mean driving white noise with a non-singular covariance matrix $Q$, and $F$ is the transition matrix, assumed to be stable. The stability of $F$ ensures the existence of a steady state variance matrix $\Pi$ which obeys a discrete-time (DT) Lyapunov equation

$$\Pi=F\Pi F^{\prime}+Q.$$

A sufficient condition for $\Pi$ to be positive definite is that $Q$ is positive definite. We introduce the steady state lag-one cross-covariance $\Pi_{10}=\operatorname{E}[y_{t}y_{t-1}^{\prime}]=F\Pi=\Pi_{01}^{\prime}$, so that $F=\Pi_{10}\Pi^{-1}$. We further introduce the associated correlation matrix $R_{F}=\Pi^{-\frac{1}{2}}\Pi_{10}\Pi^{-\frac{1}{2}}=\Pi^{-\frac{1}{2}}F\Pi^{\frac{1}{2}}$ which has the same eigenvalues as $F$.

Associated with the VAR(1) is a stationary backwards model [11]

$$y_{t-1}=F_{b}y_{t}+w_{b,t-1},\quad t=T,\dotsm,1,$$

where $F_{b}$ is the backwards transition matrix and $w_{b,t-1}$ is a white noise with covariance matrix $Q_{b}$ which is statistically independent of $y_{t}$. We note that [11]

$$F_{b}=\operatorname{E}[y_{t-1}y_{t}^{\prime}]\operatorname{E}^{-1}[y_{t}y_{t}]=\Pi_{10}^{\prime}\Pi^{-1}=\Pi F^{\prime}\Pi^{-1},$$

and

$$Q_{b} =\operatorname{E}[(y_{t-1}-F_{b}y_{t})(y_{t-1}-F_{b}y_{t})^{\prime}]$$

$$=\Pi-F_{b}\Pi_{10}-\Pi_{10}^{\prime}F_{b}^{\prime}+F_{b}\Pi F_{b}^{\prime}$$

$$=\Pi-\Pi_{01}\Pi^{-1}\Pi_{10}$$

$$=\Pi-F_{b}\Pi F_{b}^{\prime}.$$

An elementary argument shows that $F_{b}$ has the same eigenvalues of $F$. Thus $F_{b}$ is stable iff $F$ is stable. The following results are used below.

Theorem 1. Converse Lyapunov result. If $\Pi$ is positive definite and $Q$ is positive semi-definite then any eigenvalue $\lambda$ of $F$ has $|\lambda|\leq 1$. If $Q$ is positive definite then $|\lambda|<1$.

Proof. Let $v$ be left eigenvector of $F$ with eigenvalue $\lambda$. Then denoting $(v^{\ast})^{\prime}$ as $v^{H}$ we find

$$\begin{array}[]{rrcl}&v^{H}\Pi v&=&|\lambda|^{2}v^{H}\Pi v+v^{H}Qv\\ \Rightarrow&(1-|\lambda|^{2})v^{H}\Pi v&=&v^{H}Qv\geq 0\Rightarrow|\lambda|\leq 1.\end{array}$$

Clearly if $Q$ is positive definite we get $|\lambda|<1$. $\square$

Theorem 2. $Q_{b}$ has full rank iff $Q$ has full rank.

Proof. First note that $Q_{b}=\Pi Q_{c}\Pi$ where $Q_{c}=\Pi^{-1}-F^{\prime}\Pi^{-1}F$. Then clearly $Q_{b}$ has full rank iff $Q_{c}$ has full rank. Suppose for some vector $v\neq 0$, we have $Q_{b}v=0$. Then $Q_{c}\Pi v=0$. Thus, $\Pi^{-1}\Pi v=F^{\prime}\Pi^{-1}F\Pi v$, i.e. $v=F^{\prime}\Pi^{-1}F\Pi v$. Note that we cannot have $F\Pi v=0$ since then $v=0$ which is a contradiction.

Now consider that

$$\begin{array}[]{rrcl}&Q\Pi^{-1}&=&I-F\Pi F^{\prime}\Pi^{-1}\\ \Rightarrow&Q\Pi^{-1}F\Pi&=&F\Pi-F\Pi F^{\prime}\Pi^{-1}F\Pi\\ \Rightarrow&Q\Pi^{-1}F\Pi v&=&F\Pi v-F\Pi F^{\prime}\Pi^{-1}F\Pi v\\ &&=&F\Pi v-F\Pi v=0\end{array}$$

which is a contradiction and proves ‘if’; then ‘only if’ follows by running the argument in reverse. $\square$

3 Stable Estimators for full-rank VAR(1)

Here we review earlier work of [18, 22, 21] and develop a new result on the full-rank case. Given data $y_{t},t=0,\cdots,T$ set $Y_{0}=[\begin{matrix}y_{0}&\dotsm&y_{T-1}\end{matrix}]$ and $Y_{1}=[\begin{matrix}y_{1}&\dotsm&y_{T}\end{matrix}],$ and define the sample covariances $S_{ij}=\frac{1}{T}Y_{i}Y_{j}^{\prime},\quad i,j\in\{0,1\}$. We also introduce the forwards and backwards residual mean squared errors

$$S_{w,f}(F) =\frac{1}{T}(Y_{1}-FY_{0})(Y_{1}-FY_{0})^{\prime}$$

$$=S_{11}-FS_{01}-S_{10}F^{\prime}+FS_{00}F^{\prime}$$

$$S_{w,b}(F) =\frac{1}{T}(Y_{0}-F_{b}Y_{1})(Y_{0}-F_{b}Y_{1})^{\prime}$$

$$=S_{00}-F_{b}S_{10}-S_{01}F_{b}^{\prime}+F_{b}S_{11}F_{b}^{\prime},$$

where $F_{b}=PF^{\prime}P^{-1}$, where $P$ is a consistent estimator of $\Pi$ to be chosen. Here $F$ is no longer the true value.

The least squares estimator $\hat{F}_{LS}=S_{10}S_{00}^{-1}$ minimises $\operatorname{tr}(S_{w,f})$, but it is NOT guaranteed to be stable. Continuing, [18, 22] considered the minimizer of the sum of weighted forwards and backwards sample mean squared errors

$$J(F;P)=\operatorname{tr}\{P^{-1}(S_{w,f}+S_{w,b})\}.$$

This yields the following remarkable result.

Theorem 3A. [22, 21, 18].
The minimizer of $J(F;P)$ is $\hat{F}$, the solution obeys the Sylvester equation below and is stable

$$\hat{F}S_{00}P^{-1}+S_{11}P^{-1}\hat{F}=2S_{10}P^{-1}.$$

Remarks.

1. This result is not explicitly stated in [22] but is rather a special case of his results. Further the argument in [22] is very hard to follow. Partly for that reason, [21] gave a simple direct proof. Note that Theorem 4B below has nearly the same proof.

2. The result only holds with the weighting matrix $P$ used to define $F_{b}$, not for any other weighting matrix.

3. Using vec algebra we can get a closed-form solution.

We now have the further remarkable closed-form result.

Theorem 3B. Set $P=S_{11}$. Then the stable estimator

$$\hat{F}=\hat{F}_{11}=2S_{10}(S_{00}+S_{11})^{-1},$$

Proof. Follows from Theorem 3A by setting $P=S_{11}$.

Remarks.

1. Below we find that the choice $P=S_{11}$ enables our new results.

2. Note that consistency and asymptotic efficiency has been established for $\hat{F}_{LS}$ in [12]. One can then show the same for $\hat{F}$. We omit details, concentrating here on the reduced rank case.

4 Stable Reduced Rank VAR(1) Estimation

We now consider the case where $F$ has rank $m<n$. Then we can factor $F=A_{n\times m}B_{m\times n}$. We then have the following result.

Theorem 4A. [20][Section 2.3]. Reduced Rank Least Squares.

$$\min_{\operatorname{rank}(F)=m}J_{LS}:J_{LS}=\operatorname{tr}(S_{11}^{-1}S_{w,f})$$

has solution $\hat{F}_{RLS}=S_{11}^{\frac{1}{2}}\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}S_{11}^{-\frac{1}{2}}\hat{F}_{LS}$ where $\hat{V}_{*,m}=[\begin{matrix}\hat{v}_{*,1}&\dotsm&\hat{v}_{*,m}\end{matrix}]$ and $\hat{v}_{\ast,r}$ are the ‘top’ $m-$eigenvectors of $\hat{R}_{*}=S_{11}^{-\frac{1}{2}}S_{10}S_{00}^{-1}S_{01}S_{11}^{-\frac{1}{2}}.$ The noise covariance estimator is $\hat{Q}_{RLS}=S_{11}^{\frac{1}{2}}(I-\hat{V}_{*,m}\hat{D}_{*,m}^{2}\hat{V}_{*,m}^{\prime})S_{11}^{\frac{1}{2}}$ where $\hat{D}_{*,m}^{2}=\operatorname{diag}(\hat{\lambda}_{*,k})$ contains the top $m$ eigenvalues of $\hat{R}_{*}$.

Remarks.

1. $\hat{F}_{RLS}$ is not guaranteed to be stable.

2. We note for future reference that $\hat{R}_{*}$ has the same eigenvalues as $\hat{L}_{*}=S_{10}S_{00}^{-1}S_{01}S_{11}^{-1}$.

3. We can re-express the results in terms of the SVD of $\hat{G}_{*}=S_{11}^{-\frac{1}{2}}S_{10}S_{00}^{-\frac{1}{2}}=\hat{V}_{*}\hat{D}_{*}\hat{U}_{*}^{\prime}$ since $\hat{R}_{*}=\hat{G}_{*}\hat{G}_{*}^{\prime}$. It follows from the Biename-Cauchy-Schwarz inequality that the singular values of $\hat{G}_{*}$ are $<1$ w.p.1.

4. [20] gives a different formula for $\hat{Q}_{RLS}$ but it is straightforward to show it equals the one given here.

We would now like to set up a FB version of this problem. But it turns out that for general $P$ there is no simple solution. Fortunately if we choose $P=S_{11}$ we can find a simple guaranteed stable estimator.

Theorem 4B. The solution to $\min_{\operatorname{rank}(F)=m}J(F;S_{11})$ is given by $\hat{F}_{R}=S_{11}^{\frac{1}{2}}\hat{V}_{m}\hat{V}_{m}^{\prime}S_{11}^{-\frac{1}{2}}\hat{F}_{11}$ where $\hat{V}_{m}=[\begin{matrix}\hat{v}_{1}&\dotsm&\hat{v}_{m}\end{matrix}]$ and $\hat{v}_{r}$ are the ‘top’ $m$ eigenvectors (corresponding to the top $m$ eigenvalues in $\hat{D}_{m}=\operatorname{diag}(\hat{\lambda}_{k})$) of

$$\hat{R}=2S_{11}^{-\frac{1}{2}}S_{10}(S_{00}+S_{11})^{-1}S_{01}S_{11}^{-\frac{1}{2}}.$$

Further $\hat{F}_{R}$ is stable.

Proof. See the appendix.

5 Asymptotic Analysis

Here we show the consistency and central limit theorem (CLT) for $\hat{F}_{R}$ by reducing them to results for $\hat{F}_{RLS}$.

For asymptotic results, we need stronger assumptions. There are a wide range of possibilities, but to keep things simple, while still retaining reasonable generality we use:

Assumption A1. $w_{t}$ in (2.1) are iid with finite fourth moments, and $Q$ is positive definite.

Assumption A2. $F$ has rank $m$ and the non zero eigenvalues of $F$ are distinct.

This enables the following result for $\hat{F}_{RLS}$ (where $F$ is again the true value).

Theorem 5A. Assuming $m=\operatorname{rank}(F)$ is known

1. Under A1: $S_{00}\xrightarrow{p}\Pi$ and $S_{11}\xrightarrow{p}\Pi$.
   

2. Under A1: $\hat{F}_{LS}\xrightarrow{p}F\equiv S_{10}\xrightarrow{p}F\Pi$.

Under A1, A2: $R_{F}=\Pi^{-\frac{1}{2}}F\Pi^{\frac{1}{2}}$ has rank $m$ and:

1. $\hat{V}_{*,m}\xrightarrow{p}V_{m}$ the eigenvectors of $R_{F}R_{F}^{\prime}$.
   $\hat{U}_{*,m}\xrightarrow{p}U_{m}$ the eigenvectors of $R_{F}^{\prime}R_{F}$.
   $\hat{\lambda}_{k}\xrightarrow{p}\lambda_{k}\mbox{ and }\hat{\lambda}_{*,k}\xrightarrow{p}\lambda_{k}$ the $k^{th}$ largest eigenvalue of $R_{F}R_{F}^{\prime}$.

2. $\hat{G}_{*}\xrightarrow{p}R_{F}=\Pi^{-\frac{1}{2}}F\Pi^{\frac{1}{2}}$.

3. $\hat{F}_{RLS}\xrightarrow{p}F=AB$.

4. $\sqrt{T}\operatorname{vec}(\hat{F}_{RLS}^{\prime}-F^{\prime})\Rightarrow Z\sim N(0,\Sigma)$ where
   $\Sigma=MWM^{\prime}$ with $M=[(I\otimes B^{\prime}),(A\otimes I)]$ and $W$ (which is complicated) is given in [20](2.36).

Proof. (i),(ii) can be found in [8][chapter 11]. (v),(vi) are from [20][Section 2.5 (Theorem 2.4) and Section 5.2]. For (iii) first note that from (i),(ii) $\hat{L}_{*}\xrightarrow{p}L=F\Pi F^{\prime}\Pi^{-1}$ which has the same eigenvalues as $R=R_{F}R_{F}^{\prime}$. The second part follows in the same way. Next, since distinct eigenvalues and eigenvectors are continuous functions of the underlying matrix entries [13] then $\hat{V}_{*,m}\xrightarrow{p}V_{m}$ containing the ‘top’ eigenvalues of $R$. Further, the eigenvalues $\xrightarrow{p}$ as stated. (iv) now follows from (iii).

To continue for $\hat{F}_{R}$, we need some lemmas.

Lemma 5B. Let $\Gamma_{t}=\operatorname{E}[y_{t}y_{t}^{\prime}]$ and denote $\rho=$ spectral radius of $F$. Then $\|\Gamma_{t}-\Pi\|_{2}\leq c\rho^{2t}$ for some constant $c$.

Proof. Taking variances in (2.1) gives $\Gamma_{t}=F\Gamma_{t-1}F^{\prime}+Q$. Subtracting the $\Pi$ equation from this gives $\Gamma_{t}-\Pi=F(\Gamma_{t-1}-\Pi)F^{\prime}$. Iterating this gives $\Gamma_{t}-\Pi=F^{t}D(F^{\prime})^{t}$ where $D=\Gamma_{0}-\Pi$. Now recall Gelfand’s theorem [9]: $\|F^{t}\|^{1/t}\rightarrow\rho$, so that $limsup_{t}\|(\Gamma_{t}-\Pi)/\rho^{2t}\|^{1/t}\leq c$. $\square$

The result follows from this.

Lemma 5C. Under A1: $\sqrt{T}\|S_{11}-S_{00}\|\xrightarrow{p}0$.

Proof. We have

$$\begin{array}[]{rrcl}&S_{11}-S_{00}&=&\frac{1}{T}[y_{T}y_{T}^{\prime}-y_{0}y_{0}^{\prime}]\\ \Rightarrow&\|S_{11}-S_{00}\|^{2}&\leq&\frac{2}{T}(\|y_{T}\|^{2}+\|y_{0}\|^{2})\\ \Rightarrow&\sqrt{T}\operatorname{E}\|S_{11}-S_{00}\|^{2}&\leq&\frac{2}{\sqrt{T}}[\operatorname{tr}(\Gamma_{T})+\operatorname{tr}(\Gamma_{1})]\end{array}$$

The second term $\rightarrow 0$ while the first term $\rightarrow 0$ from Lemma 5B. The result then follows. $\square$

Lemma 5D. Under A1, A2: $\sqrt{T}\|\hat{F}_{LS}-\hat{F}_{11}\|\xrightarrow{p}0$.

Proof. We have

$$\begin{array}[]{rcl}&&\sqrt{T}(\hat{F}_{LS}-\hat{F}_{11})\\ &=&\sqrt{T}(S_{10}S_{00}^{-1}-2S_{10}(S_{00}+S_{11})^{-1})\\ &=&\sqrt{T}2S_{10}[(2S_{00})^{-1}-(S_{11}+S_{00})^{-1}]\\ &=&\sqrt{T}S_{10}S_{00}^{-1}[S_{00}+S_{11}-2S_{00}](S_{00}+S_{11})^{-1}\\ &=&\hat{F}_{LS}[\sqrt{T}(S_{11}-S_{00})](S_{11}+S_{00})^{-1}\end{array}$$

Now from Theorem 5A and Lemma 5C each term converges in probability, so the product also converges in probability. But $\sqrt{T}(S_{11}-S_{00})\xrightarrow{p}0$ yielding the result.

Theorem 5E. Under A1, A2,

1. $\hat{F}_{R}\xrightarrow{p}F$.

2. $\sqrt{T}(\hat{F}_{RLS}-\hat{F}_{R})\xrightarrow{p}0$.

3. $\sqrt{T}\operatorname{vec}(F_{R}^{\prime}-F^{\prime})\Rightarrow Z\sim N(0,\Sigma)$.

Proof. (i) follows from (ii) and Theorem 5A(v). (iii) follows from (ii) and Theorem 5A(vi). (ii) is proved in the appendix. $\square$

Figure 1: Pole locations of the first $50$ repeats.

Figure 2: Pole magnitude histograms of the $1000$ repeats: $*$ are the true poles. Unstable poles for LS for each $T$: $25.9\%,6.9\%,0.2\%$.

6 Simulations

We now show simulations to illustrate our stable estimator $\hat{F}_{R}$. Since there is no existing algorithm that guarantees a stable RR-VAR(1) estimator, we compare our estimator with the standard, reduced-rank, least square estimator $\hat{F}_{RLS}$. We call our method the FB method and the latter the LS method. We expect the LS method will generate unstable estimates while FB does not. Also that LS and FB will have similar computational complexity given their formulae.

We run two sets of simulations:

• Low dimensional state $n=6$.

• High dimensional state, up to $n>3000$.

The latter study represents increasingly common practical examples.

6.1 Simulation Design

Firstly, we specify $F=\left[\begin{smallmatrix}F_{0}&0\\ 0&0_{3\times 3}\end{smallmatrix}\right]\otimes I_{p}$, where $p$ is to be chosen and $F_{0}=\left[\begin{smallmatrix}0.99&-0.1&0\\ 0.1&0.99&0\\ 0&0&0.95\end{smallmatrix}\right]$. So, there are $3p$ zero eigenvalues and $3p$ non-zero eigenvalues repeating at $0.99\pm\jmath 0.1$ and $0.95$, very close to the unit circle. This is a typical case when a stability-guaranteed algorithm is needed. This is a similar setup as in [10].

Secondly, we let the length $T$ of the time series be a multiple of model order $n$, so that the computational complexity for both LS and FB is $O(n^{3})$ (assuming a standard matrix multiplication method). Also, note that in the one-dimensional case, the precision $\frac{|\hat{F}_{LS}-F|}{|F|}$ decays at a rate of $\sqrt{T}$. Thus, we consider $T/n=l^{2}\Rightarrow T=l^{2}n$. The values of $l$ are best determined by system time constants. However, here we take a trial-and-error method and let $T/n=2^{2},6^{2},10^{2}\Rightarrow T=4n,36n,100n$.

Thirdly, we choose $p=2^{k}\Rightarrow n=6\times 2^{k},k\in\mathbb{N}$, so that the logarithm of the computational time is proportional to $k$.

Throughout this section, we assume the true model order $m$ is known. Also, throughout this section, we take the noise covariance $Q=I_{n}$.

6.2 Low Dimensional Simulations

We let $k=0\Rightarrow p=1\Rightarrow n=6$, and the rank $m=3$. We consider $3$ time series lengths $T/n=2^{2},6^{2},10^{2}\Rightarrow T=24,216,600$. We simulate $N=1000$ realizations for each $T$ and get the estimates $\hat{F}_{R}$ and $\hat{F}_{RLS}$.

We first plot the estimated pole locations in Fig. 1. We only plot the first $50$ estimates due to crowdedness. It is sufficient to plot the poles with zero or positive imaginary parts. Fig. 1 gives a rough impression of the pole distributions. It is observed that as $T$ increases, the poles of the estimates from both algorithms cluster closer around the true poles. However, the LS estimates have poles outside the unit circle while all poles of our FB estimates are inside the unit circle.

We further plot the histogram of pole magnitudes from the $1000$ repeats in Fig. 2. The same conclusions can be drawn as from Fig. 1. At $T=600$, both histograms have peaks at the true pole magnitudes, indicating good performance. However, we detect $25.9\%,6.9\%,0.2\%$ unstable LS estimates for $T=24,216,600$, respectively, while our FB method always guarantees stability.

We introduce two relative estimation errors (which are plotted as percentages) (using Frobenius norms)

$$e_{FB}=\frac{\|\hat{F}_{R}-F\|}{\|F\|}\mbox{ and }e_{LS}=\frac{\|\hat{F}_{RLS}-F\|}{\|F\|}$$

and two prediction errors

$$\epsilon_{FB} =\frac{\|Y_{1}-\hat{F}_{R}Y_{0}\|-\|Y_{1}-\hat{F}_{LS}Y_{0}\|}{\|Y_{1}-\hat{F}_{LS}Y_{0}\|}$$

$$\epsilon_{LS} =\frac{\|Y_{1}-\hat{F}_{RLS}Y_{0}\|-\|Y_{1}-\hat{F}_{LS}Y_{0}\|}{\|Y_{1}-\hat{F}_{LS}Y_{0}\|}$$

Note that $\hat{F}_{LS}$ minimizes the Frobenius prediction error, so $\epsilon_{FB},\epsilon_{LS}>0$. We plot the histograms of the estimate errors $e_{LS},e_{FB}$ in Fig. 3 and of the prediction errors $\epsilon_{LS},\epsilon_{FB}$ in Fig. 4. It is observed that that FB method has competitive performance against LS even in small sample cases, whilst guaranteeing stability.

Figure 3: Estimation errors $e_{LS},e_{FB}$: the red ‘$*$’ marks and the numbers at the right upper corners are the medians and the blue ‘$|$’ marks are the upper and lower quantiles.

Figure 4: Prediction errors $\epsilon_{LS},\epsilon_{FB}$: the red ‘$*$’ marks and the numbers at the right upper corners are the medians and the blue ‘$|$’ marks are the upper and lower quantiles.

6.3 High Dimensional Simulations

While most stability-guaranteed algorithms are computationally intensive (again, these are for full-rank VAR(1)), and scale badly to high dimensional data, the FB method is not and has the same computational complexity as the LS method. We demonstrate this feature here. We also compare the estimation and prediction errors.

Figure 5: Computational times $t_{comp}$ in log scale: $t_{comp}$ for LS and FB are almost identical and $\log(t_{comp})$ is almost linear with $\log n$. The average computational time is about $49.6$s for model order $n=3072$.

We consider $k=1,2,\dotsm,9$, so that $n=12,\dotsm,3072$. We take the time series lengths $T=36n$ and $T=100n$ and simulate $N=50$ realizations for each $n$. We expect that the computational times $t_{comp}$ for LS and FB be similar and that $\log t_{comp}$ be linear with $\log n$.

We plot the computational times against the model order $n$ in Fig. 5 in log scale. We only plot the case of $T=100n$ because the plots of $T=36n$ are very similar. It is easily observed that FB is as efficient as LS. For the largest model order $n=3072$, the average computational time is only about $49.6$s. To the best of our knowledge, there exists no (full-rank) stability-guaranteed algorithm that can handle such a high model order with achievable computational resources, including time and memory.

For the estimation accuracy, we compare the estimation errors $e_{LS},e_{FB}$ and the prediction errors $\epsilon_{LS},\epsilon_{FB}$ with $T=36n$ and $T=100n$. The results are plotted in Fig. 6. It is interesting to observe that, first, LS and FB still have very similar errors in the higher-order cases, and second, the error medians converge to constants as $n$ grows for $T$ being a multiple of $n$.

Figure 6: Relative estimation and prediction errors plotted against the model order $n$ on log scale: The solid lines are the medians and the dotted lines are the standard deviations.

6.4 Simulations Summary

From the various comparisons above, we conclude that our FB estimator guarantees stability, has low computational complexity, and has competitive accuracy for both short and long length time series, and for both low and high model orders.

7 Conclusions

In this paper, we developed, for the first time, a simple, closed-form, stability-guaranteed estimator for the reduced-rank VAR(1). It is based on a forwards-backwards least squares criterion.

We also gave an asymptotic analysis showing that the new estimator is consistent and asymptotically efficient. Finally, we showed simulations demonstrating the competitive accuracy and computational efficiency of the new stable estimator.

In the future, we will develop a rank selection criterion for medium and large sample sizes.

8 Appendix: Proofs

We use the following reduced-rank least squares lemma.

Lemma R. [20][Theorem 2.2]

$$\hat{C} =\arg\min_{C_{m\times n}:\operatorname{rank}(C)=m}\operatorname{tr}\{(Z-CX)^{\prime}P^{-1}(Z-CX)\}$$

$$=P^{\frac{1}{2}}\hat{V}_{m}\hat{V}_{m}^{\prime}P^{-\frac{1}{2}}ZX^{\prime}(XX^{\prime})^{-1}.$$

where $\hat{V}_{m}=[\hat{v}_{1},\dotsm,\hat{v}_{m}]$ with $\hat{v}_{r}$ the eigenvector corresponding to the $r$-th largest eigenvalue of the matrix $\hat{R}=P^{-\frac{1}{2}}ZX^{\prime}(XX^{\prime})^{-1}XZ^{\prime}P^{-\frac{1}{2}}$.

Proof of Theorem 4B
We show that $J(F;S_{11})$ has the form specified in Lemma R. First, reorganize $J(F;P)$ to get

$$J(F;P)$$

$$= \operatorname{tr}\{P^{-1}(S_{w,f}+S_{w,b})\}$$

$$= \operatorname{tr}\{P^{-1}(-4S_{10}F^{\prime}+FS_{00}F^{\prime}+FPF^{\prime}P^{-1}S_{11})\},$$

where we dropped a constant term not dependent on $F$. Now, set $P=S_{11}$ to get

$$J(F;S_{11})$$

$$= \operatorname{tr}\{S_{11}^{-1}(-4S_{10}F^{\prime}+FS_{00}F^{\prime}+FS_{11}F^{\prime})\}$$

$$= \operatorname{tr}\{S_{11}^{-1}(-4S_{10}F^{\prime}+F(S_{00}+S_{11})F^{\prime}\}$$

$$= \|S_{11}^{-\frac{1}{2}}(2S_{10}(S_{00}+S_{11})^{-\frac{1}{2}}-F(S_{00}+S_{11})^{\frac{1}{2}})\|^{2},$$

where again we dropped a constant term not dependent on $F$. Now apply Lemma R to get the first quoted result.

We now prove that $Q_{*}=S_{11}^{-1}-\hat{F}_{R}^{\prime}S_{11}^{-1}\hat{F}_{R}$ is positive definite. Then Lyapunov’s theorem gives that $\hat{F}_{R}$ is stable.

We have

$$Q_{*} =(S_{11}^{-1}-\hat{F}_{11}^{\prime}S_{11}^{-1}\hat{F}_{11})+(\hat{F}_{11}^{\prime}S_{11}^{-1}\hat{F}_{11}-\hat{F}_{R}^{\prime}S_{11}^{-1}\hat{F}_{R})$$

$$=S_{11}^{-1}U_{b}S_{11}^{-1}+(\hat{F}_{11}^{\prime}S_{11}^{-1}\hat{F}_{11}-\hat{F}_{R}^{\prime}S_{11}^{-1}\hat{F}_{R})$$

$$U_{b} =S_{11}-\hat{F}_{b}S_{11}\hat{F}_{b}^{\prime}$$

$$\hat{F}_{b} =S_{11}\hat{F}_{11}^{\prime}S_{11}^{-1}$$

For the second term we have

$$\hat{F}_{11}^{\prime}S_{11}^{-1}\hat{F}_{11}-\hat{F}_{R}^{\prime}S_{11}^{-1}\hat{F}_{R}$$

$$= \hat{F}_{11}^{\prime}S_{11}^{-\frac{1}{2}}(I-\hat{V}_{m}\hat{V}_{m}^{\prime})S_{11}^{-\frac{1}{2}}\hat{F}_{11}\geq 0.$$

For the first term, we have

$$\begin{array}[]{rrcl}&2S_{10}&=&\hat{F}_{11}S_{00}+\hat{F}_{11}S_{11}\\ \Rightarrow&2S_{01}&=&S_{00}\hat{F}_{11}^{\prime}+S_{11}\hat{F}_{11}^{\prime}\\ \Rightarrow&2S_{01}&=&S_{00}S_{11}^{-1}\hat{F}_{b}S_{11}+\hat{F}_{b}S_{11}\\ \Rightarrow&-S_{00}S_{11}^{-1}\hat{F}_{b}S_{11}\hat{F}_{b}^{\prime}&=&-2S_{01}\hat{F}_{b}^{\prime}+\hat{F}_{b}S_{11}\hat{F}_{b}^{\prime}\\ \Rightarrow&-\Phi U_{b}&=&S_{00}-2S_{01}\hat{F}_{b}+\hat{F}_{b}S_{11}\hat{F}_{b}^{\prime}.\end{array}$$

where $\Phi=-S_{00}S_{11}^{-1}$. Now take the transpose and add the two equations together to get

$$-\Phi U_{b}+U_{b}(-\Phi)^{\prime}= 2[S_{00}-S_{01}\hat{F}_{b}^{\prime}-\hat{F}_{b}S_{10}+\hat{F}_{b}S_{11}\hat{F}_{b}^{\prime}]$$

$$= \frac{2}{T}(Y_{0}-\hat{F}_{b}Y_{1})(Y_{0}-\hat{F}_{b}Y_{1})^{\prime}=2\hat{S}_{w,b}$$

which we show below is positive definite. But note that this is a continuous time Lyapunov equation where $\Phi$ is a stability matrix - it has the same eigenvalues as $-S_{00}^{\frac{1}{2}}S_{11}^{-1}S_{00}^{\frac{1}{2}}$ which is negative definite with, therefore, negative real eigenvalues. Thus $U_{b}$ must be positive definite and the result is established.

To see that $\hat{S}_{w,b}$ is positive definite first note that from A1 and Theorem 2 $\Rightarrow Q_{b}$ is positive definite. Now the backward least squares estimator $\hat{F}_{b,LS}=S_{01}S_{11}^{-1}$ minimises $S_{w,b}(F)$. Then for any $F$, completion of squares gives

$$S_{w,b}(F)=\hat{Q}_{b,LS}+(\hat{F}_{b,LS}-F)S_{11}(\hat{F}_{b,LS}-F)^{\prime}$$

where $\hat{Q}_{b,LS}=S_{w,b}(\hat{F}_{b,LS})$. Set $F=0$ to get $\hat{Q}_{b,LS}=S_{00}-S_{01}S_{11}^{-1}S_{10}$ which $\xrightarrow{p}Q_{b}$. Now since $Q_{b}$ is full rank and the eigenvalues of $\hat{Q}_{b,LS}$ are continuous functions of the entries in $\hat{Q}_{b,LS}$ they cannot accumulate mass at $0$ i.e. the smallest eigenvalue of $\hat{Q}_{b,LS}$ is $0$ with zero probability. So $\hat{Q}_{b,LS}$ is positive definite w.p.1. But setting $F=\hat{F}_{b}$ above, we see that $\hat{S}_{w,b}\geq\hat{Q}_{b,LS}$ and so is positive definite w.p.1 as required.

Proof of Theorem 5E(i).

We have

$$\begin{array}[]{rrcl}&&&\sqrt{T}(\hat{F}_{RLS}-\hat{F}_{R})\\ &&=&\sqrt{T}S_{11}^{\frac{1}{2}}(\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}S_{11}^{-\frac{1}{2}}\hat{F}_{LS}-\hat{V}_{m}\hat{V}_{m}^{\prime}S_{11}^{-\frac{1}{2}}\hat{F}_{11})\\ &&=&S_{11}^{\frac{1}{2}}\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}S_{11}^{-\frac{1}{2}}(\hat{F}_{LS}-\hat{F}_{11})\sqrt{T}\\ &&+&S_{11}^{\frac{1}{2}}\sqrt{T}(\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}-\hat{V}_{m}\hat{V}_{m}^{\prime})S_{11}^{-\frac{1}{2}}\hat{F}_{11}\end{array}$$

From Theorem A, Lemma 5D, and the fact that $\|\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}\|^{2}=m,$ we get the first term $\xrightarrow{p}0$. For the second term, using Theorem 5A, Theorem 3B we are reduced to showing $\sqrt{T}(\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}-\hat{V}_{m}\hat{V}_{m}^{\prime})\xrightarrow{p}0$. Note that as in the proof of Theorem 5A we deduce that $\hat{V}_{m}\xrightarrow{p}V_{m}$.

Next set $\bar{S}=\frac{1}{2}(S_{00}+S_{11})$ and consider that

$$\begin{array}[]{rrcl}&&&\sqrt{T}(\hat{R}-\hat{R}_{*})\\ &&=&\sqrt{T}S_{11}^{-\frac{1}{2}}S_{10}[\bar{S}^{-1}-S_{00}^{-1}]S_{10}^{\prime}S_{11}^{-\frac{1}{2}}\\ &&=&\sqrt{T}S_{11}^{-\frac{1}{2}}S_{10}S_{00}^{-1}(S_{00}-\bar{S})\bar{S}^{-1}S_{10}^{\prime}S_{11}^{\frac{1}{2}}\\ &&=&\sqrt{T}S_{11}^{-\frac{1}{2}}S_{10}S_{00}^{-1}\frac{1}{2}(S_{00}-S_{11})\bar{S}^{-1}S_{10}^{\prime}S_{11}^{-\frac{1}{2}}\end{array}$$

and we find $\sqrt{T}(\hat{R}-\hat{R}_{*})\xrightarrow{p}0$ in view of Theorem 5A and lemma 5C.

Consider the ‘k-th’ eigenvector $\hat{v}_{k}$ and eigenvalue $\hat{\lambda}_{k}$ of $\hat{R}$ as well as the corresponding $\hat{v}_{*,k}$ and $\hat{\lambda}_{*,k}$ of $\hat{R}_{*}$.

Return now to $\sqrt{T}(\hat{V}_{*,m}\hat{V}_{*,m}^{\prime}-\hat{V}_{m}\hat{V}_{m}^{\prime})\xrightarrow{p}0$. We rewrite this as $\sum_{k=1}^{m}\sqrt{T}(\hat{v}_{k}\hat{v}_{k}^{\prime}-\hat{v}_{*,k}\hat{v}_{*,k}^{\prime})\xrightarrow{p}0$ which holds if $\sqrt{T}(\hat{v}_{k}\hat{v}_{k}^{\prime}-\hat{v}_{*,k}\hat{v}_{*,k}^{\prime})\xrightarrow{p}0$ for each $k$. This has Frobenius norm $\sqrt{T}\sqrt{2[1-(\hat{v}_{k}^{\prime}\hat{v}_{*,k})^{2}]}$. If we set $\cos(\theta_{k})=\hat{v}_{k}^{\prime}\hat{v}_{*,k}$ then the norm is $\sqrt{2T}|\sin(\theta_{k})|$.

This allows us to apply a special case of a classic result of [7] which is given in [29][the equation below their equation (1)] namely

$$|\sin(\theta_{k})| \leq\frac{1}{\hat{\delta}_{k}}\|\hat{R}-\hat{R}_{*}\|$$

$$\hat{\delta}_{k} =min(|\hat{\lambda}_{k-1}-\hat{\lambda}_{\ast,k}|,|\hat{\lambda}_{k+1}-\hat{\lambda}_{\ast,k}|)$$

Since (Theorem 5A) $\hat{\lambda}_{l}\xrightarrow{p}\lambda_{l},\hat{\lambda}_{\ast,l}\xrightarrow{p}\lambda_{l}$ for all $l$ then $\hat{\delta}_{k}\xrightarrow{p}\delta_{k}=min(|\lambda_{k-1}-\lambda_{k}|,|\lambda_{k+1}-\lambda_{k}|)$.

However, consider that

$$\sqrt{T}\|\hat{R}-\hat{R}_{*}\|=\sqrt{\operatorname{tr}[\sqrt{T}(\hat{R}-\hat{R}_{*})\sqrt{T}(\hat{R}-\hat{R}_{*})]}$$

but each term $\xrightarrow{p}0$ and so $\sqrt{T}|\sin(\theta_{k})|\xrightarrow{p}0$ and the proof is complete.