Program Impulse Response Functions for VAR

Question

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!

Answer 1

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]

Answer 2

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.