Simulate Equivalence Test I am trying to simulate a paired t-test equivalence test.  I generate x and y from the same Normal distribution thus I expect to reject the null hypothesis 95% of the time given an alpha of 0.05.
My question is what is the correct threshold (or boundary) to accurately simulate this?  I have calculated it as sqrt(2)/sqrt(n)*1.96 given my setup but this bound is too tight of a range based on the simulation.
library(equivalence)

N_sim=5000
reject=rep(0,N_sim)
n=1000
threshold=sqrt(2)/sqrt(n)*1.96

for(i in 1:N_sim){
  
  x=rnorm(n,0,1)
  y=rnorm(n,0,1)

  fit=tost(x, y, epsilon = threshold, paired = TRUE,alpha=.05)
  if(fit$tost.p.value<.05) reject[i]=1
  
}

mean(reject)

This equals about 0.25 (expect 0.95).
Update based on Ray's answer.
Below is the code to find the appropriate bound.
e_seq=seq(.14,.17,by=.001)
n=1000
N_Sim=2000
save_dat=list()
save_dat_idx=1

for(e in e_seq){
  in_interval=rep(0,N_Sim)
  for(i in 1:N_Sim){
    x=rnorm(n,0,1)
    y=rnorm(n,0,1)
    fit = t.test(x,y,paired=TRUE,conf.level = .9)$conf.int
    if(fit[1]>=-e & fit[2]<=e) in_interval[i]=1
  }
  save_dat[[save_dat_idx]]=tibble(e,mn=mean(in_interval))
  save_dat_idx=save_dat_idx+1
}

a1=do.call('rbind',save_dat)
a1

 A: It might be possible to find a threshold which will give an equivalent answer to a t-test approach, but it may not be as straighforward as it first appears. In a sense the null hypotheses of the two tests are opposites. For the t-test it is equality of means, for the equivalence test it is non-equivalence or dissimilarity. So using a theoretical CI limit as a threshold probably will not work as expected.
tost calculates a (t-) confidence interval for the mean difference between the supplied data sets, and if this interval lies within (-epsilon, epsilon), then the null hypothesisis (of non-equivalence) is rejected.
With alpha=.05, it looks like the tost function actually calculates a 90% confidence interval, rather than the expected 95% interval. So with your threshold for a 95% interval it is probably not surprising that you get a low value of 0.25 or 25% rejected. (I ran a little simulation for a standard normal variable and I got similar percentage for the 90% CI lying within the theoretical 95% CI.)
If you use the equivalent epsilon for a 90% z-interval CI (1.64*sqrt(2/n)), which would seem to be the logical choice, the rejection rate will be substantially lower.
The question is, would you expect the calculated CI to fall entirely within the theoretical CI limit? Unless the data has pretty much mean approx. 0 and SD approx. 1, this is very unlikely to happen, and hence the very low rejection rates.
I would think that the approach needed is to calculate a theoretical epsilon which would approach the 95% rejection rate based on the assumption that the data are from standard normal distributions:
Calculate $\epsilon$ such:
$P(\mbox{95% Conf. Int.} \in (-\epsilon, \epsilon) | x, y \sim N(0,1)) = 0.95$
