\chapter{Vector Autoregressions}
\label{chap:var}

Gretl provides a standard set of procedures for dealing with the
multivariate time-series models known as VARs (\emph{Vector
  AutoRegression}). More general models---such as VARMAs, nonlinear
models or multivariate GARCH models---are not provided as of now,
although it is entirely possible to estimate them by writing custom
procedures in the gretl scripting language. In this chapter, we
will briefly review gretl's VAR toolbox.

\section{Notation}
\label{sec:var-def}

A VAR is a structure whose aim is to model the time persistence of a
vector of $n$ time series, $y_t$, via a multivariate autoregression,
as in
\begin{equation}
  \label{eq:VAR}
  y_t = A_1 y_{t-1} + A_2 y_{t-2} + \cdots + A_p y_{t-p} +
  B x_t + \epsilon_t 
\end{equation}
The number of lags $p$ is called the \emph{order} of the VAR. The
vector $x_t$, if present, contains a set of exogenous variables, often
including a constant, possibly with a time trend and seasonal
dummies. The vector $\epsilon_t$ is typically assumed to be a vector
white noise, with covariance matrix $\Sigma$.

Equation \eqref{eq:VAR} can be written more compactly as
\begin{equation}
  \label{eq:VARpoly}
  A(L) y_t = B x_t + \epsilon_t
\end{equation}
where $A(L)$ is a matrix polynomial in the lag operator, or as
\begin{equation}
  \label{eq:VARcompan}
  \left[\begin{array}{c} y_t \\ y_{t-1} \\ \cdots \\ y_{t-p-1} 
    \end{array} \right] = 
  \mathbf{A}
  \left[\begin{array}{c} y_{t-1} \\ y_{t-2} \\ \cdots \\ y_{t-p} 
    \end{array} \right] +
  \left[\begin{array}{c} B \\ 0 \\ \cdots \\ 0
    \end{array} \right] x_t +
  \left[\begin{array}{c} \epsilon_t \\ 0 \\ \cdots \\ 0 
    \end{array} \right]
\end{equation}
The matrix $\mathbf{A}$ is known as the ``companion matrix'' and equals
\[
\mathbf{A} =
  \left[\begin{array}{ccccc} 
      A_1 & A_2 & \cdots & A_p \\ 
      I & 0 & \cdots & 0 \\ 
      0 & I & \cdots & 0 \\ 
      \vdots & \vdots & \ddots & \vdots
    \end{array} \right]
\]
Equation \eqref{eq:VARcompan} is known as the ``companion form'' of
the VAR.

Another representation of interest is the so-called ``VMA
representation'', which is written in terms of an infinite series of
matrices $\Theta_i$ defined as
\begin{equation}
  \label{eq:VMA}
  \Theta_i = \frac{\partial y_t}{\partial \epsilon_{t-i}}
\end{equation}
The $\Theta_i$ matrices may be derived by recursive substitution in
equation \eqref{eq:VAR}: for example, assuming for simplicity that
$B=0$ and $p=1$, equation \eqref{eq:VAR} would become
\[
  y_t = A y_{t-1} + \epsilon_t
\]
which could be rewritten as
\[
  y_t = A^{n+1} y_{t-n-1} + \epsilon_t + A \epsilon_{t-1} + A^2
  \epsilon_{t-2} + \cdots + A^n \epsilon_{t-n}
\]
In this case $\Theta_i = A^i$. In general, it is possible to compute
$\Theta_i$ as the $n \times n$ north-west block of the $i$-th power of the
companion matrix $\mathbf{A}$ (so $\Theta_0$ is always an identity matrix).

The VAR is said to be \emph{stable} if all the eigenvalues of the
companion matrix $\mathbf{A}$ are smaller than 1 in absolute value, or
equivalently, if the matrix polynomial $A(L)$ in equation
\eqref{eq:VARpoly} is such that $|A(z)| = 0$ implies $|z|>1$. If this
is the case, $\lim_{n \to \infty} \Theta_n = 0$ and the vector $y_t$ is
stationary; as a consequence, the equation
\begin{equation}
  \label{eq:VMArep}
  y_t - E(y_t) = \sum_{i=0}^{\infty} \Theta_i \epsilon_{t-i}
\end{equation}
is a legitimate Wold representation. 

