When there is ARCH effects on VAR residuals $\varepsilon_t$, we can use a GARCH model to remove them : $\zeta_t = \Sigma_{t|t-1}^{-\frac{1}{2}} \varepsilon_t$. Following [Lutkepohl, New Introduction to Multiple Time Series Analysis], coefficients of the IRF correspond to the coefficients of the well-known Moving-Average representation :
$\Delta P_{t} = J \mathbf{\Delta P}_{t} = \sum_{i=0}^{\infty} J \mathbf{\Phi}^i J^{\tau} J \mathbf{E}_{t-i} = \sum_{i=0}^{\infty} J \mathbf{\Phi}^i J^{\tau} \Sigma_{t-i|t-i-1}^{\frac{1}{2}} \Sigma_{t-i|t-i-1}^{-\frac{1}{2}} \varepsilon_{t-i} = \sum_{i=0}^{\infty} \mathbf{\Psi}_{i,t} \zeta_{t-i}$
And so the impulse response for a given $t$ after $N$ periods is the corresponding MA coefficient : $\mathbf{\Psi}_{N,t} = J \mathbf{\Phi}^N J^{\tau} \Sigma_{t-N|t-N-1}^{\frac{1}{2}}$.
I am not sure about the indices of $\Sigma$, should I use forecast of it ? Thanks for your help.
EDIT : Maybe I could directly apply IRF on VAR residuals, is it a necessary condition to remove ARCH effects to use IRF ? I don't think so, and so I could use : $\Delta P_{t} = \sum_{i=0}^{\infty} J \mathbf{\Phi}^i J^{\tau} \varepsilon_{t-i} = \sum_{i=0}^{\infty} \mathbf{\Psi}_{i,t} \varepsilon_{t-i}$.
And so the impulse response after N periods is the corresponding MA coefficient : $\mathbf{\Psi}_{N} = J \mathbf{\Phi}^N J^{\tau}$, which not depends on t anymore.