% MaxCorr_Sep2017.m : Correlation bounds for various distributions
%               Tools for computation and simulation
%               Version 1.0,
%               (c) Sergei Leonov (2017)
%
% This code is free for personal research use.
% Please cite:
%
%   Leonov, S. & Qaqish, B. (2017).
%   Correlated endpoints: simulation, modeling, and extreme correlations.
%   Stat Papers. https://doi.org/10.1007/s00362-017-0960-2
%
% For commercial use, please contact the author.
%
%       -------------------------------------
%       Testing max correlation between r.v.
%       -------------------------------------
%
%       Matlab code by Sergei Leonov, September 2017
%
%           --------
%            Inputs
%           --------      
%
%   1. s - random seed (to be defined outside, in the command window)
%   2. ExNumber - example number (to be changed below)
%   3. N - number of summands in the series for Example 9 (normal/Poisson
%                                                              or NegBin)

%   ====    September 2017 update, for "Statistical Papers"

fprintf(' 1. Maximal correlation (closed form)\n')

N = 40;     % May need changing for Example 9 (Normal - Poisson or Negative Binomial)

%   ExNumber refers to Examples' indices
ExNumber = 18;  %  9:  Normal - Poisson, 
               % 91:  Normal - NegBin, 
               % 10: Normal - Bernoulli
               % 16: two Bernoulli, 
               % 18: Bernoulli - Poisson

S = 0;      % Preliminary calculations for Examples 9/9.1              
switch ExNumber               
  case 9
     Lam1 = 1.5;
     Lam12 = sqrt(Lam1); % St. dev. of Poisson
     for i = 1:N,
        S = S + normpdf(norminv(poisscdf(i-1,Lam1)));
     end;
     fprintf(' Ex. 9, Normal-Poisson: lammda = %5.3f\n',Lam1)
     fprintf('  Number of summands in the series N = %4.0f\n',N)     
  case 91
     pS = 0.3;   
     NoS = 4;
     Lam12 = sqrt((1-pS)*NoS/(pS*pS));   % St. Dev. of Neg. binomial, bottom of p.10
     for i = 1:N,
        S = S + normpdf(norminv(nbincdf(i-1,NoS,pS)));
     end;
     fprintf(' Ex. 9, Normal-Neg. binomial: Pr(success) = %5.3f, nSucc =%5.0f\n',...
        pS,NoS);
     fprintf('  Number of summands in the seris N = %4.0f\n',N)          
  case 10    % Normal-Bernoulli
     pS = 0.3;
     Lam12 = sqrt(pS*(1-pS));       % Denominator in formula (19)
     S = normpdf(norminv(1-pS));    % Numerator in formula (19)
     fprintf(' Example 10, Normal-Bernoulli: p = %5.3f\n',pS);
  case 16    % Two Bernoulli, formula (29)
     pS = 0.7;
     pS2 = 0.4;
     Lam12 = sqrt(pS*(1-pS)*pS2*(1-pS2));
     S = min(pS,pS2) - pS*pS2;
     fprintf(' Example 16, Two Bernoulli: p1 = %5.3f, p2 = %5.3f\n',...
         pS,pS2);     
  case 18    % Bernoulli and Poisson, formula (31)
     pS = 0.2;
     Lam1 = 3;
     Lam12 = sqrt(pS*(1-pS)*Lam1);  % Product of variances
     ii = poissinv(1-pS,Lam1);
     S = ii*(poisscdf(ii,Lam1)-(1-pS));
     for i = ii+1:N,
         S = S + i*poisspdf(i,Lam1);
     end;
     S = S - Lam1*pS;
     fprintf(' Example 18, Bernoulli-Poisson, p = %5.3f, lambda = %5.3f\n',...
        pS,Lam1);    
end;

%   ----------------
    CorrP = S/Lam12;
%   ----------------

fprintf('  MaxCorr = %6.4f\n',CorrP)


%           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%   ======= Now check maximal correlation via Monte Carlo/Whitt
%           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^


%       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%       s: random seed initialization
%       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
rng(s)      % Analogous to former seed/state options
            
NN = 1000000;           % Number of simulations
zz = rand(NN,1);        % Uniform r.v.
X2 = norminv(zz);
X21 = normcdf(X2);

switch ExNumber
    case 9           % Normal-Poisson
        X1 = poissinv(zz, Lam1);     % X1 – Poisson component, inverse 
        X11 = poisscdf(X1,Lam1);     % For Spearman, use cdf of X1
    case 91          % Normal-Negative binomial
        X1 = nbininv(zz,NoS,pS); 
        X11= nbincdf(X1,NoS,pS);
    case 10          % Normal-Bernoulli
        X1 = binoinv(zz,1,pS);
        X11= binocdf(X1,1,pS);
    case 16          % Two Bernoulli
        X1 = binoinv(zz,1,pS);
        X2 = binoinv(zz,1,pS2);
        X11= binocdf(X1,1,pS);
        X21= binocdf(X2,1,pS2);        
    case 18          % Bernoulli-Poisson
        X1 = binoinv(zz,1,pS);
        X2 = poissinv(zz,Lam1);
        X11= binocdf(X1,1,pS);
        X21= poisscdf(X2,Lam1);        
end;

corrMC = corr(X1,X2);
fprintf('\n 2.    Numerical (Monte Carlo/Whitt): Corr = %6.4f\n',corrMC)

corrS = corr(X11,X21);
fprintf(' 2(a).    Spearman Corr (via cdfs) = %6.4f\n',corrS)

