1.1 Notation

To describe the use of multivariate techniques, we need to introduce new notation. In its basic form, a VAR consists of a set of \(K\) endogenous variables \({\bf{y}}_t = (y_{1t}, \ldots, y_{kt}, \ldots, y_{Kt})\) for \(k = 1, \ldots K\). After including \(p\) lags of the endogenous variables, the VAR(\(p\)) model may be defined as:

\[\begin{eqnarray} {\bf{y}}_t = A_1 {\bf{y}}_{t-1} + \ldots + A_p {\bf{y}}_{t-p} + CD_t + {\bf{u}}_t \tag{1.1} \end{eqnarray}\] where \(A_i\) are \((K \times K)\) coefficient matrices for \(i = 1, \ldots, p\) and \({\bf{u}}_t\) is a \(K\)-dimensional white noise process with time-invariant positive definite covariance matrix,⁵ \(\mathbb{E} \left[ {\bf{u}}_t {\bf{u}}_t^\prime \right] = \Sigma_{{\bf{u}}}\), such that

\[\begin{eqnarray} \nonumber \mathbb{E} \left[ {\bf{u}}_t \right] &=&0 \\ \nonumber \mathbb{E} \left[ {\bf{u}}_t {\bf{u}}_t^\prime \right] &=& \Sigma_{{\bf{u}}} \; \mathsf{which} \; \mathsf{is} \; \mathsf{positive} \; \mathsf{definite} \end{eqnarray}\]

The matrix \(C\) is the coefficient matrix of potentially deterministic regressors with dimension \((K \times M)\), and \(D_t\) is an \((M \times 1)\) column vector holding the appropriate deterministic regressors, such as a constant, trend, and dummy and/or seasonal dummy variables.

The expression in equation (1.1) clearly suggests that the VAR model relates the \(k\)’th variable in vector \({\bf{y}}_{t}\) to past values of itself and all other variables in the system. This is in contrast with the univariate autoregressive models, which only relate a given variable to past values of itself.

To consider an example, let us assume that there are no deterministic regressors, \(K=2\) and \(P=1\). Then equation (1.1) would simplify to that of a VAR(1) model that does not include any deterministic components,

\[\begin{eqnarray}\nonumber {\bf{y}}_{t}= A_{1} {\bf{y}}_{t-1} + {\bf{u}}_{t} \end{eqnarray}\] where \({\bf{y}}_{t}, {\bf{u}}_{t}\) and the \(A_{1}\) matrices that take the form

\[\begin{eqnarray} {\bf{y}}_{t}=\left[ \begin{array}{c} y_{1,t} \\ y_{2,t}% \end{array}% \right] , A_{1}=\left[ \begin{array}{cc} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22}% \end{array}% \right] \; \mathsf{ and } \; {\bf{u}}_{t}=\left[ \begin{array}{c} u_{1,t} \\ u_{2,t}% \end{array}% \right] \tag{1.2} \end{eqnarray}\]

If we are to then assume that the specific elements of the \(A_{1}\) matrix are given as,

\[\begin{eqnarray} \left[ \begin{array}{c} y_{1,t} \\ y_{2,t}% \end{array}% \right] = \left[ \begin{array}{cc} 0.5 & 0 \\ 1 & 0.2% \end{array}% \right] \left[ \begin{array}{c} y_{1,t-1} \\ y_{2,t-1}% \end{array}% \right] +\left[ \begin{array}{c} u_{1,t} \\ u_{2,t}% \end{array}% \right] \tag{1.3} \end{eqnarray}\] where after some matrix multiplications, we can write

\[\begin{eqnarray}\nonumber y_{1,t} &=& 0.5y_{1,t-1}+u_{1,t} \\ y_{2,t} &=& 1y_{1,t-1}+0.2y_{2,t-1}+u_{2,t} \nonumber \end{eqnarray}\]

With this parametrisation, \(y_{2,t}\) depends positively on past values of itself and on past values of \(y_{1,t}\). On the other hand, the dynamics of \(y_{1,t}\) depend on past values of itself only. These relationships are captured by the particular parametrisation of the model.

The variables that are to be included in a particular investigation will typically depend on the purpose of the study, where researchers may wish to incorporate those that are responsible for important dynamic interactions or a perceived causal relationships. In addition, economic theory could also be used to good effect when deciding upon the variables that are to be included. In many macroeconomic investigations it would be prudent to construct a model that does not only depend on past values of itself, but also on the past values of other variables. Hence, many VAR models typically include lags of all endogenous variables.

