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:

x <- rnorm(100)
y <- rnorm(100)
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!


2 Answers 2


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")

Notice, there is also a difference between

summary(model)$covres and

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

  • $\begingroup$ 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. $\endgroup$
    – fredrikhs
    Apr 25, 2013 at 9:35
  • $\begingroup$ 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. $\endgroup$
    – hejseb
    Apr 25, 2013 at 9:40
  • $\begingroup$ 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. $\endgroup$
    – fredrikhs
    Apr 25, 2013 at 9:47
  • $\begingroup$ So how would I obtain the structural shock? I'm still a little confused, unfortunately. $\endgroup$
    – hejseb
    Apr 25, 2013 at 11:14
  • $\begingroup$ 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) $\endgroup$
    – fredrikhs
    Apr 25, 2013 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.