If the VAR is not stable, then the inferential procedures that are
called for become somewhat more specialized, except for some simple
cases. In particular, if the number of eigenvalues of $\mathbf{A}$
with modulus 1 is between 1 and $n-1$, the canonical tool to deal with
these models is the cointegrated VAR model, discussed in chapter
\ref{chap:vecm}.

\section{Estimation}
\label{sec:var-estim}

\begin{script}[htbp]
  \scriptcaption{Estimation of a VAR via OLS}
  \label{script:var-ols}
Input:
\begin{scodebit}
open sw_ch14.gdt
series infl = 400*sdiff(log(PUNEW))

scalar p = 2
list X = LHUR infl
list Xlag = lags(p,X)

loop foreach i X
    ols $i const Xlag
endloop

var p X
\end{scodebit}

Output (selected portions):
\begin{outbit}
Model 1: OLS, using observations 1960:3-1999:4 (T = 158)
Dependent variable: LHUR

             coefficient   std. error   t-ratio   p-value 
  --------------------------------------------------------
  const       0.113673     0.0875210     1.299    0.1960  
  LHUR_1      1.54297      0.0680518    22.67     8.78e-51 ***
  LHUR_2     -0.583104     0.0645879    -9.028    7.00e-16 ***
  infl_1      0.0219040    0.00874581    2.505    0.0133   **
  infl_2     -0.0148408    0.00920536   -1.612    0.1090  

Mean dependent var   6.019198   S.D. dependent var   1.502549
Sum squared resid    8.654176   S.E. of regression   0.237830

...


VAR system, lag order 2
OLS estimates, observations 1960:3-1999:4 (T = 158)
Log-likelihood = -322.73663
Determinant of covariance matrix = 0.20382769
AIC = 4.2119
BIC = 4.4057
HQC = 4.2906
Portmanteau test: LB(39) = 226.984, df = 148 [0.0000]

Equation 1: LHUR

             coefficient   std. error   t-ratio   p-value 
  --------------------------------------------------------
  const       0.113673     0.0875210     1.299    0.1960  
  LHUR_1      1.54297      0.0680518    22.67     8.78e-51 ***
  LHUR_2     -0.583104     0.0645879    -9.028    7.00e-16 ***
  infl_1      0.0219040    0.00874581    2.505    0.0133   **
  infl_2     -0.0148408    0.00920536   -1.612    0.1090  

Mean dependent var   6.019198   S.D. dependent var   1.502549
Sum squared resid    8.654176   S.E. of regression   0.237830
\end{outbit}
\end{script}

The gretl command for estimating a VAR is \cmd{var} which, in
the command line interface, is invoked in the following manner:
\begin{flushleft}
    \texttt{[ \emph{modelname} \textless - ] var \emph{p} \emph{Ylist} [;
    \emph{Xlist}]}
\end{flushleft}
where \texttt{p} is a scalar (the VAR order) and \texttt{Ylist} is a
list of variables specifying the content of $y_t$.  The optional
\texttt{Xlist} argument can be used to specify a set of exogenous
variables. If this argument is omitted, the vector $x_t$ is taken to
contain a constant (only); if present, it must be separated from
\texttt{Ylist} by a semicolon. Note, however, that a few common
choices can be obtained in a simpler way: the options \option{trend}
and \option{seasonals} call for inclusion of a linear trend and a set
of seasonal dummies respectively. In addition the \option{nc} option
(no constant) can be used to suppress the standard inclusion of a
constant.

The ``\texttt{\textless -}'' construct can be used to store the model
under a name (see section \ref{sect-script-objects}), if so
desired. To estimate a VAR using the graphical interface, choose
``Time Series, Vector Autoregression'', under the Model menu.

The parameters in eq. \eqref{eq:VAR} are typically free from
restrictions, which implies that multivariate OLS provides a
consistent and asymptotically efficient estimator of all the
parameters.\footnote{In fact, under normality of $\epsilon_t$ OLS is
  indeed the conditional ML estimator. You may want to use other
  methods if you need to estimate a VAR in which some parameters are
  constrained.}  Given the simplicity of OLS, this is what every