1.2 The companion-form

It is possible to express a VAR(\(p\)) model in a VAR(1) form, which turns out to be very useful in many practical situations (and for some technical derivations). This reformulation is done by expressing the VAR(\(p\)) model in the companion-form, which is accomplished as follows. Consider the model,

\[\begin{eqnarray}\nonumber Z_{t}=\Gamma_{0}+\Gamma_{1}Z_{t-1}+ \Upsilon_{t} \end{eqnarray}\]

where we have defined

\[\begin{eqnarray}\nonumber Z_{t}=\left[ \begin{array}{c} {\bf{y}}_{t} \\ {\bf{y}}_{t-1} \\ \vdots \\ {\bf{y}}_{t-p+1} \end{array} \right] , \hspace{1cm} \Gamma_0=\left[ \begin{array}{c} \mu \\ 0 \\ \vdots \\ 0 \end{array} \right] , \hspace{1cm} \Upsilon_{t} =\left[ \begin{array}{c} {\bf{u}}_{t} \\ 0 \\ \vdots \\ 0 \end{array} \right] \end{eqnarray}\]

So that the matrix notation is

\[\begin{eqnarray}\nonumber \left[ \begin{array}{c} {\bf{y}}_{t} \\ {\bf{y}}_{t-1} \\ {\bf{y}}_{t-2} \\ \vdots \\ {\bf{y}}_{t-p+1} \end{array} \right] =\left[ \begin{array}{c} \mu \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right] +\left[ \begin{array}{ccccccc} A_{1} & A_{2} & \cdots & A_{p-1} & A_{p} \\ I & 0 & \cdots & 0 & 0 \\ 0 & I & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I & 0 \end{array} \right] \left[ \begin{array}{c} {\bf{y}}_{t-1} \\ {\bf{y}}_{t-2} \\ {\bf{y}}_{t-3} \\ \vdots \\ {\bf{y}}_{t-p} \end{array} \right] + \left[ \begin{array}{c} {\bf{u}}_{t} \\ 0 \\ 0 \\ \vdots \\ 0 \end{array} \right] \end{eqnarray}\]

The dimensions of the vectors \(Z_{t}\), \(\Gamma_{0}\) and \(\Upsilon_{t}\) are of dimension \(Kp\times 1\). The coefficient matrix, \(A_{j}\) for \(j=1,\ldots, p\) will be of dimension \(K\times K\), and \(\Gamma_{1}\) is \(Kp\times Kp\). In this case, \(\Gamma_{1}\) is called the companion-form matrix.⁶

1.3 Stability in a VAR model

As in the case of the univariate autoregressive model, the VAR model is covariance-stationary when the effect of the shocks, \({\bf{u}}_t\), dissipate over time. This will be the case if the eigenvalues of the companion-form matrix are all less than one in absolute value. The eigenvalues of the matrix \(\Gamma_1\) are represented by \(\lambda\) in the expression,

\[\begin{eqnarray} |\Gamma_{1}-\lambda I|=0 \tag{1.4} \end{eqnarray}\]

To show how one could derive these eigenvalues, we make use of the simple VAR(1) example that we employed in (1.2), where the \(A_{1}\) coefficient matrix is represented by what is contained in (1.3). In this case, \(\Gamma_{1} = A_{1}\), which may be inserted into expression (1.4) to provide,

\[\begin{eqnarray} \nonumber \det \left[ \left[ \begin{array}{cc} 0.5 & 0 \\ 1 & 0.2% \end{array}% \right] -\lambda \left[ \begin{array}{cc} 1 & 0 \\ 0 & 1% \end{array}% \right] \right] &=&\det \left[ \left[ \begin{array}{cc} 0.5-\lambda & 0 \\ 1 & 0.2-\lambda \end{array}% \right] \right] (0.5-\lambda )(0.2-\lambda ) &=&0 \nonumber \end{eqnarray}\]

Such that,

\[\begin{eqnarray} \nonumber \lambda_{1} &=&0.5,\mathsf{ \ }\lambda_{2}=0.2 \notag \end{eqnarray}\]

Hence, the eigenvalues are \(\lambda_{1}=0.5\) and \(\lambda_{2}=0.2\), which are both less than \(1\), which would suggest that this VAR(1) process is stable (i.e. covariance-stationary). In the more general cases of VAR(\(p\)) models, the eigenvalues of \(\Gamma_{1}\) may be easily computed with the aid of most standard software.

