I had posted that older question. I only noticed the issue because so many papers use an analytic Bayesian approach. The analytic formula was often more likely to be stationary than the simulated set of parameters I was dealing with.
In retrospect, it's a hard problem that I shouldn't have expected an answer from on the site. Gary Koop has some papers that you might find informative and can check out the references therein.
I don't much use Gibbs sampling these days. Regardless, here are probably the things I would focus on
Fit the Bayesian VAR with stationary data. This means that there
won't be any coefficients near one and it is a little easier to
think about setting priors tightly near zero.
A corollary of 1 is to use a VECM structure when dealing with cointegration, rather than a regression in levels. Again, keeping the coefficients near zero is helpful for thinking about priors.
Reduce the dimension. I don't mean throw out
some variables. But if you have a large VAR, you can instead convert it to an factor-augmented VAR with PCA and only do the VAR on the smaller subset. This will reduce the computational load.
Use something like a Minnesota prior on the
I'm using Stan more recently, though I haven't fit VARs in it. I would probably use some of the above tricks. Stan would probably do a good job with considering the correlation among the parameters. Really, Stan was created to deal with hierarchical parameters, where the problem is correlations between them. For a given parameter, I would expect the coefficients on its own lags to have stronger correlation than the coefficient of the other variables. Nevertheless, anytime I think of something like rejecting non-stationary sets of parameters I just imagine that it would be incredibly slow.