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!
 A: 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]

A: 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. 