Some researchers consider the stability of a model by making use of the definition for characteristic roots rather than eigenvalues. In such cases, we could find the roots, \(z\), of a matrix \(\Gamma_{1}\), from

\[\begin{eqnarray}\nonumber |I-\Gamma_{1}z| \end{eqnarray}\]

When using this formulation, the interpretation is reversed, such that one would suggest that the VAR model is stable when the roots of the polynomial, det\((I-\Gamma_{1}z)\), lie outside of the unit circle (i.e. are bigger than one in absolute value).⁷ In what follows, we will seek to refer to the eigenvalues when investigating the stability of the variables in the model.

To simulate processes that may be described by a stable VAR(1) model, we could make use of equation (1.3), with \(y_{k,0}=0\), \(\mu_{k}=1\) for \(k=[1,2]\) and

\[\begin{eqnarray} \nonumber {\bf{u}}_t \sim \mathcal{N}\left( \left[ \begin{array}{c} 0 \\ 0% \end{array}% \right] ,\left[ \begin{array}{cc} 1 & 0.2 \\ 0.2 & 1% \end{array}% \right] \right) \end{eqnarray}\]

The result of this simulation is shown in Figure 1. As the process is stable by definition, both series \(y_{1,t}\) and \(y_{2,t}\) fluctuate around a constant mean. Note also that their variability does not change over time.

1.4 Wold decomposition

Just as the stable AR(\(p\)) model has a MA representation, the stable VAR(\(p\)) has a VMA representation. This is proposition is provided by the Wold (1938) decomposition theorem, which states that every covariance-stationary time series \(y_{t}\) can be written as the sum of two uncorrelated processes, a deterministic component (which could be the mean) \(\kappa_{t}\) and a component that has an infinite moving average representation of \(\sum_{j=0}^{\infty }\theta_{j}\varepsilon_{t-j}\)

\[\begin{eqnarray} y_{t}=\sum_{j=0}^{\infty }\theta_{j}\varepsilon_{t-j}+\kappa_{t} \tag{1.5} \end{eqnarray}\]

where we assume that \(\theta_{0}=1\), \(\sum_{j=0}^{\infty }|\theta_{j}| <\infty\) and the term \(\varepsilon_{t}\) is Gaussian white noise.

Finding the Wold representation involves fitting an infinite number of parameters \(( \theta_{1}, \theta_{2}, \theta_{3}, \ldots )\) to the data. With a finite number of observations, this is never possible. However, one can approximate an infinite order \(\theta (L)\), by using models that have a finite number of parameters.

We already have implicitly used the Wold’s decomposition theorem when deriving the moving average form of an AR(\(p\)) model. In what follows we will now derive the moving average form of a stable VAR(\(p\)).

2.1 Finding the MA coefficients

To identify the moving average coefficients, \(B_{j}\), we could make use of the relationship \(I=B(L)A(L)\), which would hold due to the result in (2.5).⁸ This allows for the following derivation,

\[\begin{eqnarray} \nonumber I&=&B(L)A(L) \\ \nonumber I &=&(B_{0}+B_{1}L+B_{2}L^{2}+ \ldots )(I-A_{1}L-A_{2}L^{2}- \ldots - A_{p}L^{p}) \\ \nonumber &=&[B_{0}+B_{1}L+B_{2}L^{2}+ \ldots ] \\ \nonumber && -[B_{0}A_{1}L+B_{1}A_{1}L^{2}+B_{2}A_{1}L^{3}+ \ldots ] \\ \nonumber && -[B_{0}A_{2}L^{2}+B_{1}A_{2}L^{3}+B_{2}A_{2}L^{4}+ \ldots] - \ldots \\ \nonumber &=&B_{0}+(B_{1}-B_{0}A_{1})L+(B_{2}-B_{1}A_{1}-B_{0}A_{2})L^{2}+ \ldots \\ \nonumber && +\left(B_{p}-\sum_{j=1}^{p}B_{p-j}A_{j}\right) L^{p}+ \ldots \end{eqnarray}\]

Solving for the relevant lags (noting that \(A_{1}=0\) for \(j>p\)), we get,

