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!