Interpreting the Output From R Package vars

There’s a useful package on CRAN for working with VAR models, going by the innovative name vars. This brief post provides an overview of its basic usage and interpretation of the results. Note: The package comes with a comprehensive manual. Consult it to learn about more advanced topics than are covered here. The user interface may leave you confused when you start using it. That’s where this post comes in.

Let’s start by generating a couple of time series variables, y1 and y2.

set.seed(100)
y <- ts(matrix(rnorm(300), ncol=2))
colnames(y) <- c("y1", "y2")

Load the package and estimate a reduced form VAR(2) model:

library(vars)
varfit <- VAR(y, p=2)

You can print out the estimated coefficients and residuals:

varfit
residuals(varfit)

It’s important to note that the residuals do not have time series properties, so you’ll have to add them in as needed:

class(residuals(varfit)) # matrix, not mts!

As long as the variables in y have no omitted values (i.e., they are all complete and cover the same time period), you can add time series properties to the residuals by doing this:

ts(residuals(varfit), end=end(y), frequency=frequency(y))

It’s important to specify the end date, not the start date, since the initial observations drop out prior to estimation when you take lags.

One of the more common operations for a VAR model is to compute impulse response functions. Identification of impulse response functions is a large and growing topic. If you can justify recursiveness, you can apply the irf function. The default is to take a Choleski of the residual covariance matrix.

Let’s take a look at the inner workings of the irf function. Since it’s not an exported function of that package, you can print it out with

vars:::irf.varest

(Note the three : characters. That’s how you get a function from a package if it’s not exported.) Scrolling through the code, we see this:

irs <- .irf(x = x, impulse = impulse, response = response, 
        y.names = y.names, n.ahead = n.ahead, ortho = ortho, 
        cumulative = cumulative)

Okay, a call to .irf. What’s that function look like?

vars:::.irf

That in turn contains this (relevant) code:

if (ortho) {
  irf <- Psi(x, nstep = n.ahead)
}

Looking up that function,

vars:::Psi.varest

we finally have the code we’re after:

params <- ncol(x$datamat[, -c(1:x$K)])
sigma.u <- crossprod(resid(x))/(x$obs - params)
P <- t(chol(sigma.u))

P is the matrix of impact multipliers (the contemporaneous responses of each of the variables to each of the shocks). For our data, we have

params <- ncol(varfit$datamat[, -c(1:varfit$K)])
sigma.u <- crossprod(resid(varfit))/(varfit$obs - params)
P <- t(chol(sigma.u))

We get this for the output:

> P
           y1        y2
y1 0.94005058 0.0000000
y2 0.03015228 0.9483967

That’s only slightly different from the asymptotic version, where you don’t have to mess around with degrees of freedom corrections:

sigma.u <- cov(residuals(varfit))
t(chol(sigma.u))

           y1        y2
y1 0.92717256 0.0000000
y2 0.02973921 0.9354043

The way to interpret this output is that the column name is the name of the shock and the row name is the name of the variable responding to the shock. You can confirm that by plotting the IRFs, which the vars package author has labeled:

plot(irf(varest))

The initial response of y1 to shocks to y2 is equal to zero. The secret to interpreting the output is that you should have an upper triangular matrix after taking the Choleski. The first row is the response of the variable in the first column of your dataset, the second row is the response of the variable in the second column of your dataset, and so on.