\[\begin{eqnarray} \notag B_{0} &=&I \\ B_{1} &=&B_{0}A_{1} \notag \\ B_{2} &=&B_{1}A_{1}+B_{0}A_{2} \notag \\ \vdots & & \vdots \notag \\ B_{i} &=&\sum_{j=1}^{i}B_{i-j}A_{j}\mathsf{ \ \ \ for }i=1,2, \ldots \tag{2.6} \end{eqnarray}\]

Hence, the \(B_{j}\)’s can be computed recursively from (2.6). To show how this may be applied in a practical example, consider the simple VAR(1) example in equation (2.1),

\[\begin{eqnarray*} {\bf{y}}_{t} &=&\mu +A_{1}{\bf{y}}_{t-1}+{\bf{u}}_{t} \\ \end{eqnarray*}\]

Since we know from (2.6) that the coefficients will take the form,

\[\begin{eqnarray*} B_{0} &=&I \\ B_{1} &=&B_{0}A_{1}=A_{1} \\ B_{2} &=&B_{1}A_{1}+B_{0}A_{2}=A_{1}^{2} \\ \vdots && \vdots \\ B_{p} &=&A_{1}^{p} \end{eqnarray*}\]

We are then able to plug in the known parameters from (1.3), to get the moving average coefficients,

\[\begin{eqnarray*} B_{0}&=&I \\ B_{1}&=&A_{1}=\left[ \begin{array}{cc} 0.5 & 0 \\ 1 & 0.2% \end{array} \right] \\ B_{2}&=&A_{1}^{2}=\left[ \begin{array}{cc} 0.25 & 0 \\ 0.70 & 0.04 \end{array} \right] \\ \vdots && \vdots \end{eqnarray*}\]

In this way, \(B_{i}\) approaches zero as \(p \rightarrow \infty\), which would arise if the processes in the VAR(\(p\)) model are stable (i.e. they do not have infinite memory).

2.2 Mean, variance and autocovariance of a VAR

Just as in the univariate case, the first-and second-order moments of the VAR(\(p\)) can easily be derived from the moving average representation in equation (2.4). Hence, given

\[\begin{eqnarray*} {\bf{y}}_{t} &=& \varphi+\sum_{j=0}^{\infty }B_{j} {\bf{u}}_{t-j} \end{eqnarray*}\]

with \(\varphi=\left( \sum_{j=0}^{\infty }B_{j}\right) \mu =(I-A_{1})^{-1}\mu\). Firstly note that, due to Gaussian white noise assumption about the error terms, the expected mean value of the process is simply,

\[\begin{eqnarray*} \mathbb{E} \left[{\bf{y}}_{t}\right]= \varphi =(I-A_{1})^{-1}\mu \end{eqnarray*}\]

Sometimes this mean is described as the steady-state of the system. We will denote the covariance and autocovariances of the process by \(\Psi\), which will be derived with the aid of Yule-Walker equations.

To complete this procedure we firstly write the VAR(\(p\)) in the mean-adjusted form, where we subtract the expected value from both sides of the equality sign,

\[\begin{eqnarray*} {\bf{y}}_{t}-\varphi = A_{1}({\bf{y}}_{t-1}-\varphi)+{\bf{u}}_{t} \end{eqnarray*}\]

Postmultiplying by \(({\bf{y}}_{t-s}-\varphi)^{\prime }\) and taking the expectation gives,

\[\begin{eqnarray*} \mathbb{E} \left[ \left({\bf{y}}_{t}-\varphi \right) \left( {\bf{y}}_{t-s}-\varphi \right)^\prime \right] &=& A_{1} \mathbb{E} \left[ \left({\bf{y}}_{t-1}-\varphi \right) \left( {\bf{y}}_{t-s}-\varphi \right)^\prime \right] \\ &&+ \; \mathbb{E} \left[{\bf{u}}_{t} \left({\bf{y}}_{t-s}-\varphi \right)^{\prime }\right] \end{eqnarray*}\]

Thus for \(s=0\),

\[\begin{eqnarray} \Psi_{s}=A_{1} \Psi_{-1} + \Sigma_{{\bf{u}}}= A_{1} \Psi_{1}^{\prime}+ \Sigma_{{\bf{u}}} \tag{2.7} \end{eqnarray}\]

where after the second equality sign, we have used the fact that \(\Psi_{-1}=\Psi_{1}^{\prime }\). For \(s>0\) we would then have,

\[\begin{eqnarray} \Psi_{s}=A_{1}\Psi_{s-1} \tag{2.8} \end{eqnarray}\]

