Simulation Based on Michel Fodje's epr-simple simulation translated from 
Python to Mathematica by John Reed 13 Nov 2013 
Modified by Fred Diether for Completely Local-Realistic Sep 2021
Some parts by Bill Nelson.

Set Run Time Parameters, Initialize Arrays and Tables

m=10000;
trialDeg=721;
s1=ConstantArray[0,m];
s2=ConstantArray[0,m];
λ1=ConstantArray[0,m];
A12=Table[{0,0,0,0},m];
A22=Table[{0,0,0,0},m];
B12=Table[{0,0,0,0},m];
B22=Table[{0,0,0,0},m];
a1=ConstantArray[0,m];
b1=ConstantArray[0,m];
A=ConstantArray[0, m];
B=ConstantArray[0, m];
nPP=ConstantArray[0,trialDeg];
nNN=ConstantArray[0,trialDeg];
nPN=ConstantArray[0,trialDeg];
nNP=ConstantArray[0,trialDeg];
nAP=ConstantArray[0,trialDeg];
nBP=ConstantArray[0,trialDeg];
nAN=ConstantArray[0,trialDeg];
nBN=ConstantArray[0,trialDeg];
ϕ=3.1; β=0.32; ξ=-16;  (*Adustable parameters for fine tuning*)

Generating Particle Data with Three Independent Do-Loops

Do[θ=RandomReal[{0,360}];  (*Singlet vector angle*)   (*Hidden Variable*)
λ1[[i]]=β(Cos[(θ/ϕ) ]^2);
s1[[i]]=θ;
s2[[i]]=θ+180,{i,m}]          (*Conservation of angular momentum*)

Do[a=RandomInteger[{0,360}];  (*Detector vector angle 1 degree increments*)
If[Abs[Cos[(a-s1[[i]])Degree]]>λ1[[i]],Aa1=-Sign[Cos[(a-s1[[i]])Degree]],Aa1=0];
If[Abs[Cos[(a-s1[[i]])Degree]]<λ1[[i]],Aa2=-Sign[Sin[(a-s1[[i]]+ξ)Degree]],Aa2=0];
A0=-Sign[Sin[(a-s1[[i]]+ξ)Degree]];
A12[[i]]={a,Aa1,i,A0};
A22[[i]]={a,Aa2,i,A0},{i,m}]
outA1=DeleteCases[A12,{_,0,_,_}];  (*Split into outA1 and outA2*)
outA2=DeleteCases[A22,{_,0,_,_}];

Do[b=RandomInteger[{0,360}];  (*Detector vector angle 1 degree increments*)
If[Abs[Cos[(b-s2[[i]])Degree]]>λ1[[i]],Bb1=-Sign[Cos[(b-s2[[i]])Degree]],Bb1=0];
If[Abs[Cos[(b-s2[[i]])Degree]]<λ1[[i]],Bb2=-Sign[Sin[(b-s2[[i]]+ξ)Degree]],Bb2=0];
B0=-Sign[Sin[(b-s2[[i]]+ξ)Degree]];
B12[[i]]={b,Bb1,i,B0};
B22[[i]]={b,Bb2,i,B0},{i,m}]
outB1=DeleteCases[B12,{_,0,_,_}];      (*Split into outB1 and outB2*)
outB2=DeleteCases[B22,{_,0,_,_}];

Matching Events Observed by Alice and Bob using Trial Numbers

(Subscript[k, B1])=outA1[[All,3]];  (*Two lists of only trial numbers used for matching*)
(Subscript[k, A1])=outB1[[All,3]];

Local Detection Analysis of the Events Observed by Alice

listA3=Select[outA1,Intersection[{#[[3]]},(Subscript[k, A1])]!={#[[3]]}&];
M1=Length[listA3];
tna=listA3[[All,3]];
outA4=Select[outA1,Intersection[{#[[3]]},tna]!={#[[3]]}&];
outA5=Table[{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]]},{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]]},(Subscript[k, B1])]!={#[[3]]}&];
M2=Length[listB3];
tnb=listB3[[All,3]];
outB4=Select[outB1,Intersection[{#[[3]]},tnb]!={#[[3]]}&];
outB5=Table[{0,0,0,0},M2];
sscb=ConstantArray[0,M2];

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]]},{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]];

m4=Length[outA5];
listA6=Table[{0,0,0,0},m4];
Do[If[ssca[[i]]==-1,listA6[[i]]=outA5[[i]],listA6[[i]]=0],{i,m4}]
listA7=DeleteCases[listA6,0];
listA83=listA7[[All,3]];
listB9=Select[outB,Intersection[{#[[3]]},listA83]=={#[[3]]}&];
m5=Length[listB9];
listB9=Table[{0,0,0,0},m5];
bb1=ConstantArray[0,m5];
Do[If[Abs[Cos[(listB9[[i]][[1]]-s2[[i]])Degree]]<λ1[[i]],bb1[[i]]=1,bb1[[i]]=0],{i,m5}]
N[Total[bb1]/m]*100

So that is the percentage of possible sign changes in A due to a change of angle b.  
This is “possible” events.  We don’t know for sure if they actually are due to a change
 of angle b.  Now, we need to average 10 of these.

(0.35+0.61+0.52+0.53+0.58+0.5+0.59+0.51+0.59+0.61)/10

100-0.539

So, we are for sure about 99.461 percent local for those that don’t accept the spinorial
sign changes as the reason.  The changes in B are probably about the same.