What would be the most appropriate method for constructing a bootstrap percentile confidence region?

Question

The situation: let suppose we have a random variable that follows some parametric distribution $X\sim F_{\boldsymbol\theta}$, $\boldsymbol\theta\in\mathbb{R}^2$. We are provided with the bootstrap distribution (parametric, non-parametric, ...) of the estimates $$ \hat{\boldsymbol{\theta}}_{1},\hat{\boldsymbol{\theta}}_{2}, \dots, \hat{\boldsymbol{\theta}}_{B} \sim \hat{f}(\hat{\boldsymbol\theta}), $$ where $B$ represents the number of replicates. We want to measure the uncertainty around $\boldsymbol\theta$ from $\hat{f}(\hat{\boldsymbol\theta})$. Assume that if $\theta_1$ and $\theta_2$ are independent, we can build two confidence intervals from the marginal percentiles of $\hat{f}(\hat{\boldsymbol\theta})$ that leads to the exact overall coverage of $1-\alpha$ (for example by applying Dunn-Sidak correction https://en.wikipedia.org/wiki/%C5%A0id%C3%A1k_correction).

The problem: if $\theta_1$ and $\theta_2$ are ''dependent enough'', the coverage will not be exact and this dependence should be taken into account somehow.

My question: Since we have at our disposal the distribution $\hat{f}(\hat{\boldsymbol\theta})$, it is tempting to use it. How can I construct an equivalent of the percentile method in 2-dimensions, that is, a 2-dimensional percentile confidence region?

Example: We observe a sample of size $n=35$ from a Lomax random variable (https://en.wikipedia.org/wiki/Lomax_distribution). This distribution has two parameters. From my experience, both parameters are highly correlated (see the graph below).

The red points are the estimates $\hat{\boldsymbol{\theta}}_{1},\hat{\boldsymbol{\theta}}_{2}, \dots$ (on the log-scale!)

The black point is the true value $\boldsymbol\theta = (2, 2.3)^T$, on the log-scale.

The black lines are the percentile marginal 95% confidence intervals (no correction, otherwise too wide for the graph).

The green area is the 95% Highest Density Region (see What is a Highest Density Region (HDR)?).

The gray-shaded area was built using the post https://stackoverflow.com/questions/31893559/r-adding-alpha-bags-to-a-2d-or-3d-scatterplot It seems to be ''closer'' to what ''should be'' the confidence region.

Are there other methods? Which one seems more appropriate?

Answer 1

First, I think there may be an error in your calculation. Your green area does not make sense as a highest-density-region (HDR), no matter how it is calculated.

Second, in your case I would propose either a correct HDR. Or alternatively a bagplot, which is a bivariate generalization of the boxplot. As the boxplot contains 50% of the observations in the "box", a bagplot contains 50% of a set of bivariate observations in the "bag", and, importantly, the 50% with the greatest bivariate depth. It also contains "fences" that play the role of the whiskers in the boxplot. The original publication is Rousseeuw, Ruts & Tukey (1999, The American Statistician).

Answer 2

After some thoughts and experiment, I think I have found a valid answer! I simply need to find the solutions in $(x,y)^T$ such that $F_n(x,y)$ equals $0.025$ or $0.975$, where $F_n(\cdot)$ is the bivariate bootstrap empirical cumulative distribution. In the example I find the quantiles numerically. See the blue lines on the new graph. I think I have never seen that before, certainly because quantiles are not easy to obtain (or at least from my short experience) in dimension larger than one. I was expecting HDR method (see the question) to cover the area between the two blue lines. So maybe HDR is not appropriate in the situation reported (look at the graph in the question).