When \(A_{1}\) and \(\Sigma_{{\bf{u}}}\) are known, we can compute the autocovariance for \(s=0,\ldots, S\) using equations (2.7) and (2.8). Hence, for \(s=1\), we use (2.8) to show that \(\Psi_{1} = A_{1} \Psi_{0}\). Substituting this into equation (2.7) and noting that \((A_{1}\Psi_{0})^{\prime }= \Psi_{0}^{\prime } A_{1}^{\prime }\) by the rules of matrix algebra, we get,

\[\begin{eqnarray*} \Psi_{0}=A_{1}\Psi_{0} A_{1}^{\prime} + \Sigma_{{\bf{u}}} \end{eqnarray*}\]

This equation can be solved for \(\Psi_{0}\), by using the definition of the Kronecker product and the \(vec\) operator to derive,⁹

\[\begin{eqnarray*} vec \Psi_{0} =(I-A_{1}\otimes A_{1})^{-1}vec\Sigma_{{\bf{u}}} \end{eqnarray*}\]

where \(vec\) denotes the column stacking operator that stacks the columns of a matrix in a column vector. It follows that once \(\Psi_{0}\) has been derived in this manner, we can use equation (2.8) recursively to get the autocovariance for \(s>0\).

To get the autocorrelation function, we would then need to normalise the autocovariances such that they have ones on the diagonal at \(s=0.\) Thus, we define the diagonal matrix, \(\vartheta\), whose diagonal elements are the square roots of the diagonal elements of \(\Psi_{0}\). Then the autocorrelation function for the VAR(\(p\)) is simply,

\[\begin{eqnarray*} R_s=\vartheta^{-1}\Psi_{h} \vartheta^{-1} \end{eqnarray*}\]

To keep the notations as simple as possible, the above derivation assumed that the underlying model was a VAR(1). However, as noted earlier, a VAR(\(p\)) can always be written as a VAR(1) using companion-form of matrix. For the VAR(1), the companion-form matrix is of course equal to the \(A_{1}\) matrix, and nothing in the previous derivations changes. For higher order VAR(\(p\)) models the companion-form, all of the expressions for the moving average representation, expected value and autocovariances do not change; but the size of the matrices will of course change. In addition, we would also need to make use a selection matrix to extract the values of interest. For a more detailed explanation of how this is done, see Lutkepohl (2005).

3.1 Choice of variables and lags

When deciding on the variables and lags that should be included in the model, it is important to note that VAR can become heavily parameterised. This could result in problems where we have too many parameters, relative to observations in the data (i.e. degrees of freedom problem). For example, after including 6 variables in a model with 8 lags (i.e. two years of lags for quarterly data), then each equation in the VAR will contain (6 \(\times\) 8)+1=49 coefficients (including a constant). In many macroeconomic or financial investigations, we often find ourself in situations with limited observations, which would imply that our parameter estimates are imprecise (at best). Therefore, in most cases, including much more than 6 variables in the VAR that is estimated with classical techniques could result in several difficulties.

Therefore, since one cannot include all variables of potential interest, one has to make a careful choice by referring to economic theory or a prior knowledge when choosing variables. Having specified the model, the appropriate lag length of the VAR model also has to be chosen. As in the univariate case, one can use information criterion for this purpose (which would include the BIC or AIC).

4.1 Point forecast

To generate the \(h\)-step ahead forecast, we could iterate the VAR model forward by a number of observations. Replacing the notation used for the univariate case with the VAR notation, the multivariate expression is,

\[\begin{eqnarray} {\bf{y}}_{t+h}=A_{1}^{h} {\bf{y}}_{t}+\sum_{i=0}^{h-1}A_{1}^{i} {\bf{u}}_{t+h-i} \tag{4.1} \end{eqnarray}\]

While this forecasting framework has been described within the context of a VAR(\(1\)) model, we could express the VAR(\(p\)) model as a VAR(\(1\)) by using the companion-form, where the formula would remain largely unchanged. That is, \(Z_{t+h}=\Gamma_{1}^{h}Z_{t}+\sum_{i=0}^{h-1}\Gamma_{1}^{i} \Upsilon_{t+h-i}\) . Using the selector matrix \(\Xi=[I,0, \ldots 0]\), which is of \(K\times (Kp)\) dimension, the appropriate forecasts for \({\bf{y}}_{t+h}\) can be extracted (i.e.\({\bf{y}}_{t+h}=\Xi Z_{t+h}\)). Moreover, adding a constant to the forecasting relationship is simply a multivariate extension of the univariate case. Therefore, continuing with the VAR(1) example, but adding a constant to the VAR, results in the following \(h\)-step ahead forecast formula,

