# Parametric bootstrap samples include impossible values - should they be removed?

Suppose I am fitting a generalised pareto distribution (GPD) to data, and wish to use a parametric bootstrap to estimate uncertainty on the fitted parameters and quantiles.

The GPD depends on 3 parameters. It always has a lower bound, and depending on the parameter values it may also have an upper bound.

If I have a dataset sampled from a GPD distribution with unknown parameters, I know for sure that the upper bound of the parent distribution must be greater than (or equal) the maximum value in the dataset.

My problem is that this fact is not respected with a naive parametric bootstrap of the GPD parameters. See code below which illustrates this (in that example, 10% of the samples correspond to impossible values, and confidence intervals are likewise affected).

The fitting is done with L-moments but I think the same problem should apply with other fitting methods (e.g. standard maximum likelihood) if a parametric bootstrap is used.

Intuitively, the problem seems to be that when we draw a random sample from the distribution which is fit to the data (test_data_fit in the code below), the resulting sample does not contain the 'extra' information about the upper bound in the original data.

Question: Is there a valid way to 'remove' such unrealistic bootstrap parameter estimates?? For example, can I just delete them, and compute confidence intervals from the bootstrap samples which have realistic parameters? Or is there something else that can be done?

I know there are non-bootstrap approaches to deriving confidence intervals that will not have this problem (e.g. maximum likelihood + profile likelihood intervals). However, I have an application where a GPD fit like this is part of a much more complex model, for which a bootstrap type approach is much more straightforward.

# Fit with L-moments
library(lmomco)

#############################################
#
# MAKE TEST DATA
#

# Construct the 'known' distribution parameters
gpd_par = c(1.0, 0.56, 0.206)
names(gpd_par) = c('xi', 'alpha', 'kappa')
true_distribution = list(
type = 'gpa',
para = gpd_par,
zeta = 1,
source = 'pargpa')

# Make a dataset to work with
set.seed(123456)
test_data = rlmomco(200, true_distribution)

# In practical situations we don't know the true_distribution, BUT from the
# data, we still know that the TRUE distribution must have an upper bound >=
# data_maximum
data_maximum = max(test_data)
data_maximum
# [1] 2.765597

###########################################
#
# PARAMETRIC BOOTSTRAP
#

# Step1: Fit the GPD distribution to the data ('gpa' in lmomco)
test_data_fit = lmom2par(lmom.ub(test_data), 'gpa')

# Step2: Compute gpd fits from random samples drawn from 'test_data_fit'
test_data_bootstrap_par = replicate(1000,
lmom2par(lmom.ub(rlmomco(length(test_data), test_data_fit)), 'gpa'),
simplify=FALSE
)

# The problem: I know some of the bootstrap parameter combinations are
# impossible, because they imply distributions with an upper bound  < data_maximum
bootstrap_upper_bound = unlist( lapply(test_data_bootstrap_par, f<-function(x) par2qua(1, x)))

# We get a significant fraction of results which do not conform to the limit
fraction_impossible_values = sum(bootstrap_upper_bound < data_maximum) / length(bootstrap_upper_bound)
fraction_impossible_values
# [1] 0.102

# Clearly then confidence intervals will be affected, for example a 95% percentile interval includes impossible values.
quantile(bootstrap_upper_bound, p=c(0.025, 0.975))
#    2.5%    97.5%
# 2.579955 5.160431

# Similarly, estimates of upper quantiles will presumably be affected.