The following is a simple proof of the (perhaps obvious) fact that a Bayesian credible interval (or set) always has average frequentist coverage equal to the nominal level (assuming a proper prior). Some authors (e.g. Brown et al. 2001 and Fagerland, Lydersen & Laake 2017, section 1.4.2), in choosing between the large number of approximate frequentist confidence intervals, argue that such average coverage may be more relevant than ensuring that the coverage at least has the nominal level for all $\theta$. If one accepts this, an average weighted by prior belief in different values of $\theta$ is arguably a sensible way to assess overall coverage (as done below).
By definition, a Bayesian $(1-\alpha)$-credible interval (or set) $C(x)$ (some function of the data $x$) includes $\theta$ with probability
$$
P(\theta \in C(x)|x)=1-\alpha. \tag{1}
$$
conditional on the data $x$. Thinking of both $\theta$ and $x$ as random variables, it is also true, by the law of total probablity, that
$$
P(\theta \in C(x))=\int P(\theta \in C(x)|x)f(x)dx=1-\alpha \tag{2}
$$
where $f(x)$ is the marginal density of $x$.
For the frequentist coverage of the Bayesian credible interval, we instead condition on only the parameter $\theta$ such that $x$ and $C(x)$ are random. This gives the frequentist coverage
$$
P(\theta \in C(x)|\theta) \tag{3}
$$
which in general will depend on $\theta$. But the average frequentist coverage, weighted by the prior $\pi(\theta)$, is
$$
\int P(\theta \in C(x)|\theta) \pi(\theta) d\theta=P(\theta \in C(x)) =1-\alpha, \tag{4}
$$
by the law of total probability and (2), and so always perfectly matches the nominal level.
The following plot and R-code illustrates (3) and (4) for a binomial model with a uniform prior on $p$ and $n=20$:
R code
m <- 1e+5
p <- seq(0,1,len=m)
coverage <- numeric(m)
alpha <- .05
a <- 1
b <- 1
n <- 20
x <- 0:n # possible outcomes
lower <- qbeta(alpha/2, a+x, b+n-x) # lower and upper credible limits
upper <- qbeta(1-alpha/2, a+x, b+n-x) # for each value of x
for (i in 1:length(p)) {
px <- dbinom(x,size=n,prob=p[i]) # probability of possible outcomes given p
coverage[i] <- sum(px*ifelse(lower<=p[i] & p[i]<=upper, 1,0)) # coverage
}
plot(p,coverage,type="l",ylim=c(.85,1))
abline(h=1-alpha, lty=3)
> mean(coverage)
[1] 0.9499928