\[\begin{eqnarray} {\bf{y}}_{t+h}=(I+A_{1}+ \ldots A_{1}^{h-1})\mu +A_{1}^{h}{\bf{y}}_{t}+\sum_{i=0}^{h-1}A_{1}^{i} {\bf{u}}_{t+h-i} \tag{4.2} \end{eqnarray}\]

By employing the conditional expectation to (4.2), we get the VAR point forecast,

\[\begin{eqnarray*} \mathbb{E}\left[{\bf{y}}_{t+h}|{\bf{y}}_{t}\right] = (I+A_{1}+ \ldots +A_{1}^{h-1})\mu +A_{1}^{h}{\bf{y}}_{t} \end{eqnarray*}\]

Which under the maintained assumption about stationarity, converges to the unconditional mean of the process when \(h\) goes to infinity,

\[\begin{eqnarray*} \mathbb{E}\left[{\bf{y}}_{t+h}|{\bf{y}}_{t}\right]=\frac{\mu }{I-A_{1}} \hspace{1cm}\mathsf{when} \; h\rightarrow \infty \end{eqnarray*}\]

Once again, these equations are just the multivariate extensions of the formulas that were derived previously.

4.2 Mean Squared Errors

The Mean Squared Errors (MSE) are derived from equation (4.2). Using the notation from the previous exposition, let \(\acute{{\bf{y}}}(h)=\mathbb{E}\left[{\bf{y}}_{t+h}|{\bf{y}}_{t}\right]\) denote the predictor. Then, the VAR forecast error at horizon \(h\) is,

\[\begin{eqnarray*} {\bf{y}}_{t+h} - \acute{{\bf{y}}}(h) =\sum_{i=0}^{h-1}A_{1}^{i} {\bf{u}}_{t+h-i} \end{eqnarray*}\]

As before, the expectation of this expression is the expected forecast error. Under the current assumptions about \(\mathbb{E}\left[{\bf{u}}_t\right]\), this expectation is zero. Hence,

\[\begin{eqnarray*} \mathbb{E}\left[{\bf{y}}_{t+h} - \acute{{\bf{y}}}(h) \right] = \mathbb{E}\left[{\bf{y}}_{t+h}\right] - \mathbb{E}\left[\acute{{\bf{y}}}(h) \right]=0 \end{eqnarray*}\]

Thus the predictor \(\acute{{\bf{y}}}(h)\) is unbiased. The MSE is simply the forecast error variance. Therefore, in the multivariate VAR setting the MSE, which we denote, \(\acute{{\bf{\sigma}}}_{y}(h)\), may be given as

\[\begin{eqnarray*} \acute{{\bf{\sigma}}}_{y}(h)&=& \mathbb{E}\left[\left({\bf{y}}_{t+h} - \acute{{\bf{y}}}(h)\right)^2 \right] \\ &=& \mathbb{E}\left[ \left( \sum_{i=0}^{h-1}A_{1}^{i} {\bf{u}}_{t+h-i}\right) \left( \sum_{i=0}^{h-1}A_{1}^{i} {\bf{u}}_{t+h-i}\right) \right] \end{eqnarray*}\]

By noting that the \(A_{1}\) terms are constant, which would allow us to move them outside of the expectation, and that \(\mathbb{E}({\bf{u}}_{t+h-i},{\bf{u}}_{t+h-i})=\Sigma_{{\bf{u}}}\) for all \(h\), we are able to derive

\[\begin{eqnarray} \acute{{\bf{\sigma}}}_{y}(h)=\sum_{i=0}^{h-1} \left( A_{1}^{j}\Sigma_{{\bf{u}}}A_{1}^{j^{\prime}} \right) \tag{4.3} \end{eqnarray}\]

Given this expression, it is worth noting that there is an interesting relationship between the MSE and the moving average form of the VAR. If we expand equation (2.4), and move it \(h\)-periods forward we get,

\[\begin{eqnarray*} {\bf{y}}_{t+h}= \varphi + {\bf{u}}_{t+h}+B_{1} {\bf{u}}_{t+h-1}+B_{2} {\bf{u}}_{t+h-2}+ \ldots +B_{h} {\bf{u}}_{t}+B_{h+1} {\bf{u}}_{t-1}+ \ldots \\ \end{eqnarray*}\]

