# How to calculate the impulse response function of a VAR(1)? (With example)

How to calculate:

1) Simple IRF

2) Orthological IRF (Y2 -> Y1)

Of an unrestricted VAR(1) model:

$Y_{1, t} = A_{11}Y_{1, t-1} + A_{12} Y_{2, t-1} + e_{1,t}$ ,

$Y_{2, t} = A_{21}Y_{1, t-1} + A_{22} Y_{2, t-1}+e_{2,t}$

Let's just say that $A_{11} = 0.8$, $A_{12} = 0.4$, $A_{21} = -0.3$, $A_{22} = 1.2$

And the shock size is 1 to both residuals. You don't have to use the provided values as long as the point gets across. Let's also say that the IRF length is 4. I think this should be enough info but let me know if something else is needed.

Bonus question: How does the response change in a structural VAR (any structure)?

• @RichardHardy This question was motivated by the lack of detail to the process in the manuals of statistical packages or any internet source. As such I don't think it classifies for self-study tag. – Dole Mar 24 '15 at 20:42

For a VAR(1), we write the model as $$y_t=\Pi y_{t-1}+\epsilon_t$$ where $y$ and $\epsilon$ are $p\times 1$ vectors. If you have more lags, the idea of extension is the same (and it is particularly straight-forward using the companion form).

The impulse response is the derivative with respect to the shocks. So the impulse response at horizon $h$ of the variables to an exogenous shock to variable $j$ is $$\frac{\partial y_{t+h}}{\partial \epsilon_{j, t}}=\frac{\partial }{\partial \epsilon_{j, t}}\left(\Pi y_{t+h-1}+\epsilon_{t+h-1}\right)=\cdots=\frac{\partial }{\partial \epsilon_{j, t}}\left(\Pi^{h+1} y_{t}+\sum_{i=0}^h\Pi^i\epsilon_{t+h-i}\right).$$ This derivative will eliminate all terms but one, namely the term in the sum which is $\Pi^h\epsilon_t$, for which we get $$\frac{\partial y_{t+h}}{\partial \epsilon_{j, t}}=\frac{\partial }{\partial \epsilon_{j, t}}\left(\Pi^{h+1} y_{t}+\sum_{i=0}^h\Pi^i\epsilon_{t+h-i}\right)=\frac{\partial }{\partial \epsilon_{j, t}}\Pi^h\epsilon_{t}=\Pi^he_j$$ where $e_j$ is the $j$th row of the $p\times p$ identity matrix. That is, the response of all $p$ variables at horizon $h$ to a shock to variable $j$ is the $j$th column of $\Pi^h$. If you take the derivative with respect to the matrix $\epsilon_t$ instead, the result will be a matrix which is just $\Pi^h$, since the selection vectors all taken together will give you the identity matrix.

That is the non-orthogonalized case without identification, which I believe is not so common in the literature. What people usually use is either some sophisticated identification scheme, or more often a Cholesky decomposition. To study this, it is more convenient to work with the vector moving average form of the model (which exists if it is stationary) $$y_t=\sum_{s=0}^\infty\Psi_s\epsilon_{t-s}.$$ The problem for interpretation is when the error terms are correlated, because then an exogenous shock to variable $j$ is simultaneously correlated with a shock to variable $k$, for example. To eliminate this, you can use a Cholesky decomposition which orthogonalizes the innovations. Let's suppose that the covariance matrix of the errors is $\Omega$. We decompose it as $\Omega=PP'$ and introduce $v_t=P^{-1}\epsilon_t$ which are error terms with the identity matrix as covariance matrix. Do some manipulation: $$y_t=\sum_{s=0}^\infty\Psi_s\epsilon_{t-s}=\sum_{s=0}^\infty\Psi_sPP^{-1}\epsilon_{t-s}=\sum_{s=0}^\infty\Psi_s^*v_{t-s}.$$ where $\Psi_s^*=\Psi_sP$. Consider now the response to an orthogonalized shock: $$\frac{\partial y_{t+h}}{\partial v_{j, t}}=\frac{\partial }{\partial v_{j, t}}\left(\sum_{s=0}^\infty\Psi_s^*v_{t+h-s}\right)=\Psi_h^*e_j.$$

