The calculations for this are actually quite simple. The h step responses are computed by regressing Y(t) on Y(t-h), Y(t-(h+1)), ... and extracting the coefficient vector on the leading lag (that is, the Y(t-h)) in that. That n x n matrix gives the local projection responses of the n variables (rows) to unit shocks in the n variables (columns). Different shock patterns can be obtained by post-multiplying that matrix by a set of weights on the unit shocks (in this case, by the Cholesky factor of the residuals from a standard VAR). (Standard IRF's effectively do the same thing for multi-step responses except the multi-step projection is a known function of the standard VAR coefficients rather than being estimated directly).

The regression with skipped lags will have serially correlated errors. The standard errors of the IRF's are computed using some form of HAC covariance estimates for that lead coefficient matrix. Because there are different ways of computing those, the standard errors will differ as well.

A few things to note.

- The graphs in the RATS program label the impact shock as 0 (as is more typical) while Jorda uses 1.
- The RATS graphs are organized in the more typical way with variables across and shocks down. They also use a common scale going across for all responses of a given variable.
- This does not do the cubic projections. That's not dramatically more complicated, but it's not clear that there's a great advantage to it, and the cubic responses have to be linearized about some value for the lagged y's anyway (which the linear projection responses are independent of those).