Coda, and other R packages, allow to calculate HPD intervals for multiple posterior chains. However, the HPD intervals are always calculated marginally per chain.
Is there a way in R to find the joint HPD interval, taking into account the correlation structure of the chains?