To calculate this in practice, you will need to find the moving average matrices $\Psi$. This you do recursively. If you have $K$ lags: $$\Psi_s=0, \quad (s=-K+1, -K+2, \dots, -1)\\ \Psi_0=I\\ \Psi_s=\sum_{i=1}^K\Pi_i\Psi_{s-i}, \quad (s=1, 2, \dots).$$ With estimates, you just put hats on the $\Pi$ matrices and proceed. $P$ we find from using a Cholesky decomposition of the estimated error covariance matrix, $\hat\Omega$.

Edit: In univariate time series analysis, one standard result is that every AR process can be written as an MA($\infty$) process. You have the same result for multivariate time series, meaning that we can always rewrite a stationary VAR($p$) as a VMA($\infty$). This is central to impulse response analysis. The case with only one lag is the easiest. In this case, we may write $$y_t=\Pi y_{t-1}+\epsilon_t=\Pi(\Pi y_{t-2}+\epsilon_{t-1})+\epsilon_t=\cdots=\sum_{s=0}^\infty \Pi^i\epsilon_{t-s}.$$ The implied steps in the $\cdots$ part might not be obvious, but there is just a repeated substitution going on using the recursive nature of the model. So for the VAR(1), the moving average coefficients $\Psi_s$ are just $\Psi_s=\Pi^s$. For more lags, it gets a little more complicated, but above you will find the recursive relations.

In impulse response analysis, the moving average form of the model is particularly convenient. The reason is that if you want to find the response of $y_{t+h}$ to a shock to $\epsilon_{j, t}$, then if you start with the usual VAR(1) form $$y_{t+h}=\Pi y_{t+h-1}+\epsilon_{t+h},$$ then there is no $\epsilon_t$ in your model as it stands, but you will have to do recursive substitution until you get to it (as I did in the beginning). But, if you have the moving average form of the model, you have it immediately on the right hand side. So for the VAR(1), you will find that $$\frac{\partial y_{t+h}}{\partial \epsilon_{j, t}}=\frac{\partial}{\partial \epsilon_{j, t}}\left(\sum_{s=0}^\infty\Psi_s\epsilon_{t+h-s}\right)=\Psi_he_j=\Pi^he_j,$$ where $e_j$ again is the $j$th column of the $p\times p$ identity matrix. As you see, this is the same result as we found in the beginning, but here we used the moving average form of the model to do it. But the two representations are just two sides of the same coin.