software package, including gretl, uses; example script
\ref{script:var-ols} exemplifies the fact that the \cmd{var} command
gives you exactly the output you would have from a battery of OLS
regressions. The advantage of using the dedicated command is that,
after estimation is done, it makes it much easier to access certain
quantities and manage certain tasks. For example, the \dollar{coeff}
accessor returns the estimated coefficients as a matrix with $n$
columns and \dollar{sigma} returns an estimate of the matrix $\Sigma$,
the covariance matrix of $\epsilon_t$. 

Moreover, for each variable in the system an $F$ test is automatically
performed, in which the null hypothesis is that no lags of variable
$j$ are significant in the equation for variable $i$. This is commonly
known as a \textbf{Granger causality} test.

\begin{table}[htbp]
  \centering
  \begin{tabular}{rl}
    \hline
    Periodicity & horizon \\
    \hline
    Quarterly & 20 (5 years) \\
    Monthly & 24 (2 years) \\
    Daily & 3 weeks \\
    All other cases & 10 \\
    \hline
  \end{tabular}
  \caption{VMA horizon as a function of the dataset periodicity}
  \label{tab:var-horizon}
\end{table}

In addition, two accessors become available for the companion matrix
(\dollar{compan}) and the VMA representation (\dollar{vma}). The
latter deserves a detailed description: since the VMA representation
\eqref{eq:VMArep} is of infinite order, gretl defines a
\emph{horizon} up to which the $\Theta_i$ matrices are computed
automatically. By default, this is a function of the periodicity of
the data (see table \ref{tab:var-horizon}), but it can be set by the
user to any desired value via the \cmd{set} command with the
\cmd{horizon} parameter, as in
\begin{code}
set horizon 30
\end{code}
Calling the horizon $h$, the \dollar{vma} accessor returns an $(h+1)
\times n^2$ matrix, in which the $(i+1)$-th row is the vectorized form
of $\Theta_i$.

\subsection{VAR lag-order selection}

In order to help the user choose the most appropriate VAR order,
gretl provides a special variant of the \texttt{var} command:
\begin{flushleft}
  \texttt{var \emph{p} \emph{Ylist} [; \emph{Xlist}]} \option{lagselect}
\end{flushleft}
When the \option{lagselect} option is given, estimation is performed
for all lags up to \texttt{\emph{p}} and a table is printed: it
displays, for each order, a Likelihood Ratio test for the order $p$
versus $p-1$, plus an array of information criteria (see chapter
\ref{chap:criteria}). For each information criterion in the table, a
star indicates what appears to be the ``best'' choice. The same output
can be obtained through the graphical interface via the ``Time Series,
VAR lag selection'' entry under the Model menu.

\begin{script}[htbp]
  \caption{VAR lag selection via Information Criteria}
  \label{script:var-lagselect}
Input:
\begin{scodebit}
open denmark
list Y = 1 2 3 4
var 4 Y --lagselect
var 6 Y --lagselect
\end{scodebit}
% $
Output (selected portions):
\begin{outbit}
VAR system, maximum lag order 4

The asterisks below indicate the best (that is, minimized) values
of the respective information criteria, AIC = Akaike criterion,
BIC = Schwarz Bayesian criterion and HQC = Hannan-Quinn criterion.

lags        loglik    p(LR)       AIC          BIC          HQC

   1     609.15315           -23.104045   -22.346466*  -22.814552 
   2     631.70153  0.00013  -23.360844*  -21.997203   -22.839757*
   3     642.38574  0.16478  -23.152382   -21.182677   -22.399699 
   4     653.22564  0.15383  -22.950025   -20.374257   -21.965748 

VAR system, maximum lag order 6

The asterisks below indicate the best (that is, minimized) values
of the respective information criteria, AIC = Akaike criterion,
BIC = Schwarz Bayesian criterion and HQC = Hannan-Quinn criterion.

lags        loglik    p(LR)       AIC          BIC          HQC

   1     594.38410           -23.444249   -22.672078*  -23.151288*
   2     615.43480  0.00038  -23.650400*  -22.260491   -23.123070 
   3     624.97613  0.26440  -23.386781   -21.379135   -22.625083 
   4     636.03766  0.13926  -23.185210   -20.559827   -22.189144 
   5     658.36014  0.00016  -23.443271   -20.200150   -22.212836 
   6     669.88472  0.11243  -23.260601   -19.399743   -21.795797 
\end{outbit}
\end{script}

