Simulation Based on Michel Fodje's epr-simple simulation translated from 
Python to Mathematica by John Reed 13 Nov 2013 and Quaternions
Modified by Fred Diether for Completely Local-Realistic Sep 2021
Some parts by Bill Nelson.  Includes Joy’s (S^3) Quaternion Model.

Load Quaternion Package, Set Run Time Parameters, Initialize Arrays and Tables

<<Quaternions`;
β0=Quaternion[1,0,0,0];
β1=Quaternion[0,1,0,0];
β2=Quaternion[0,0,1,0];
β3=Quaternion[0,0,0,1];
Qcoordinates={β1,β2,β3};
m=20000;
trialDeg=721;
Ls1=ConstantArray[0,m];
Ls2=ConstantArray[0,m];
λ1=ConstantArray[0,m];
λ2=ConstantArray[0,m];
outAa=Table[{0,0,0,0,0,0},m];
outBb=Table[{0,0,0,0,0,0},m];
outA12=Table[{0,0,0,0,0,0},m];
outA22=Table[{0,0,0,0,0,0},m];
outB12=Table[{0,0,0,0,0,0},m];
outB22=Table[{0,0,0,0,0,0},m];
plotq=Table[{0,0},m];
a1=ConstantArray[0,m];
b1=ConstantArray[0,m];
A=ConstantArray[0, m];
B=ConstantArray[0, m];
ϕ=3.1; β=0.32; ξ=-((4π)/45);  (*Adustable parameters for fine tuning*)

Generating Particle Data with Three Independent Do-Loops

Do[θ=RandomReal[{-179,180}];  (*Singlet vector angle*)   (*Hidden Variable*)
λ1[[i]]=β(Cos[(θ/ϕ) ]^2);
λ2[[i]]=Sign[θ];  (*Orientation Hidden Variable*)
θθ=N[Flatten[{FromPolarCoordinates[{1,θ *π/180}],0}]];
Ls1[[i]]=λ2[[i]]*θθ.Qcoordinates;  (*Quaternion particle for the A side*)
Ls2[[i]]=-Ls1[[i]],{i,m}]          (*Quaternion particle for the B side*)(*Conservation of angular momentum*)

Do[a=RandomInteger[{-179,180}];  (*Detector vector angle 1 degree increments*)
aa=N[Flatten[{FromPolarCoordinates[{1,a *π/180}],0}]];
Da=aa.Qcoordinates;  (*Convert to quaternion coordinates*)
qa=Da**Ls1[[i]];
aq=-Ls1[[i]]**Da;
If[Abs[Re[qa]]>λ1[[i]],Aa=Re[Da**Limit[ Ls1[[i]], Ls1[[i]]->Sign[Re[Da**Ls1[[i]]]]Da]],
Aa=λ2[[i]]*Sign[qa[[4]]+ξ]];
A0=λ2[[i]]*Sign[qa[[4]]+ξ];
outAa[[i]]={a,Aa,i,A0,qa,aq};
If[Abs[Re[qa]]>λ1[[i]],outA12[[i]]=outAa[[i]],outA22[[i]]=outAa[[i]]],{i,m}]
outA1=DeleteCases[outA12,{0,0,0,0,0,0}];  (*Split outAa into outA1 and outA2*)
outA2=DeleteCases[outA22,{0,0,0,0,0,0}];

Do[b=RandomInteger[{-179,180}];  (*Detector vector angle 1 degree increments*)
bb=N[Flatten[{FromPolarCoordinates[{1,b *π/180}],0}]];
Db=bb.Qcoordinates;  (*Convert to quaternion coordinates*)
qb=Ls2[[i]]**Db;
bq=-Db**Ls2[[i]];
If[Abs[Re[qb]]>λ1[[i]],B=Re[Db**Limit[ Ls2[[i]], Ls2[[i]]->Sign[Re[Db**Ls2[[i]]]]Db]],
B=-λ2[[i]]*Sign[qb[[4]]+ξ]];
B0=-λ2[[i]]*Sign[qb[[4]]+ξ];
outBb[[i]]={b,B,i,B0,qb,bq};
If[Abs[Re[qb]]>λ1[[i]],outB12[[i]]=outBb[[i]],outB22[[i]]=outBb[[i]]],{i,m}]
outB1=DeleteCases[outB12,{0,0,0,0,0,0}];      (*Split outBb into outB1 and outB2*)
outB2=DeleteCases[outB22,{0,0,0,0,0,0}];

Matching Events Observed by Alice and Bob using Trial Numbers

list13=outA1[[All,3]];  (*Two lists of only trial numbers used for matching*)
list23=outB1[[All,3]];

Local Detection Analysis of the Events Observed by Alice

listA3=Select[outA1,Intersection[{#[[3]]},list23]!={#[[3]]}&];
M1=Length[listA3];
tna=listA3[[All,3]];
outA4=Select[outA1,Intersection[{#[[3]]},tna]!={#[[3]]}&];
outA5=Table[{0,0,0,0,0,0},M1];
ssca=ConstantArray[0,M1];

For the spinorial sign changes we will need eq. (12).

Do[If[listA3[[i]][[2]]==listA3[[i]][[4]],δ=0,δ=1];
If[δ==0,ssca[[i]]=1,ssca[[i]]=-1];  (*spinorial sign change*)outA5[[i]]={listA3[[i]][[1]],ssca[[i]]*listA3[[i]][[2]],listA3[[i]][[3]],listA3[[i]][[4]],listA3[[i]][[5]],listA3[[i]][[6]]},{i,M1}]
outA=Sort[Catenate[{outA2,outA4,outA5}],#1[[3]]<#2[[3]]&];  (*Combine lists and sort*)
a1= outA[[All,1]];  (*These results are what Alice observes as defined in Eq.(??)*)
A =  outA[[All,2]];

Local Detection Analysis of the Events Observed by Bob

listB3=Select[outB1,Intersection[{#[[3]]},list13]!={#[[3]]}&];
M2=Length[listB3];
tnb=listB3[[All,3]];
outB4=Select[outB1,Intersection[{#[[3]]},tnb]!={#[[3]]}&];
outB5=Table[{0,0,0,0,0,0},M2];
sscb=ConstantArray[0,M2];
listB35=listB3[[All,5]];
listB36=listB3[[All,6]];
Do[If[listB3[[i]][[2]]==listB3[[i]][[4]],δ=0,δ=1];
If[δ==0,sscb[[i]]=1,sscb[[i]]=-1];  (*spinorial sign change*)outB5[[i]]={listB3[[i]][[1]],sscb[[i]]*listB3[[i]][[2]],listB3[[i]][[3]],listB3[[i]][[4]],listB3[[i]][[5]],listB3[[i]][[6]]},{i,M2}]
outB=Sort[Catenate[{outB2,outB4,outB5}],#1[[3]]<#2[[3]]&];  (*Combine lists and sort*)
b1= outB[[All,1]];  (*These results are what Bob observes as defined in Eq.(??)*)
B =  outB[[All,2]];

All the spinorial sign changes happen in outA5 and outB5 so we will get the percentage of events first from that.

N[100*Length[outA5]/m]
N[100*Length[outB5]/m]

Now we need to know  what the percentage is of actual sign changes.

ssca1=ConstantArray[0,M1];
Do[If[ssca[[i]]==-1,ssca1[[i]]=1,ssca1[[i]]=0],{i,M1}]
N[100*Total[ssca1]/m]

sscb1=ConstantArray[0,M2];
Do[If[sscb[[i]]==-1,sscb1[[i]]=1,sscb1[[i]]=0],{i,M2}]
N[100*Total[sscb1]/m]

So, there you have it.  We have 4.37 percent sign changes for A and 4.325 percent sign changes for B.
About the same as doing it the other way.