• I really dropped out at the part where the equation was converted to moving average form. Is the orthogonal IRF equation (using cholesky decomposition) even the same equation being estimated via OLS anymore? (IE does the VAR equation and thus coefficients actually change?) In addition, is the error matrix purposely written as $e$ in the first equation or is it supposed to be $e_t$? – Dole Mar 25 '15 at 21:28
• @Dole The IRFs are not estimated per se, they are functions of the parameter matrices, which in turn are estimated. If you don't do orthogonalization, you can still compute them using the moving average way (but you use $P=I$ in the equations above). I'll edit my post to make it clearer. And yes, that is well spotted, that should be $\epsilon_t$. – hejseb Mar 26 '15 at 6:17
• Thanks, I definitely understand the point of the moving average transformation now. However, I always thought that using the Cholesky decomposition for an orthogonalized IRF adds a [1, 0, // B, 1) matrix to the left side of the equation (// marking a change of column). Am I conflating the concept of orthogonal IRF with some other concept here? – Dole Mar 27 '15 at 16:16
• @Dole Yes, I think you might be confusing it with something else. I'm not sure what, though. I guess that you could just as well work with the transformed model which you'd obtain by premultiplying by $P$, i.e. $P y_t=P\Pi y_{t-1}+P\epsilon_t$ since that would have orthogonal errors, but I'm not sure that is what you're thinking about. – hejseb Mar 30 '15 at 5:20
• The following VAR presentation has the equation in the form I spoke about earlier, slightly past the 3 minute mark: "youtube.com/watch?v=DUuIqXbaUAc". I thought that using OLS/GMM on that equation was the way to obtain orthogonal IRFs (the leftmost matrix being the Cholesky decomposition), but it appears I was wrong. Do you know if the resulting orthogonalized(described by you) IRFs and the structural var(the equation in the video) IRFs will still be the same though, given that the idea behind them seems basically the same? – Dole Apr 4 '15 at 12:17

Here is how to calculate the simple IRF.

The idea is to compare a base case where the innovations are

$$(\varepsilon_{1,t+1},\varepsilon_{1,t+2},...)=(0,0,...)$$ $$(\varepsilon_{2,t+1},\varepsilon_{2,t+2},...)=(0,0,...)$$

to an alternative case where the innovations are

$$(\varepsilon_{1,t+1},\varepsilon_{1,t+2},...)=(1,0,...)$$ $$(\varepsilon_{2,t+1},\varepsilon_{2,t+2},...)=(0,0,...)$$

for example (corresponding to a one-time shock of size 1 to $y_1$).

In a VAR(1) system, the $y_1$'s corresponding to the base case will be

$y_{1,t+1} = a_{11} y_{1,t} + a_{12} y_{2,t} + 0$
$y_{1,t+2} = a_{11} y_{1,t+1} + a_{12} y_{2,t+1} + 0 = a_{11} (a_{11} y_{1,t} + a_{12} y_{2,t} + 0) + a_{12} (a_{21} y_{1,t} + a_{22} y_{2,t} + 0) + 0$
$y_{1,t+3} = ...$

The $y_1$'s corresponding to the alternative case will be

$y_{1,t+1} = a_{11} y_{1,t} + a_{12} y_{2,t} + 1$
$y_{1,t+2} = a_{11} y_{1,t+1} + a_{12} y_{2,t+1} + 0 = a_{11} (a_{11} y_{1,t} + a_{12} y_{2,t} + 1) + a_{12} (a_{21} y_{1,t} + a_{22} y_{2,t} + 0) + 0$
$y_{1,t+3} = ...$

The impulse-responses for $y_1$ will be the difference between the alternative case and the base case, that is,

$ir_{1,t+1} = 1$
$ir_{1,t+2} = a_{11}$
$ir_{1,t+3} = ...$

Analogously, you could obtain the impulse responses of a one-time shock of size 1 to $y_1$ on $y_2$. They would be

$ir_{2,t+1} = 0$
$ir_{2,t+2} = a_{21}$
$ir_{2,t+3} = ...$

There must be a more compact way of writing it out, but I wanted to be clear and show it step by step.

Extending this to different kinds of shocks (e.g. unit shock to both $y_1$ and $y_2$ at time $t+1$ followed by zero shocks afterwards) should be straightforward.

Note: it might be more common to consider a shock at time $t$ rather than $t+1$, but that does not change the essence.

• Thanks, perfect answer for the simple IRF case! For some reason eviews prints out IRFs with just slightly different values to what I get calculating by hand. Must be an interpolation issue or something. – Dole Mar 24 '15 at 22:20
• @Dole IIRC, the default option in EViews is to use a Cholesky decomposition. Are you sure you're comparing the same numbers (i.e. non-orthogonalized)? – hejseb Mar 25 '15 at 7:10
• @hejseb That's correct, I did change the IRF to simple one unit shock. To be clear I did not export the values but rather looked at the IRF graphs where eviews prints the "precise" values if the navigator is hovered over the graph long enough. – Dole Mar 25 '15 at 21:17