Simulation Based on Michel Fodje’s  epr-simple and
Joy  Christian’s Updated 3-Sphere Model.
Parts of Quaternion and Matching Code by John Reed.
Assembled by Fred Diether, Jan. 2022, for CHSH Analysis.

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;
λ=ConstantArray[0,m];
ss1=ConstantArray[0,m];
Ls1=ConstantArray[0,m];
Ls2=ConstantArray[0,m];
outqA=Table[{0,0,0,0},m];
outqB=Table[{0,0,0,0},m];
outA12=Table[{0,0,0,0},m];
outB12=Table[{0,0,0,0},m];
outA22=Table[{0,0,0,0},m];
outB22=Table[{0,0,0,0},m];
a2=ConstantArray[0,m];
b2=ConstantArray[0,m];
h1=ConstantArray[0,m];
h2=ConstantArray[0,m];
ϕ=3; β=0.25; ξ=-15;

Generating Particle Data with Three Independent Do-Loops

Do[s=RandomPoint[Sphere[]];  (*Singlet 3D vector*)   (*Hidden Variable*)
θ1=ToSphericalCoordinates[s][[3]]*180/π;
ss1[[i]]=θ1;
λ[[i]]=β(Cos[((θ1)/ϕ) ]^2);  (*Hidden variable mechanism*)
Ls1[[i]]=s.Qcoordinates;  (*Convert to quaternion coordinates*)(*A particle*)
Ls2[[i]]=-s.Qcoordinates,{i,m}]      (*B particle*)

Do[a=RandomChoice[{0,90}];  (*Detector 2D vector angle*)
a2[[i]]=a;
aa=N[Flatten[{FromPolarCoordinates[{1,a *π/180}],0.000001}]];
Da=aa.Qcoordinates;  (*Convert to quaternion coordinates*)
qa=Da**Ls1[[i]];
If[Abs[Re[qa]]>λ[[i]],qA=Sign[Re[qa]],qA=Sign[Sin[(a-ss1[[i]]+ξ)Degree]]];
A0=Sign[Sin[(a-ss1[[i]]+ξ)Degree]];
outqA[[i]]={a,qA,i,A0};
If[Abs[Re[qa]]>λ[[i]],outA12[[i]]=outqA[[i]],outA22[[i]]=outqA[[i]]],{i,m}]
outA2=DeleteCases[outA22,{0,0,0,0}];
outA3=outqA;
Do[A2tn=outA2[[i]][[3]];  (*Trial numbers from outA2*)
h1[[A2tn]]=1,{i,Length[outA2]}]

Do[b=RandomChoice[{45,135}];  (*Detector 2D vector angle*)
b2[[i]]=b;
bb=N[Flatten[{FromPolarCoordinates[{1,b *π/180}],0.000001}]];
Db=bb.Qcoordinates;  (*Convert to quaternion coordinates*)
qb=Ls2[[i]]**Db;
If[Abs[Re[qb]]>λ[[i]],qB=Sign[Re[qb]],qB=-Sign[Sin[(b-ss1[[i]]+ξ)Degree]]];
B0=-Sign[Sin[(b-ss1[[i]]+ξ)Degree]];
outqB[[i]]={b,qB,i,B0};
If[Abs[Re[qb]]>λ[[i]],outB12[[i]]=outqB[[i]],outB22[[i]]=outqB[[i]]],{i,m}]
outB2=DeleteCases[outB22,{0,0,0,0}];
outB3=outqB;
Do[B2tn=outB2[[i]][[3]];  (*Trial numbers from outB2*)
h2[[B2tn]]=1,{i,Length[outB2]}]

Spinorial Sign Changes in A and B
For the spinorial sign changes we will need,

Do[If[h2[[i]]==1&&outqA[[i]][[2]]≠outqA[[i]][[4]],outA3[[i]][[2]]=outqA[[i]][[2]]*-1],{i,m}]  (*Spinorial sign change*)
A=outA3[[All,2]];

Do[If[h1[[i]]==1&&outqB[[i]][[2]]≠outqB[[i]][[4]],outB3[[i]][[2]]=outqB[[i]][[2]]*-1],{i,m}]   (*Spinorial sign change*)
B=outB3[[All,2]];

CHSH Analysis of the Particle Data Received from Alice and Bob

nP1=0; nN1=0; nP2=0; nN2=0; nP3=0; nN3=0;nP4=0;  nN4=0;
Do[a1=a2[[i]];b1=b2[[i]];
aliceD=A[[i]];bobD=B[[i]];
If[(b1==45)&&(a1-b1==-45)&&aliceD*bobD==1,nP1++];
If[(b1==45)&&(a1-b1==-45)&&aliceD*bobD==-1,nN1++];If[(a1-b1)==-135&&aliceD*bobD==1,nP2++];
If[(a1-b1)==-135&&aliceD*bobD==-1,nN2++];If[(a1-b1)==45&&aliceD*bobD==1,nP3++];
If[(a1-b1)==45&&aliceD*bobD==-1,nN3++];
If[a1==90&&(a1-b1)==-45&&aliceD*bobD==1,nP4++];
If[a1==90&&(a1-b1)==-45&&aliceD*bobD==-1,nN4++],{i,m}]
E1=N[(nP1-nN1)/(nP1+nN1)];
E2=N[(nP2-nN2)/(nP2+nN2)];
E3=N[(nP3-nN3)/(nP3+nN3)];
E4=N[(nP4-nN4)/(nP4+nN4)];
tot1=Total[nP1+nP2+nP3+nP4+nN1+nN2+nN3+nN4];
CHSH=N[Abs[E1-E2+E3+E4]];
Print["CHSH = ",CHSH]
Print["Total Events Detected = ",tot1]

(2.778545010598531+2.8292076083106292+2.8335071953985898+2.8260985341778437+2.8254732463559082+2.8043705677375788+2.800333098716658+2.7828705966805343+2.7640294414490336+2.8108636228331267)/10  
(*Average of 10 Runs of 20000 Trials Each to Machine Precision*)