You can use the BayesianVARMA model in pymc_extras.statespace. Example notebooks here and here, and I also talk about them in this video.
You have full control over the coefficient matrix and covariance matrix when doing inference (so you can impose zeros in either), and over the covariance matrix used when doing impulse response function. There’s an option to do cholesky decomposition on the covariance matrix out of the box, but it needs a bit of work. PRs welcome 