Warning: in finite samples the choice of the maximum lag, $p$, may
affect the outcome of the procedure. \textit{This is not a bug}, but
rather an unavoidable side effect of the way these comparisons should
be made. If your sample contains $T$ observations and you invoke the
lag selection procedure with maximum order $p$, gretl examines
all VARs of order ranging form 1 to $p$, estimated on a uniform sample
of $T-p$ observations. In other words, the comparison procedure does
not use all the available data when estimating VARs of order less than
$p$, so as to ensure that all the models in the comparison are
estimated on the same data range. Choosing a different value of $p$
may therefore alter the results, although this is unlikely to happen
if your sample size is reasonably large.

An example of this unpleasant phenomenon is given in example script
\ref{script:var-lagselect}. As can be seen, according to the
Hannan-Quinn criterion, order 2 seems preferable to order 1 if the
maximum tested order is 4, but the situation is reversed if the
maximum tested order is 6.

\section{Structural VARs}
\label{sec:svar}

Gretl's built-in \cmd{var} command does not support the general class
of models known as ``Structural VARs''---though it does support the
Cholesky decomposition-based approach, the classic and most popular
structural VAR variant. If you wish to go beyond that there is a gretl ``addon''
named \app{SVAR} which will likely meet your needs. \app{SVAR} is
supplied as part of the gretl package, you can find its documentation
(which is quite detailed) as follows: under the \textsf{Tools} menu in
the gretl main window, go to ``Function packages/On local machine.''
(Or use the ``$fx$'' button on the toolbar at the foot of the main
window.) In the function packages window either scroll down or use the
search box to find \textsf{SVAR}. Then right-click and select
``Info.'' This opens a window which gives basic information on the
package, including a link to \texttt{SVAR.pdf}, the full
documentation.

The remainder of this section will thus only deal with the Cholesky-based
recursive shock identification used by the native \cmd{var} command.

\subsection{IRF and FEVD}

Assume that the disturbance in equation \eqref{eq:VAR} can be thought
of as a linear function of a vector of \emph{structural shocks} $u_t$,
which are assumed to have unit variance and to be mutually
unncorrelated, so $V(u_t) = I$. If $\epsilon_t = K u_t$, it follows
that $\Sigma = V(\epsilon_t) = KK'$.

The main object of interest in this setting is the sequence of
matrices
\begin{equation}
  \label{eq:svma}
  C_k = \frac{\partial y_t}{\partial u_{t-i}} = \Theta_k K, 
\end{equation}
known as the structural VMA representation. From the $C_k$ matrices
defined in equation \eqref{eq:svma} two quantities of interest may be
derived: the \textbf{Impulse Response Function} (IRF) and the
\textbf{Forecast Error Variance Decomposition} (FEVD).

The IRF of variable $i$ to shock $j$ is simply the sequence of the
elements in row $i$ and column $j$ of the $C_k$ matrices. In symbols:
\[
  \mathcal{I}_{i,j,k} = \pder{y_{i,t}}{u_{j, t-k}}
\]
As a rule, Impulse Response Functions are plotted as a function of
$k$, and are interpreted as the effect that a shock has on an
observable variable through time. Of course, what we observe are the
estimated IRFs, so it is natural to endow them with confidence
intervals: following common practice, gretl computes the
confidence intervals by using the bootstrap;\footnote{It is possible,
  in principle, to compute analytical confidence intervals via an
  asymptotic approximation, but this is not a very popular choice:
  asymptotic formulae are known to often give a very poor
  approximation of the finite-sample properties.} details are given
later in this section.

