# Is there a convenient form for this large covariance matrix?

Consider the following bivariate vector autoregression: $$X_t=\mu +X_{t-1}A+\varepsilon_t,\ \varepsilon_t \overset{iid}{\sim} MVN(0, V),\ X_t=(X_{1,t},X_{2,t})',$$

where the assumptions on the coefficient matrix $A$ are such that the process $\{X_t\}$ is stationary, i.e. there are no unit roots.

The goal is to find a nice (from a computational perspective) expression for the inverse covariance matrix for the vector $y=(X_{1,1},X_{1,2},\dots , X_{1,T}, X_{2,1} \dots X_{2,T})'$.

In the univariate case, i.e. when $\{X_t\}$ is an AR(1) process such as:

$$X_t=\mu + aX_{t-1}+\varepsilon_t,\ \varepsilon_t \sim N(0, v)$$

then Chen and Deo (2009) gives the following convenient expression: $\Omega^{-1}=B'B,$ where

$$B= \ \left( \begin{array}{ccc} \sqrt{1-a^2} & 0 & 0 & \dots & 0\\ -a & 1 & 0 & \dots & 0\\ 0 & -a & 1 &\dots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & -a & 1 \end{array} \right)\$$

and $v\Omega$ is the $T\times T$ covariance matrix for $(X_1,X_2,\dots,X_T)'$.

If it helps, and I think it does, we may assume that $A$ is diagonal so that the bivariate autoregression is really two separate AR(1)-processes with possibly correlated innovations.

The problem is research relate and I need to to speed up numerical computations of bilinear forms such as $w'\Sigma^{-1} u$, for some vectors $w$ and $u$, where $\Sigma$ is the covariance matrix of $y$.

I have a feeling this may be a standard problem solved in some old paper -- I was indeed expecting to find something in a textbook even -- but I can't find any good references. All help, even if just hints and not full answers, is greatly appreciated.

Update In the case with diagonal $A$, all the elements of $\Sigma$ are easily found. We have $$\mathrm{Cov}(X_{i,t},X_{j,t+h})=A_{j,j}^h \frac{\mathrm{Cov}(\varepsilon_{i,t},\varepsilon_{j,t})}{1-A_{i,i}A_{j,j}}=A_{j,j}^h \frac{V_{i,j}}{1-A_{i,i}A_{j,j}}.$$ Note that the fraction only takes 3 different values, one for $i=j=1$, one for $i=j=2$ and one for $i\neq j$. These values corresponds to different blocks of $\Sigma$. Each of these blocks' inverse can be found directly using Chen and Deo's expression (I think).

Another observation is that the power of the $A_{j,j}$ term is just incremented by $\pm 1$ as we move up/down or left/right within these blocks of the matrix.

It looks like that one should assume that the system has already evolved for a sufficiently long time and reached its stationary point. Then, one considers $T-1$ steps of subsequent evolution and the covariance in this interval.

Hence, let us first find the stationary point. Suppose that we started with an initial value $\tilde{X}_{1}$. Iterating the autoregression gives $$\tilde{X}_{t} = (I-A)^{-1}(I-A^{t-1})\mu + A^{t-1}\tilde{X}_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}$$ Assuming that this is a stationary process, large powers of $A$ can be ignored, and we obtain the following result: $$\begin{split} \tilde{X}_{\infty} &= (I-A)^{-1}\mu + \sum_{s=0}^{\infty} A^{s}\varepsilon_{1-s}\\ &\equiv X_{1}. \end{split}$$ This will be our starting point for further evolution. From the fact that $\varepsilon_{t}$ are i.i.d. $\rm{MVN}(0,V)$, we can deduce that $$\rm{E}(X_{1}) = (I-A)^{-1}\mu,$$ $$\rm{Cov}(X_{1}) = \sum_{s=0}^{\infty} A^{s}V A^{\prime\,s} \equiv V_{1}.$$ The covariance $V_{1}$ satisfies the following discrete Lyapunov equation: $$AV_{1}A^{\prime} - V_{1} +V = 0.$$

