Logically, your exact approach would depend on the audience for your paper, the types of problems they routinely compute, and the specific software to which they have access and know how to use.
Also,
your paper may draw attention to seldom-used software that is freely available and very accurate. Or warn against popular software that often makes
serious computational errors for particular tasks.
So the useful scope of your paper may not be
entirely clear until some of your results are
known.
If there has not been a paper along these lines
since McCullough (1999), you might begin by
checking whether any important shortcomings noted there have been overcome in later releases of the software.
You begin your Question with, "[W]e are trying to establish numerical equivalence (within reasonable precision) for selected statistical models across programming languages such as SAS, R & Python or even different packages within same programming language. This will allow one to confidently use programming language/Package of their choice,..."
which makes good sense.
Note: I am pretty sure that McCullough (1999) is the paper I vaguely remembered. Glad you found it.
Addenda:
(1) Do R and Minitab give the same P-values
for the Welch t test, for (essentially)
the same two normal samples.
The Welch two-sample t test is widely used and included in many computer programs. Unlike the pooled 2-sample t test it
does not assume that the two samples come from populations with equal variances.
For the same data, do two computer programs give essentially the same results.
Consider fictitious data generated in R as follows; R first:
set.seed(1234)
x1 = rnorm(10, 50, 5)
x2 = rnorm(100, 55, 2)
summary(x1); sd(x1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
38.27 45.94 47.22 48.08 51.96 55.42
[1] 4.978938
summary(x2); sd(x2)
Min. 1st Qu. Median Mean 3rd Qu. Max.
50.64 53.32 54.39 54.73 55.74 60.10
[1] 1.932823
t.test(x1, x2)$p.val
[1] 0.002190574
Now Minitab:
Now, I input the sample sizes, means, and standard deviations
into a recent release of Minitab, which happens to be installed
on my computer as I type this (you should use more decimal places of accuracy and the latest release of Minitab). Also, Minitab is one of the few
programs that will directly accept such summarized data.]
Sample N Mean StDev SE Mean
1 10 48.08 4.98 1.6
2 100 54.73 1.93 0.19
Difference = μ (1) - μ (2)
Estimate for difference: -6.65
95% CI for difference: (-10.24, -3.06)
T-Test of difference = 0 (vs ≠):
T-Value = -4.19 P-Value = 0.002 DF = 9
The P-value agrees with the tha P-value from R to two places. (You could also compare the T statistics and DF.)
Difference = μ (1) - μ (2)
Estimate for difference: -6.65
T-Test of difference = 0 (vs ≠):
T-Value = -4.22 P-Value = 0.002 DF = 9
One might also discuss whether it is a good
idea to allow the use of summarized data, which cannot be tested for normality (or other assumptions), and Minitab's frequent use of 'sample' to mean observation.
(2) Bootstrapping. There are many styles of bootstrap confidence intervals for parameters and many of them have gained popularity since 2000. Even if two programs are based on the same method, you can't expect exactly the same result because bootstraps use random re-sampling of data. Moreover as this Answer suggests, the best bootstrap CI for a parameter (here the sample variance) may depend on what is known and what is assumed. What would be your standard of judging whether two bootstrap CIs are (essentially) "the same"?