Another quantity of interest that may be computed from the structural
VMA representation is the Forecast Error Variance Decomposition
(FEVD). The forecast error variance after $h$ steps is given by
\[
  \Omega_h = \sum_{k=0}^h C_k C_k'
\]
hence the variance for variable $i$ is
\[
  \omega^2_i = \left[ \Omega_h \right]_{i,i} = \sum_{k=0}^h
  \mathrm{diag}(C_k C_k')_i =
  \sum_{k=0}^h \sum_{l=1}^n ({}_kc_{i.l})^2 
\]
where ${}_kc_{i.l}$ is, trivially, the $i,l$ element of $C_k$. As a
consequence, the share of uncertainty on variable $i$ that can be
attributed to the $j$-th shock after $h$ periods equals
\[
  \mathcal{VD}_{i,j,h} =
  \frac{\sum_{k=0}^h ({}_kc_{i.j})^2 }{  \sum_{k=0}^h \sum_{l=1}^n
    ({}_kc_{i.l})^2 } .
\]
This makes it possible to quantify which shocks are most important to
determine a certain variable in the short and/or in the long run.

\subsection{Triangularization}

The formula \ref{eq:svma} takes $K$ as known, while of course it has
to be estimated. The estimation problem has been the subject of an
enormous body of literature we will not even attempt to summarize
here: see for example \cite[chapter 9]{Lutkepohl05}.

Suffice it to say that the most popular choice dates back to
\cite{sims80}, and consists in assuming that $K$ is lower triangular,
so its estimate is simply the Cholesky decomposition of the estimate
of $\Sigma$. The main consequence of this choice is that the ordering
of variables within the vector $y_t$ becomes meaningful: since $K$ is
also the matrix of Impulse Response Functions at lag 0, the
triangularity assumption means that the first variable in the ordering
responds instantaneously only to shock number 1, the second one only
to shocks 1 and 2, and so forth. For this reason, each variable is
thought to ``own'' one shock: variable 1 owns shock number 1,
and so on.

In this sort of exercise, therefore, the ordering of the $y$ variables
is important.
% , and the applied literature has developed the ``most
% exogenous first'' mantra---where, in this setting, ``exogenous'' really
% means ``instantaneously insensitive to structural
% shocks''.\footnote{The word ``exogenous'' has caught on in this
%   context, but it's a rather unfortunate choice: for a start, each
%   shock impacts on every variable after one lag, so nothing is really
%   exogenous here. A better choice of words would probably have
%   been something like ``sturdy'', but it's too late now.} 
To put it
differently, if variable \texttt{foo} comes before variable
\texttt{bar} in the $Y$ list, it follows that the shock owned by
\texttt{foo} affects \texttt{bar} instantaneously, but not
vice versa.

Impulse Response Functions and the FEVD can be printed out via the
command line interface by using the \option{impulse-responses} and
\option{variance-decomp} options, respectively. If you need to store
them into matrices, you could compute the structural VMA and proceed
from there. For example, the following code snippet shows you how to
manually compute a matrix containing the IRFs:
\begin{code}
open denmark
list Y = 1 2 3 4
scalar n = nelem(Y)
var 2 Y --quiet --impulse-responses

matrix K = cholesky($sigma)
matrix V = $vma
matrix IRF = V * (K ** I(n))
print IRF
\end{code}
in which the equality
\[
\mathrm{vec}(C_k) = \mathrm{vec}(\Theta_k K) = (K' \otimes I)
\mathrm{vec} (\Theta_k)
\]
was used.

A more convenient way of obtaining the desired quantities is to use the 
\cmd{irf} and \cmd{fevd} functions which can be used in scripts after
a VAR (or VECM, see the next chapter) has been estimated. In these functions
you must specify the number of the responding (target) variable and the number 
of the analyzed shock to get the corresponding results as a column vector. The 
choice of how many periods should be calculated --and thus how long the result
vector will be-- is determined by previously invoking \cmd{set horizon x}, where
x is a non-negative integer and the first response concerns the impact effect.
As always, it is recommended to consult the function reference under the 
help menu, where in the case of the \cmd{irf} function it is also explained that
the implicit shock size is such that the impact response in the same equation
is one standard deviation (of the corresponding error term).

\subsection{IRF bootstrap}

The IRFs obtained above are estimates and as such they are uncertain. 
Mostly due to the fact that they are nonlinear functions of the VAR parameters
the standard way of assessing this estimation uncertainty and to derive 
confidence intervals or bands is to use a bootstrap approach. Again, more
advanced options are available with the \app{SVAR} addon, but the \cmd{irf}
function used after the built-in \cmd{var} (or \cmd{vecm}) command also 
provides the option to run a bootstrap based on resampling from the residuals.
(The number of bootstrap iterations can be adjusted through \cmd{set boot\_iters x},
where x must be larger than 499.) 
The desired nominal confidence level must be specified after the target and 
shock numbers as the third argument, and in that case the return vector becomes
a three-column matrix where the lower and upper bounds of the confidence intervals
are given in the extra two columns.

\subsection{Menu-driven usage}

Almost all the functionality related to the described (recursively identified) structural 
VARs is also available under the menus in the model window that appears 
after a VAR is estimated in the GUI.\footnote{Note that you cannot directly invoke 
the \app{SVAR} addon from the model window of an estimated VAR; that menu 
entry is only present in gretl's \emph{main} window under the Model menu and 
multivariate time series sub-menu.}

\begin{itemize}
\item In the ``Plots'' menu there are a number of menu entries relating to the
impulse responses as well as one entry for the forecast error variance decomposition.
Selecting any of these will bring up a little specification window where the ordering
for the Cholesky decomposition must be chosen, and in case of IRFs the intended 
bootstrap coverage can be set. 

\item In the ``Analysis'' menu there are also entries for IRF and FEVD, which may 
sometimes be a little confusing. The point is that here the numbers (of the point 
estimates) will be printed out in a tabular format instead of being plotted.

\end{itemize}

\section{Residual-based diagnostic tests}

Three diagnostic tests based on residuals are available after
estimating a VAR---for normality, autocorrelation and ARCH
(Autoregressive Conditional Heteroskedasticity). These are implemented
by the \cmd{modtest} command, using the options \option{normality},
\option{autocorr} and \option{arch}, respectively.

The (multivariate) normality test is that of \cite{doornik-hansen94};
it is based on the skewness and kurtosis of the VAR residuals.

The autocorrelation and ARCH test are also by default multivariate;
they are described in detail by \cite{Lutkepohl05} (see sections 4.4.4
and 16.5.1). Both tests are of the LM type, although the
autocorrelation test statistic is referred to a Rao $F$ distribution
\citep{Rao73}. These tests may involve estimation of a large number of
parameters, depending on the lag horizon chosen, and can fail for lack
of degrees of freedom in small samples. As a fallback, the
\option{univariate} option can be used to specify that the tests
be run per-equation rather than in multivariate mode.

Listing~\ref{script:var-autocorr} illustrates the VAR autocorrelation
tests, replicating an example given by \citet[p.\ 174]{Lutkepohl05}.
Note the difference in the interpretation of the \textsl{order}
argument to \cmd{modtest} with the \option{autocorr} option (this also
applies to the ARCH test): in the multivariate version \textsl{order}
is taken as the \textit{maximum} lag order and tests are run from lag
1 up to the maximum; but in the univariate version a single test is
run for each equation using just the specified lag order. The example
also exposes what exactly is returned by the \dollar{test} and
\dollar{pvalue} accessors in the two variants.

\begin{script}[htbp]
  \scriptcaption{VAR autocorrelation test from L\"utkepohl}
  \label{script:var-autocorr}
Input:
\begin{scodebit}
open wgmacro.gdt --quiet
list Y = investment income consumption
list dlnY = ldiff(Y)
smpl 1960:4 1978:4
var 2 dlnY
modtest 4 --autocorr
eval $test ~ $pvalue
modtest 4 --autocorr --univariate
eval $test ~ $pvalue
\end{scodebit}

Output from tests:
\begin{outbit}
? modtest 4 --autocorr
Test for autocorrelation of order up to 4

          Rao F   Approx dist.  p-value
lag 1     0.615    F(9, 148)     0.7827
lag 2     0.754    F(18, 164)    0.7507
lag 3     1.143    F(27, 161)    0.2982
lag 4     1.254    F(36, 154)    0.1743

? eval $test ~ $pvalue
     0.61524      0.78269 
     0.75397      0.75067 
      1.1429      0.29820 
      1.2544      0.17431 

? modtest 4 --autocorr --univariate
Test for autocorrelation of order 4

Equation 1:
Ljung-Box Q' = 6.11506 with p-value = P(Chi-square(4) > 6.11506) = 0.191

Equation 2:
Ljung-Box Q' = 1.67136 with p-value = P(Chi-square(4) > 1.67136) = 0.796

Equation 3:
Ljung-Box Q' = 1.59931 with p-value = P(Chi-square(4) > 1.59931) = 0.809

? eval $test ~ $pvalue
      6.1151      0.19072 
      1.6714      0.79591 
      1.5993      0.80892 
\end{outbit}
\end{script}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "gretl-guide"
%%% End: 
