# Program Impulse Response Functions for VAR

I'm trying to program impulse response functions for a VAR model using Cholesky decomposition. The thing is I do not completely understand how I should do this when I read in the literature. Suppose I have:

$$\begin{bmatrix}x_t\\y_t\\z_t\end{bmatrix}=\begin{bmatrix}\alpha_1\\\alpha_2\\\alpha_3\end{bmatrix}+\begin{bmatrix}b_{11}&b_{12}&b_{13}\\b_{21}&b_{22}&b_{23}\\b_{31}&b_{32}&b_{33}\end{bmatrix}\begin{bmatrix}x_{t-1}\\y_{t-1}\\z_{t-1}\end{bmatrix}+\begin{bmatrix}u_{1t}\\ u_{2t}\\ u_{3t}\end{bmatrix}$$ which we can write as $$\mathbf{x}_t=\mathbf{a}+\mathbf{B}\mathbf{x}_{t-1}+\mathbf{u}_t.$$

Further, suppose the covariance matrix of $\mathbf{u}_t$ is $$Cov(\mathbf{u}_t)=\Sigma_u=PP^\prime.$$

Now, let's say I want the impulse responses to a unit shock in $u_{1t}$. I want the effects then on say $\mathbf{x}_t, \, \mathbf{x}_{t+1}, \, \mathbf{x}_{t+2}$ and $\mathbf{x}_{t+3}$. As I understand it, the orthogonalization is done by multiplying the error vector with $P$. Let's call the responses in period $t$ to the shock $\mathbf{x}_t^*$. Would what I am interested in then be (assume unit shock in $u_{1t}$ such that $\mathbf{u}_t^*=\begin{bmatrix}1&0&0\end{bmatrix}^\prime$):

$$\mathbf{x}_t^*=P\mathbf{u}^*_t=P\begin{bmatrix}1\\0\\0\end{bmatrix}\\ \mathbf{x}_{t+1}^*=\mathbf{B}\mathbf{x}^*_t=\mathbf{B}P\mathbf{u}^*_t\\ \mathbf{x}_{t+2}^*=\mathbf{B}\mathbf{x}^*_{t+1}=\mathbf{B}\mathbf{B}P\mathbf{u}^*_t\\ \mathbf{x}_{t+3}^*=\mathbf{B}\mathbf{x}^*_{t+2}=\mathbf{B}\mathbf{B}\mathbf{B}P\mathbf{u}^*_t$$

Now, extend this to include more lags (for example 4). The model is then $$\mathbf{x}_t=\mathbf{a}+\sum_{k=1}^4\mathbf{B}_k\mathbf{x}_{t-k}+\mathbf{u}_t.$$

Thus the impulses are:

$$\mathbf{x}_t^*=P\mathbf{u}^*_t=P\begin{bmatrix}1\\0\\0\end{bmatrix}\\ \mathbf{x}_{t+1}^*=\mathbf{B}_1\mathbf{x}^*_t=\mathbf{B}_1P\mathbf{u}^*_t\\ \mathbf{x}_{t+2}^*=\mathbf{B}_1\mathbf{x}^*_{t+1}+\mathbf{B}_2\mathbf{x}^*_t=\mathbf{B}_1\mathbf{B}_1P\mathbf{u}^*_t + \mathbf{B}_2P\mathbf{u}^*_t\\ \mathbf{x}_{t+3}^*=\mathbf{B}_1\mathbf{x}^*_{t+2}+\mathbf{B}_2\mathbf{x}^*_{t+1}+\mathbf{B}_3\mathbf{x}^*_{t}=\mathbf{B}_1\mathbf{B}_1\mathbf{B}_1P\mathbf{u}^*_t + \mathbf{B}_1\mathbf{B}_2P\mathbf{u}^*_t+\mathbf{B}_2\mathbf{B}_1P\mathbf{u}^*_t+\mathbf{B}_3P\mathbf{u}^*_t$$

Is this line of thinking correct? If so, then this simple R code should be fine:

library(vars)
set.seed(1)
x <- rnorm(100)
set.seed(2)
y <- rnorm(100)
set.seed(3)
z <- rnorm(100)
data <- cbind(x, y, z)

model <- VAR(data, p=4, type = "const")
u <- matrix(c(1, 0, 0), ncol=1)
P <- chol(cov(residuals(model)))
B1 <- Acoef(model)[[1]]
B2 <- Acoef(model)[[2]]
B3 <- Acoef(model)[[3]]
B4 <- Acoef(model)[[4]]

xt <- P %*% u
xt1 <- B1 %*% xt
xt2 <- B1 %*% xt1 + B2 %*% xt
xt3 <- B1 %*% xt2 + B2 %*% xt1 + B1 %*% xt


Any input would be very much appreciated!

Maybe we can help eachother, I am working on a similar problem. It is my understanding that the impulse response is the response of one variable to a structural shock in another, the structural shock given by

$x_t=P*u_t$ . Why do you write that $$u_t=\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}$$, and not the actual residuals?

Have you compared you findings to the irf() and Psi() function of library(vars)? the code thought was correct for finding the $x_t$'s are:

model <- VAR(data, p=4, type = "const")
u<-residuals(model)
P=t(chol(summary(model)$covres)) x<-P%*%t(u)  Notice, there is also a difference between summary(model)$covres and
cov(residuals(model))


I don't know where this difference comes from to be honest, but if you compare you results with the Psi(model)-function which is used by irf() you find that the first covariance matrix is the one used to obtain the

Psi(model)[,,1]

• the VAR model creates a set of residuals for each regression, one for each point in time t left, after you've used the 4 lags in the VAR(4). t=1,2,3,4,n-p residuals. residuals(model) will give the $u_t$'s. Apr 25 '13 at 9:35
• Yeah, and I tried it now but x<-P%*%t(u) yields 3 by 396 matrix. And then what? Super thanks for your input though, feels like we're close. Apr 25 '13 at 9:40
• x is the structural errors at each point in time, it's kind of awkward that we called it x, when x is also one of your variables thought. But x is the structural error, and the impulse response is $\frac {\delta y_{i,t+s}} { \delta x_{j,t}} = P_{i,j}^s$ where $y$ is the dependent variable and $x$ is the structural shock. Apr 25 '13 at 9:47
• So how would I obtain the structural shock? I'm still a little confused, unfortunately. Apr 25 '13 at 11:14
• P is the coefficient matrix, it is the change in $y_{j,t}$ from a one unit change in $x_{i,t}$. The after what I've read on the subject, the structural shocks are obtained by multiplying $P$ and $u$, like I wrote earlier; x<-P%*%t(u) Apr 25 '13 at 12:02

I just came across this question. Have you tried the irf() function in the vars package for R? I'm pretty sure that's what you're asking for.