At the start, you have a fundamental decision to make:
Are you trying to learn something from the data? Or are you
trying to teach the data to behave as you suppose they should?
This answer is mainly oriented toward the former approach.
Usually, it is a mistake to remove an 'outlier' from a dataset
unless you can establish that the observation in question arose from a documentable error (equipment failure, data entry error, etc.)
or you know for sure that its value is impossible (person's age above 140, negative height, etc.)
Anecdote: Where I live, the highest heating bill arises in December and January. For personal and family reasons unlikely ever to be
repeated, my energy bills for 12/2019 and 1/2020 were extraordinarily high. By your criterion, I might be removed from your list, which I think would be a mistake. I can't foresee the same circumstances
ever coming around for me again, but others may well be surprised by similar temporary periods of high energy usage in the future for very similar reasons.
One way to stabilize averages without tampering with data is to use trimmed means. To find a trimmed mean, data are sorted, a certain percentage of the very lowest and of the very highest
observations are ignored, and the mean of the more central remaining observations is taken. Depending on circumstances, typical trimming percentages may be 2% to 20% (sometimes higher), leaving the central 96% to 60% (sometimes fewer) to be averaged.
Consider data with $n=1000$ observations from
a gamma distribution with shape parameter 10 (perhaps waiting
times for finishing multi-phase projects.) Here is a graph of
its density curve--made in R.
curve(dgamma(x,10,1), 0, 25, lwd=2, ylab="PDF",
main="Density of GAMMA(10, 1)")
abline(v=0, col="green2"); abline(h=0, col="green2")
Almost all such samples have at least one boxplot outlier and the average number of outliers in a sample of 1000 is about 14.
nr.out = replicate(10^5,
Let's take a look at boxplots of 20 samples of size 1000 from this distribution in order to see the outliers.
m = 20; n=1000
x = rgamma(m*n,10,1); g = rep(1:m, n)
main="GAMMA(10,1) Population: Boxplots of 20 Samples of 1000")
It seems that 2% trimmed means of the 1000 observations in each sample should allow us to ignore boxplot outliers in finding the means. (But the trimmed values are are not removed, so the ordinary mean
and quartiles are not affected.)
My simulated gamma observations have $\mu = 10, \sigma^2 = 10,$ so samples of
1000 have ordinary means average about $10$ with variances about $0.01$ (from theory). By contrast
2% trimmed means of samples average about $9.93$ with variances about $0.01$ (from simulation).
a.02 = replicate(10^5, mean(rgamma(1000,10,1),trim=.02))
By using trimmed means we have retained all of the
data. In a fair and systematic way, we have mainly avoided using boxplot outliers to estimate means. Very roughly speaking we have
computed trimmed means by ignoring values that are more than double the
ordinary mean. And at the same time we have ignored values that are less than half the ordinary mean. Perhaps we find that
temporarily ignoring values that are proportionately far from
the ordinary mean (still the best estimate of the population mean), we can make better judgments from our data.
qgamma(c(.02,.98), 10, 1)
 4.618349 17.509813
However, over time we may come to realize that all of the observations have
a legitimate role to play in understanding how to use the data
to best advantage. In that case, the data are intact and we can do so.
Note: There are distributions with such heavy tails that
a trimmed sample mean is a better estimate of the location of the population than an ordinary sample mean.
The Cauchy is one such distribution. In that case tails are so heavy that
a 38% trimmed mean seems optimal. See a brief discussion
and further information at its links.