Looking back at the equation at the top, it is obvious that subsequent evolution from $X_{1}$ is simply given by $$X_{t} = (I-A)^{-1}(I-A^{t-1})\mu + A^{t-1}X_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}$$ As we are ultimately concerned with covariance, considering the deviation makes life easier: $$\delta X_{t}\equiv X_{t} - \rm{E}(X_{t}),$$ whose evolution is as follows: $$\delta X_{t} = A^{t-1}\delta X_{1} +\sum_{s=0}^{t-2} A^{s}\varepsilon_{t-s}.$$ Notice that $\delta X_{1}$, $\epsilon_{2}$, $\epsilon_{3}$, ... are independent (vector) variables of zero mean and covariance $V_{1}$, $V$, $V$, .... That is, the covariance matrix of $\boldsymbol{\varepsilon} \equiv (\delta X_{1},\varepsilon_{2},\ldots,\varepsilon_{T})^{\prime}$ is $$\textbf{V} \equiv \left(\begin{array}{cccccc} V_{1}&0&0&\cdots&0&0\\ 0&V&0&\cdots&0&0\\ 0&0&V&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&V&0\\0&0&0&\cdots&0&V\end{array}\right).$$ We need to turn this into the covariance of $\delta\textbf{X} \equiv (\delta X_{1},\delta X_{2},\ldots,\delta X_{T})^{\prime}$. To do this, we note that $$\delta\textbf{X} = \textbf{M} \boldsymbol{\varepsilon},$$ where $$\textbf{M} \equiv \left(\begin{array}{cccccc} I&0&0&\cdots&0&0\\ A&I&0&\cdots&0&0\\ A^2&A&I&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\A^{T-2}&A^{T-3}&A^{T-4}&\cdots&I&0\\A^{T-1}&A^{T-2}&A^{T-3}&\cdots&A&I\end{array}\right).$$ This follows from the expression for $\delta X_{t}$ given above. Therefore, the covariance matrix of the transformed variable $\delta\textbf{X}$, and hence that of $\textbf{X} \equiv (X_{1},X_{2},\ldots,X_{T})^{\prime}$ itself, is given by $$\textbf{W} = \textbf{M}\textbf{V}\textbf{M}^{\prime}.$$ Its inverse is $$\textbf{W}^{-1} \equiv \textbf{B}^{\prime}\textbf{V}^{-1}\textbf{B},$$ where $$\textbf{B}\equiv\textbf{M}^{-1} = \left(\begin{array}{cccccc} I&0&0&\cdots&0&0\\ -A&I&0&\cdots&0&0\\ 0&-A&I&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&I&0\\0&0&0&\cdots&-A&I\end{array}\right).$$ As a final remark, notice that in the single-variable case, $V_{1}$ is simply given by $$V_{1} = \frac{V}{1-A^{2}},$$ which leads to the formula OP quoted.

Update: To make it easier to work in a rearranged basis, one can explicitly carry out the matrix multiplication in $\textbf{W}^{-1} \equiv \textbf{B}^{\prime}\textbf{V}^{-1}\textbf{B}$. The result is $$\textbf{W}^{-1} = \left(\begin{array}{cccccc} V_{1}^{-1}+P&-Q^{\prime}&0&\cdots&0&0\\ -Q&V^{-1}+P&-Q^{\prime}&\cdots&0&0\\ 0&-Q&V^{-1}+P&\cdots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&V^{-1}+P&-Q^{\prime}\\0&0&0&\cdots&-Q&V^{-1}\end{array}\right),$$ where $P\equiv A^{\prime}V^{-1}A$ and $Q\equiv V^{-1}A$. One can now rearrange the rows and columns of this matrix to transform to a basis that is differently ordered.

• Thanks for the answer! I'm on my phone until this weekend, so I haven't been able to check it out seriously. I'll award the bounty before it runs out as soon as I can check it and unless something even better comes in. – ekvall Aug 28 '14 at 19:18
• @KarlOskar You're welcome. Please let me know if there is anything unclear so that I can improve my answer. – higgsss Aug 29 '14 at 0:24
• It looks right to me, though I'm not very used to working with block matrices. One question: you write that $\boldsymbol W$ is the covariance matrix of $(X_1, \dots , X_T)'$, but is this really the same as the covariance matrix for $(X_{1,1}, \dots X_{1,T}, \dots, X_{2,T})'$ ? – ekvall Aug 29 '14 at 1:28
• @KarlOskar It will be a covariance matrix of $(X_{1,1},X_{2,1}, \ldots, X_{1,T}, X_{2,T})^{\prime}$. In the notation I used, the "element" $(m,n)$ in this covariance matrix really denotes the 2x2 covariance matrix of two vector variables $X_{m}$ and $X_{n}$. – higgsss Aug 29 '14 at 1:42
• Ok! That's what I thought. Obviously, the covariance matrix you derived contains the same information as the one I seek. It may give me a route to speed up computations. – ekvall Aug 29 '14 at 2:04