Doubt about level of confidence I run 100 simulations with CI at 99% and different seeds. For each simulation I created a report that contains different performance indexes (means) with their CI. However, I noticed that one index over 5 simulations didn't contain the expected mean. Maybe I'm wrong but when I choose 99% confidence level does it mean that 1 simulation over 100 doesn't contain expected mean?
Edit
I can't upload all source code because is too long in fact I realized a discrete-event simulator. So I can upload a simple code to show my doubt. In essence for every simulation (100 altogether) this code samples an exponential distribution (or another it is not important) with mean 2 until precision (half-width CI/average) is 3% or below. For this part of algorithm I adapted code from chapter 8 of Leemis-Park. Then it increments an accumulator if CI is out from expected mean.
import umontreal.ssj.probdist.StudentDist;
import umontreal.ssj.randvar.ExponentialGen;
import umontreal.ssj.randvar.RandomVariateGen;
import umontreal.ssj.rng.MRG32k3a;

import java.util.concurrent.ThreadLocalRandom;

public class Estimate {
    public static void main(String[] args) {
        long outlier = 0;
        for (int simulation = 0; simulation < 100; simulation++) {
            int numObs = 0;
            double level = .99;
            double currentAvg = 0, currentSumVariance = 0;
            long[] seed = new long[6];
            for (int i = 0; i < 3; i++) seed[i] = ThreadLocalRandom.current().nextLong(4294967087L);
            for (int i = 3; i < 6; i++) seed[i] = ThreadLocalRandom.current().nextLong(4294944443L);
            MRG32k3a.setPackageSeed(seed);
            RandomVariateGen randVar = new ExponentialGen(new MRG32k3a(), .5);

            do {
                double diff = randVar.nextDouble() - currentAvg;
                currentSumVariance += diff * diff * numObs / (double) ++numObs;
                currentAvg += diff / numObs;
            } while (numObs < 2 || StudentDist.inverseF(numObs - 1, .5 * (level + 1)) * Math.sqrt(currentSumVariance / (numObs - 1)) > (.03 * currentAvg) * Math.sqrt(numObs));
            double t = StudentDist.inverseF(numObs - 1, .5 * (level + 1));
            double ci = t * Math.sqrt((currentSumVariance / (numObs - 1)) / numObs);

            System.out.println(currentAvg + " ± " + ci);
            System.out.println("Precision (accuracy) = " + (ci / currentAvg * 100) + "%");
            System.out.println("numObs = " + numObs);
            System.out.println("expected mean = " + randVar.getDistribution().getMean());
            System.out.println("===========");
            if (!(currentAvg - ci <= randVar.getDistribution().getMean() && randVar.getDistribution().getMean() <= currentAvg + ci))
                outlier++;

        }
        System.out.println("outlier = " + outlier);
    }
}

Another thing I used a library for RNG and random variate SSJ by Pierre L'Ecuyer.
 A: First, I want to reinforce two caveats from comments: (a) Per @Dave: $n = 100$ may be too small to see a
clear pattern. If you have one defective part among 99 good parts in a sample of size 100, then the number $X$ of bad parts you get in 100 individual draws with replacement would be $X \sim\mathsf{Binom}(100, .01).$ Thus the number of bad parts in one experiment with 100 draws might
be that I never see the bad part.
set.seed(1234)
x = rbinom(1, 100, .01)
x
[1] 0

But in 100,000 such draws of 100 parts,
the average number $Y$ of bad parts seen in $100\,000$ draws is $1.$  Then I might see
a number near to $1.$
set.seed(1235)
y = mean(rbinom(10^5, 100, .01)); y
[1] 1.00104 

Second, per @MichaelLew: your goal is to get a confidence interval
for the population mean $\mu,$ not for its
estimate $\bar X.$
Finally, to look at a particularly simple case, suppose you are taking samples of size $n = 100$ from the normal population $\mathsf{Norm}(\mu, \sigma=1).$ Then, the z confidence interval for $\mu$ based on $n = 100, \sigma=10$ is of the form
$\bar X \pm 1.96\sigma/\sqrt{n}$ or
$\bar X \pm 1.96.$
Then you would hope to find that the CI covers (includes) the population mean $\mu = 85$ in about 95% of the $100\,000$ experiments, as illustrated by the simulation (using R) below.
set.seed(2022)
mu = 85
a = replicate(10^5, mean(rnorm(100, mu, 10)))
summary(a)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  80.29   84.33   85.00   85.00   85.67   89.46 
 mean((a - 1.96 < 85) & (a + 1.96 > 85))
 [1] 0.95092   # aprx. coverage prob of 95% z CI

A: Caveats about the correct terminology aside, you have implemented Algorithm 8.2.3, pp. 364 in "Discrete-Event Simulation" incorrectly. The relevant line is this while predicate:
while ((n < 40) or (t * sqrt(v / n) > w * sqrt(n - 1)))

where w is the required margin of error (half-width of the confidence interval).
In your implementation, you compare the current margin of error to (a function of) the current average:
t * sqrt(currentSumVariance / (n - 1)) > c * currentAvg * sqrt(n)

The details about the constants are not important (when do we use n-1 vs n?) What matters is that you have currentSumVariance on the left side and currentAvg on the right side of the inequality. So both sides of the comparison are stochastic. currentAvg is the running sample mean, so effectively, the target precision varies from iteration to iteration, esp. at the start of the simulation. Most likely your simulations happened to exit the while loop too early due to the unintended randomness in the required precision / margin of error.
You should keep w fixed.
Aside: This method for computing the sample variance one step at a time, rather than with a big sum over the entire sample, is known as Welford accumulator.
To verify, I coded the algorithm in python. (For reproducibility, I fixed the seed.) There are 13 outliers in 1,000 simulations; this agrees well with the significance level α = 0.01.
[1] L. M. Leemis and S. K. Park. Discrete-Event Simulation: A First Course (2006)
import numpy as np
import scipy.stats as stats

np.random.seed(seed=20220608)

def rexp(rate):
    return stats.expon.rvs(scale=1 / rate, size=1)[0]


def qt(p, df):
    return stats.t.ppf(1 - alpha / 2, df - 1)


def t_margin_of_error(n, alpha, sigma2):
    return qt(1 - alpha / 2, n - 1) * np.sqrt(sigma2 / (n - 1))


def run_simulation(alpha, delta, rate, verbose = False):

    currentAvg = 0.0  # estimate
    currentSumVariance = 0.0  # residual sum of squares

    n = 0

    while n < 100 or t_margin_of_error(n, alpha, currentSumVariance / n) > delta:

        diff = rexp(rate) - currentAvg

        currentAvg += diff / (n + 1)
        currentSumVariance += diff * diff * n / (n + 1)

        n += 1

    half_width = t_margin_of_error(n, alpha, currentSumVariance / n)

    if verbose:
        print(f"CI = {currentAvg:.3f} ± {half_width:.3f}")
        print(f"Precision = {half_width / currentAvg * 100:.2f}%")
        print(f"Iterations = {n}")
        print(f"Population mean = {1/rate:.3f}")
        print("===========")

    return np.abs(1 / rate - currentAvg) > half_width


alpha = 0.01
delta = 0.3
rate = 0.5

outliers = 0
draws = 1000

for _ in range(draws):
    outliers += run_simulation(alpha, delta, rate)

outliers, draws
#> (13, 1000)

