% g ~ uniform on (0,1);  
% e = a ~ exp(cbar)
% x = step of size 1 ... see text .. x has mean zero
%

clear all;
small = 1.0e-10;
rng(3418);

% -- File Directories   
figdir = 'fig/';
outdir = 'out/';
matdir = 'mat/';

% global variables
global GQxw;
GQxw = setGQxw;

ci_level = 0.95;
qmax = 30;       % Maximum value of q

n = 250;          % Number of obs per draw
nrep_e = 50000;    % Number of draws for e for each draws of s
nrep_sxz = 500;

rhobar = 0.03;

mod_str = ['mc_unitline_model_3_dm_' num2str(n)];

parms_misc.qmax = qmax;
parms_misc.ci_level = ci_level;
cv_normal = norminv(1-(1-ci_level)/2);

bvec1 = linspace(0.01,0.1,10)';
bvec2 = linspace(0.11,0.25,10)';
bvec3 = linspace(0.27,0.99,10)';
bvec = [bvec1;bvec2;bvec3];
nbvec = size(bvec,1);

rejfreq_scpc = NaN(nrep_sxz,1);
rejfreq_cscpc = NaN(nrep_sxz,1);
rejfreq_hr = NaN(nrep_sxz,1);
rejfreq_conley_mat = NaN(nrep_sxz,nbvec);

for is = 1:nrep_sxz
  % Draw s and compute distmat
  s_draw = rand(n,1);
  distmat = getdistmat_euclidean(s_draw);
  distmat = distmat/max(max(distmat));
  % Construct Conley Bartlet Matrices
  CMat = NaN(n,n,nbvec);
  for ib = 1:nbvec
      Q = dist_to_bartlet(distmat,bvec(ib));
      CMat(:,:,ib) = chol(Q);
  end
  cbar = getcbar(rhobar,distmat);
  cmax = getcbar(0.0001,distmat);
  % Draw a
  a = draw_a(cbar,distmat,nrep_e);
  % Draw X,Z,e
  [X,Z,e] = draw_XZe(a,s_draw);
  % Get results for SCPC
  [q,cv_scpc,r] = scpc_parms(distmat,cbar,cmax,qmax,1-ci_level);
 % Get results for cscpc
  if size(Z,1) > 1
    bZ = Z\X;
    x = X-Z*bZ;
   else
    x = X;
   end
   scl_X = sqrt(mean(x.^2));
   x = x/scl_X;  % Normalize so that MS = 1
   cv_eta = get_cv_eta(x,r,Z,cbar,cmax,distmat,1-ci_level);
   cv_cscpc = max([cv_scpc cv_eta]);
       
   % Compute tstats for all e draws
       
   V = x;
   if size(Z,1) > 1
     V = [x Z];
   end
   bV = V\e;
   ehat = e-V*bV;
   bx = (e'*x)/n;
   xehat = repmat(x,1,nrep_e).*ehat;
   pc_mat = r'*xehat;
   s_scpc = sqrt(mean(pc_mat.^2))';
   tstat = sqrt(n)*bx./s_scpc;
   sig_scpc = s_scpc/sqrt(n);
   s_hr = sqrt(mean(xehat.^2))';
   tstat_hr = sqrt(n)*bx./s_hr;
        
   % Save results
   rejfreq_scpc(is) = mean(abs(tstat) > cv_scpc); 
   rejfreq_cscpc(is) = mean(abs(tstat) > cv_cscpc);
   rejfreq_hr(is) = mean(abs(tstat_hr) > cv_normal);    
   
   % Compute Conley 
   for ib = 1:nbvec
       C = squeeze(CMat(:,:,ib));
       Cxe = C*xehat;
       s_conley = sqrt(mean(Cxe.^2))';
       tstat_conley = sqrt(n)*bx./s_conley;
       rejfreq_conley(is,ib) = mean(abs(tstat_conley) > cv_normal);
   end
    
       is 
end

% Save results
save([matdir mod_str '_rejfreq_scpc'],'rejfreq_scpc');
save([matdir mod_str '_rejfreq_cscpc'],'rejfreq_cscpc');
save([matdir mod_str '_rejfreq_hr'],'rejfreq_hr');
save([matdir mod_str '_rejfreq_conley'],'rejfreq_conley');

% ----- Functions -----
function[a,cbar] = draw_a(cbar,distmat,nrep_e)
  n = size(distmat,1);
  Sigma_a = exp(-cbar*distmat);
  S = sqrt_psd(Sigma_a);
  a = S*randn(n,nrep_e);
end


function[X,Z,e] = draw_XZe(a,s_draw)
  n = size(a,1);
  Z = NaN;
  e = a;
  [~,ii]=sort(s_draw);
  X = zeros(n,1);
  n1 = floor(0.85*n);
  X(n1:n) = 1;
  X = -(X - mean(X));
  X(ii) = X;
end


        

