NNNNNNNNNNNNNNNN code NNNNNNNNNNNNNNNNNNN

H[p_]:=-Chop[p.Log[p+$MinMachineNumber]]; (* entropy of the distribution p *)
RH[p_,q_]:=Chop[p.(Log[p+$MinMachineNumber]-Log[q+$MinMachineNumber])]; (*relative entropy p vs q*)
RHB[p_,α_]:= (*relative entropy p vs Binomial(α)*)p.(Log[p+$MinMachineNumber]-Table[Log[Binomial[Length[p]-1,k]]+k Log[Min[α,1-α]]+(Length[p]-1-k)Log[Max[α,1-α]],{k,0,Length[p]-1}]);

MIs[α_,l_,NE_:$MaxMachineNumber,quantiles_:{},print_:False]:=Module[{i,T,dim,z,a,c},
(* Computes the (system) CTI: If NE ≤ 0 it computes the CTI κ(α), if NE ≥ sample size T it computes the system CTI (OverTilde[κ])(α)*)
z=If[Length[quantiles]==0,Quantile[#,α,{{1/2,0},{0,1}}]& /@ Transpose[l],quantiles];T=Length[l];dim=Length[l[[1]]];a=α;
z=Transpose[z-#&/@ l];If[α> .5,a=1-α;z=-z];c=Clip[((Log[T]-Log[Max[NE,$MinMachineNumber]]+dim Log[1-a])/(Log[1/a-1.]+$MinMachineNumber)),{0,T}];
z=Clip[ Sign[#],{0,1}]&/@ z;z=Transpose[z]; 
z=If[c<Total[#],{Total[#]},#] &/@z; 
z=Tally[(#&/@z),#1==#2&];z={Total[#[[1]]],#[[2]]/T} &/@z;
If[print,Print["T = ",T,", α = ",a,", dim = ",dim,", c= ",c,", z= ",1.z]];
z=If[c<#[[1]],{#[[2]], Log[Binomial[dim,#[[1]]]]+#[[1]]Log[a]+(dim-#[[1]])Log[1-a]},
{#[[2]],#[[1]]Log[a]+(dim-#[[1]])Log[1-a]}] &/@z;z=Transpose[z];
Return[(( z[[1]].(Log[z[[1]]]-z[[2]]))/1) ] (*/(dim-1)H[{α,1-α}]*)];

ExceedL[α_,l_,NE_:0,quantiles_:{},print_:False]:=Module[{T,dim,z,a},
(* Transforms the list l into a list of α-JTs indicators if NE ≤ 0 or a list of # of α-exceedences if NE ≥ T*)
z=If[Length[quantiles]==0,Quantile[#,α,{{1/2,0},{0,1}}]& /@ Transpose[l],quantiles];
T=Length[l];dim=Length[l[[1]]];a=α;
z=Transpose[z-#&/@ l];If[α> .5,a=1-α;z=-z];
If[print,Print["T = ",T,", α = ",α,", dim = ",dim]];
z=Clip[ Sign[#],{0,1}]&/@ z;z=Transpose[z]; 
z=If[T*(a^(Total[#]))(((1-a))^(dim-Total[#]))<NE,{Total[#]},#] &/@z; 
Return[z]];

(**************************** EWMA(λ)  ********************************)
JTPrλs[α_,l_,λ_:1,NE_:0,quantiles_:{},print_:False]:=Module[{T,dim,z,λT},
(* Computes the EWMA(λ) probabilities of JTs if NE ≤ 0 and of # of exceedances if NE ≥ T*)
T=Length[l];dim=Length[l[[1]]];λT=If[λ≥ 1,(1./T),((1.-λ)/(1-(λ^T)))];
z=ExceedL[α,l,NE,quantiles,print];
z=Normal[PositionIndex[z]]; z={#[[1]],λT Chop[ Total[ ((Max[λ,$MachineEpsilon])^(#[[2]]-1))]]}&/@z;
Return[Select[z,#[[2]]>0&]]];

MIλs[α_,l_,λ_:1,NE_:0,quantiles_:{},print_:True]:=Module[{T,dim,z,a,λT},
(* Computes first the EWMA(λ) probabilities, where all cells with less than NE expected obs are lumped together, and then the MI
 NE ≤ 0 computes κ(α), NE ≥ T computes (OverTilde[κ])(α)*)
T=Length[l];dim=Length[l[[1]]];a=Min[1-α,α];
z=JTPrλs[α,l,λ,NE,quantiles,print]; z={Total[#[[1]]],#[[2]]}&/@z;
z=If[ T*(a^(#[[1]]))(((1-a))^(dim-#[[1]]))<NE,{#[[2]], Log[Binomial[dim,#[[1]]]]+#[[1]]Log[a]+(dim-#[[1]])Log[1-a]},
{#[[2]],#[[1]]Log[a]+(dim-#[[1]])Log[1-a]}] &/@z;z=Transpose[z];
Return[ z[[1]].(Log[z[[1]]]-z[[2]])] (*/(dim-1)H[{a,1-a}]*)];

MLλ[α_,l_,Minλ_:0.9,Maxλ_:1,ini_:$MaxMachineNumber,quantiles_:{},print_:True]:=Module[{i,T,dim,z,λ,p0,p,pM,LogL,MinLL,λλ,f},
(* ML estimate of λ in EWMA p(t+1)= λp(t)+(1-λ)u(t) *)
T=Length[l];dim=Length[l[[1]]];λλ=Minλ;MinLL=-$MaxMachineNumber;
p0=JTPrλs[α,l[[;;Min[ini,T]]],1,T,{},print];p0=Flatten[#]&/@p0;
p0=Sort[p0~Join~Table[{i,0},{i,0,dim}],#1[[1]]<#2[[1]]&];
p0=GatherBy[p0,First];p0=Total[#]&/@p0;p0=#[[2]]&/@ p0;pM=p0;

If[Length[p0]≠ dim+1,Print["Error: Length[p]≠ dim+1 "];Abort[]];z=Flatten[ExceedL[α,l,T,{},False]];
Monitor[
For[λ=Minλ,λ≤Maxλ,λ=λ+Max[$MachineEpsilon,(Maxλ-Minλ)/100],LogL=0;p={p0};
For[i=1,i≤T,i++,
LogL=LogL+Log[$MachineEpsilon+(Last[p][[z[[i]]+1]])];
(*f=(((dim-1)α)/(Max[z[[i]]-1,(dim-1)α])); f={1-f+f(((1-(a/f)))^dim)}~Join~Table[f PDF[BinomialDistribution[dim,(a/f)],s],{s,1,dim}];factor probability f= (((dim-1)a)/(k-1)) ∈ [a,1]  ML estimated for k exceedances*)
p=Append[p,λ Last[p] +(1-λ)UnitVector[dim+1,z[[i]]+1]]];
If[LogL>MinLL,MinLL=LogL;λλ=λ,pM=Last[p]]],{λ,λλ,MinLL}];
If[Minλ==Maxλ,Return[((RHB[#,α])/((dim-1)H[{α,1-α}]))&/@ p]];
Return[{λλ,MinLL,LogL,Range[0,dim].pM,Range[0,dim].p0}]];

Dynλs[α_,l_,λ_,υ_:1,ini_:$MaxMachineNumber,quantiles_:{},print_:True]:=Module[{t,T,dim,z,p,κ,v},
(* List of dynamic (OverTilde[κ])(α) and Pr[# α-exceed > dim-v] for EWMA(λ) PDF of the # of exceedances*)
T=Length[l];dim=Length[l[[1]]];κ=v={};
p=JTPrλs[α,l[[;;Min[ini,T]]],1,T,{},print];p=Flatten[#]&/@p;
p=Sort[p~Join~Table[{i,0},{i,0,dim}],#1[[1]]<#2[[1]]&];
p=GatherBy[p,First];p=Total[#]&/@p;p=#[[2]]&/@ p;If[print,Print["p= ",p]];
If[Length[p]≠ dim+1,Print["Error: Length[p]≠ dim+1 "];Abort[]];z=Flatten[ExceedL[α,l,T,{},print]];
For[t=1,t≤ T,t++,
κ=AppendTo[κ,RHB[p,α]]; (*(OverTilde[κ])(α) forecast for period t*)
v=AppendTo[v,Total[p[[-υ;;]]]];(* Pr[# α-exceed > dim-v] forecast for period t*)
p= p +(1-λ)UnitVector[dim+1,z[[t]]+1];
p=p Min[ ((dim α)/(Range[0,dim].p)),1];p[[1]]=Max[1-Total[Rest[p]],0]];
Return[(({κ,v})/((dim-1)H[{α,1-α}]))]];

(**************************** EWMA(λ)  ********************************)
SysCI[α_,l_,iter_:1000,minbin_:0,print_:False]:=Module[{x,i,j,T,dim,z,Bi,k},
(* Computes the bootstrap 5%-CI for the sys α-CTI *)
T=Length[l];dim=Length[l[[1]]];x=l;
Bi=Table[PDF[BinomialDistribution[dim,Min[α,1-α]],i],{i,0,dim}];j=0;k={};
i=LengthWhile[Sort[Bi,Greater],#>minbin/T&];
If[print,Print["T = ",T,", α = ",α,", dim = ",dim,", Pr cutoff at = ",i]];
Monitor[
While[j<iter,
z=Flatten[#]&/@JTPrλs[α,x,1,T];
z=Sort[z~Join~Table[{i,0},{i,0,dim}],#1[[1]]<#2[[1]]&];
z=GatherBy[z,First];z=Total[#]&/@z;z=#[[2]]&/@z;z=Transpose[{z,Bi}];z=Sort[z,#1[[2]]>#2[[2]]&][[;;i ]];
If[i≤ dim,z=Append[z,Chop[1.-Total[z]]]];z=Transpose[z];
k=Append[k, ((z[[1]].(Log[$MachineEpsilon+z[[1]]]-Log[z[[2]]]))/((dim-1)H[{α,1-α}])) ];x=RandomChoice[l,T];j++],{j,k[[j]]}];
z={k[[1]],Quantile[k,.025,{{1/2,0},{0,1}}],Quantile[k,1-.025,{{1/2,0},{0,1}}]};
Return[z~Join~{z[[3]]-z[[1]]}]];

CRX2[T_,dim_:2,q_:0.01,αlist_:ticks50]:=((Quantile[ChiSquareDistribution[NrNonEmptyCells[T,dim,#]-dim-1],1.-q])/(2T*(dim-1) H[{#,1.-#}]))&/@αlist;
(* (Chi^2) q% critical region for the κ(α) statistic for sample size T, dimension dim and all α in αlist. For example, CR[10000]*)

CRX2compacted[T_,dim_:2,q_:0.01,αlist_:ticks50]:=((Quantile[ChiSquareDistribution[dim-1],1.-q])/(2T*(dim-1) H[{#,1.-#}]))&/@αlist;
(* (Chi^2) q% critical region for the compacted κ(α) statistic for sample size T, dimension dim and all α in αlist. For example, CR[10000]*)


Tallyz[α_,l_,quantiles_:{}]:=Module[{z},(* counts obs in the list l falling into different JTs *) 
z=If[Length[quantiles]==0,Quantile[#,α,{{1/2,0},{0,1}}]& /@ Transpose[l],quantiles];z=Transpose[z-#&/@ l];
If[α≤ .5,z=Total[Clip[#,{0,0},{0,1}]&/@ z],z=Total[Clip[#,{0,0},{0,1}]&/@ -z]];
z=Tally[z~Join~Range[0,Length[l[[1]]]],Equal];z=Sort[z,#1[[1]]<#2[[1]]&];z=Transpose[z][[2]]-1;
If[Total[z]≠Length[l],Print["****** ∑z= ",Total[z]," ≠ Length[l]= ",Length[l]," .... Aborting..."];Abort[],Return[z]]];

CondSymBin[α_,l1_,l2_,max_:1000]:=Module[{dim,p,len1,len2,z1,z2,z,x,randMH,chi2},
(* conditional (generalized Fisher) test of asymmetry between two tail interdep structures (TIS). Usage: CondAsym[.4,x,y]
To test the upper/lower tail symmetry: CondAsym[.4,x,-x] *) 
dim=Length[l1[[1]]];If[dim≠ Length[l2[[1]]],Print["******* dim l1 ≠ dim l2 *******"];Abort[]];
z1=Tallyz[α,l1];z2=Tallyz[α,l2];z=z1+z2;If[Total[z1]≠ Total[z2],Print["***** ∑l1= ",Total[z1]," ≠ ∑l2 = ",Total[z2]," *****"]];
If[Min[z]≤ 0.,z1=Pick[z1, Sign[z],1];z2=Pick[z2, Sign[z],1];z=Pick[z, Sign[z],1] ];
len1=Total[z1];len2=Total[z2];p=1.z/Total[z];

chi2=Total[(((z1-len1 p))^2)/(len1 p)]+Total[(((z2-len2 p))^2)/(len2 p)];

x=PDF[MultivariateHypergeometricDistribution[len1,z],z1];
randMH=RandomVariate[MultivariateHypergeometricDistribution[len1,z],max];randMH=PDF[MultivariateHypergeometricDistribution[len1,z],#]& /@ randMH;
Return[{1.BinCounts[-randMH,{{-max,-x,max}}][[2]]/max,{chi2,1.-Quiet[CDF[ChiSquareDistribution[Max[(2-1)Length[z]-1,1]],chi2]]}}]];

NNNNNNNNNNNNNNNNNNNN        Synthetic data         NNNNNNNNNNNNNNNNNNNNNNNNNNN

ρ=0.7;dim=4;T=2500;x=RandomReal[MultinormalDistribution[Table[0,{dim}],ρ Table[1,{dim},{dim}]+(1-ρ)IdentityMatrix[dim]],{T}];

MIs[.2,x,0]

MIs[.2,x,Length[x]]

(*NNNNNNNNNNNNNNNNNNNN     SP500   NNNNNNNNNNNNNNNNNNNNNNNNNNN*)

r=Import["C:\UEA\projects\Tail Interdependence\data\SP500\\RAWRETS_SP500.csv","CSV"][[2;;]];(*RAWRETS_SP500, ARGARCHResids_SP500, FFCResids_SP500, FFCGARCHResids_SP500*)

Mean[r][[1;;10]]

(*********************************  Contributions  ***************************************)

a=.1;zc={};maxi=Dimensions[r][[2]];
Monitor[For[i=1,i≤maxi,i++,
AppendTo[zc,MIs[a,r[[All,Delete[Range[maxi],i]]],Length[r]]]],{i}]; 
zc=MIs[a,r,Length[r]]-zc;

ListPlot[Sort[zc,Greater],PlotLegends->Placed[{"α=0.1","α=0.5","α=0.9"},Above],Joined-> True,PlotStyle->{Black,{Black,Dotted},Gray},AxesLabel->{"i","(OverTilde[MI])(α)-(Subscript[(OverTilde[MI]), -i])(α)"},PlotRange->Automatic]

cr=Import["C:\UEA\projects\Tail Interdependence\JE\Results & Diagrams SP500\\MI contributions to FFCGARCH.csv","CSV"][[2;;]];Dimensions[cr]

G1=ListPlot[Sort[#,Greater]&/@Transpose[cr],PlotLegends->Placed[{"α=0.1","α=0.5","α=0.9"},Above],Joined-> True,PlotStyle->{Black,{Black,Dotted},Gray},AxesLabel->{"i","(OverTilde[MI])(α)-(Subscript[(OverTilde[MI]), -i])(α)"},PlotRange->Automatic]

x[[Ordering[zc]]]