Employing the conditional expectation to this expression, and computing the forecast error gives,

\[\begin{eqnarray*} \mathbb{E} [{\bf{y}}_{t+h}|t]=\varphi+B_{h}{\bf{u}}_{t}+B_{h+1}{\bf{u}}_{t-1}+B_{h+2}{\bf{u}}_{t-2}+ \ldots \end{eqnarray*}\]

and

\[\begin{eqnarray*} {\bf{y}}_{t+h}- \mathbb{E} [{\bf{y}}_{t+h}|t]={\bf{u}}_{t+h}+B_{1}{\bf{u}}_{t+h-1}+ \ldots +B_{h-1}{\bf{u}}_{t+1} \end{eqnarray*}\]

Which has variance equal to,

\[\begin{eqnarray*} \acute{{\bf{\sigma}}}_{y}(h)=\mathbb{E}\left[ \left( {\bf{y}}_{t+h}-\mathbb{E}[{\bf{y}}_{t+h}|t] \right) \left( {\bf{y}}_{t+h}-\mathbb{E}[{\bf{y}}_{t+h}|t] \right)^{\prime}\right] =\sum_{j=0}^{h-1}\left( B_{j}\Sigma_{{\bf{u}}} {B}_{j}^{\prime} \right) \end{eqnarray*}\]

This expression is of course equal to (4.3), but now stated using the moving average coefficients. Moreover, when \(h \rightarrow \infty\) it is equal to \(\Psi (0)\), the unconditional variance of the VAR process. In mathematical terms,

\[\begin{eqnarray*} \acute{{\bf{\sigma}}}_{y}(h)=\Psi (0) \hspace{1cm} \; \mathsf{when} \; h\rightarrow \infty \end{eqnarray*}\]

Thus just as in the univariate case, the forecast variance converges to the unconditional variance of the process.

4.3 Uncertainty

When we assume that the errors of the VAR model are Gaussian and independent across time ( i.e. \({\bf{u}}_{t}\sim \mathsf{i.i.d} \mathcal{N}(0,\Sigma_{{\bf{u}}})\) ), it implies that the forecast errors for the \(k\) variables in the model are also normally distributed, such that,

\[\begin{eqnarray*} \frac{{\bf{y}}_{k,t+h} - \acute{{\bf{y}}_k}(h)}{\acute{{\bf{\sigma}}}_{y,k}(h)}\sim \mathcal{N}(0,1) \end{eqnarray*}\]

Where \({\bf{y}}_{k,t+h}\) and \(\acute{{\bf{y}}_k}(h)\) are the \(k\)’th elements of the actual and predictor values, and \(\acute{{\bf{\sigma}}}_{y,k}(h)\) is the square root of the variance for the \(k\)’th equation. That is, the square root of the \(k\)’th element on the diagonal on the \(\acute{{\bf{\sigma}}}_{y}(h)\) matrix. As such, forecast intervals can be generated around the VAR point forecasts using,

\[\begin{eqnarray*} \big[ \acute{{\bf{y}}_k}(h) - z_{\alpha /2} \; \acute{{\bf{\sigma}}}_{y,k}(h)\; , \; \acute{{\bf{y}}_k}(h) + z_{\alpha /2} \; \acute{{\bf{\sigma}}}_{y,k}(h) \big] \end{eqnarray*}\]

which is equivalent to what was derived in the previous discussion on forecasts

4.4 Forecast failure in macroeconomics

Small-scale VARs are widely used in macroeconomics. Since Sims (1980) pioneering work, examples of VARs used to forecast key economic variables such as output, prices, and the interest rates have been numerous. Yet, a body of recent work suggests that VAR models may be prone to instabilities. In the face of such instabilities, a variety of estimation or forecasting methods might be used to improve the accuracy of forecasts from a VAR. These methods include using an intercept correction, time varying-parameter, differencing the data and model averaging, to mention but a few.¹¹

While many different structural factors could lead to instabilities in macroeconomic variables, breaks are regarded as the primary cause of poor forecast performance. The best way to deal with this problem often varies with the variable being forecasted, but some features emerge. According to Clements and Hendry (2011), forecasting in differences rather the in levels can provide some protection against mean shifts in the dependent variable.

Furthermore, Clark and McCracken (2008) argue that when the VAR in difference is estimated with Bayesian method or using some form of modelling averaging, the forecast consistency performs among the best methods.