% mc_break_1_discrete.m
% August 31, 2015, MWW
% Monte Carlo Experiment 2
% ----------- Step 1: Set up parameters and so forth for run
clear all;
small = 1.0e-10;
  % -- File Directories 
  
  outdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/out/';
  figdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/fig/';
  matdir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/mat/';;
  mat_save_dir = '/Users/mwatson/Dropbox/VIC/ddisk_Sept_2015/Matlab/mat_save/';;  % Location for Saving Various Files
   
  str_out = [outdir 'mc_break_2_discrete.out'];
  fid_out = fopen(str_out,'w');
  
  q_vec = [12]';
  
  for iq = 1: size(q_vec,1);
  
  q = q_vec(iq);  % Value of q to use 
  
  % Fix the Seed for the random number generator
  s_seed = 265497;
  rng(s_seed);
  
  % -- Parameter
  clevel_vec = [0.67 0.90]';      % Confidence Level
  rho = 0.50;                     % Horizon fraction
  rho_ar1 = 0.98;                 % AR1 coef
  scl = 0.46;                     % LR STD of u-rho_ar*u(-1)
  
  
  % -- DGP Parameters 
  T_mc = 260;  % In-sample Size
  N1_mc = 5000;  % Number of draws of in-sample data
  N2_mc = 500; % Number of draws of OOS data, given each draw of in-sample data 
  
 % -- Parameters to govern calculation of covariance matrices, grid of BCD values and so forth
  i_first = 0;  % = 1, the first time this program is used.  In this case, the program sets up various matrices
                %      (cov. matrices, BCD values, horizon fractions, etc.) that can be used for future calls.
               
% ---------- Step 2: Set up BCD parameters, Sigma matrices and so forth
  if i_first == 1;
    lr_main_setup_matrices;
  end; 
  
% ---------- Step 3: ---- Read in variables and sigma matrices
 %  ----- Read in variables
      str_tmp = 'hfrac_vec'; load([matdir str_tmp '_' int2str(q)]);    
      str_tmp = 'bcd_mat'; load([mat_save_dir str_tmp '_' int2str(q)]);    
      str_tmp = 'gdist'; load([mat_save_dir str_tmp '_' int2str(q)]);  
      str_tmp = 'cvalue_mat'; load([mat_save_dir str_tmp '_' int2str(q)]); 
      str_tmp = 'lfd_mat'; load([mat_save_dir str_tmp '_' int2str(q)]); 
      str_tmp = 'w_mat'; load([mat_save_dir str_tmp '_' int2str(q)]);
      str_tmp = 'sig_matlab_save'; load([mat_save_dir str_tmp '_' int2str(q)]);  
   % --- Some Dimensions
   n_bcd = size(bcd_mat,1);
   n_clevels = size(clevel_vec,1);
   n_horizons = size(hfrac_vec,1);
   
    % Sigma Matrices;
      sigma_save = sig_matlab_save;
      
      % Save components relevant for horizon
      ii_rho = hfrac_vec == rho;
      trend = (1:1:size(hfrac_vec,1))';
      irho = trend(ii_rho==1);
      cvalue_mat = cvalue_mat(irho,:);
      lfd_mat = lfd_mat(:,irho,:);
      w_mat = w_mat(:,irho,:);
      sig_tmp = sigma_save;
      sigma_save = zeros(q+1,q+1,n_bcd);
      sigma_save(1:q,1:q,:) = sig_tmp(1:q,1:q,:);
      sigma_save(1:q,q+1,:) = sig_tmp(1:q,q+irho,:);
      sigma_save(q+1,1:q,:) = sig_tmp(q+irho,1:q,:);
      sigma_save(q+1,q+1,:) = sig_tmp(q+irho,q+irho,:);
      hfrac_vec = hfrac_vec(irho);
   
   % ---------
   T_oos = ceil(rho*T_mc);
   T_total = T_mc + T_oos; 
  
   p1 = 1/40;
   sz_vec1 = [0.50;1.0];
   p2 = 1/260;
   sz_vec2 = [1.0;2.0];
   
   p_vec = [p1;p2];
   sz_mat = [sz_vec1 sz_vec2];
  
   s = AR1_cov(T_total,rho_ar1);
   s11 = s(1:T_mc,1:T_mc);
   s21 = s(T_mc+1:T_total,1:T_mc);
   s22 = s(T_mc+1:T_total,T_mc+1:T_total);
   reg_coef = s21*inv(s11);
   s22_1 = s22-reg_coef*s21';
   ch11 = chol(s11);
   ch22_1 = chol(s22_1);
  
   
for ip = 1:size(p_vec,1);
    p = p_vec(ip,1);
    sz_vec = sz_mat(:,ip);
 
  for isz = 1:size(sz_vec);
   sz = sz_vec(isz);
   sz_scl = sz/scl;
   
   
   n0 = 0;
   n0_0 = 0;
   n0_1 = 0;
   n1 = 0;
   n1_0 = 0;
   n1_1 = 0;
   
   
    % -- Matrices to Save Result
   % -- Overall --
   bayes_save = zeros(n_clevels,1);
   MN_save = zeros(n_clevels,1);
   
   % No in-sample break
   bayes_save_0 = zeros(n_clevels,1);
   MN_save_0 = zeros(n_clevels,1);
   
   % In-sample Break
   bayes_save_1 = zeros(n_clevels,1);
   MN_save_1 = zeros(n_clevels,1);
   
   % In-sample break and no OOS
   bayes_save_1_0 = zeros(n_clevels,1);
   MN_save_1_0 = zeros(n_clevels,1);
   
   % In-sample break and OOS
   bayes_save_1_1 = zeros(n_clevels,1);
   MN_save_1_1 = zeros(n_clevels,1);
   
   % No In-sample break and no OOS
   bayes_save_0_0 = zeros(n_clevels,1);
   MN_save_0_0 = zeros(n_clevels,1);
   
   % No In-sample break and OOS
   bayes_save_0_1 = zeros(n_clevels,1);
   MN_save_0_1 = zeros(n_clevels,1);
   
   
   
   for i_mc = 1:N1_mc
     % Select d
     % Generate Data
     [u,uhat_oos,bis,boosvec,mu_iis_T,muv_oos_T] = dgp_gen_I1b_discrete(ch11,ch22_1,reg_coef,p,sz,N2_mc);
     uhat_0 = uhat_oos(boosvec == 0);
     uhat_1 = uhat_oos(boosvec == 1);
     noos_0 = sum(boosvec == 0);
     noos_1 = sum(boosvec == 1);
     ubar = NaN;
     if noos_0 > 0;
      ubar = [ubar;uhat_0];
     end;
     if noos_1 > 0;
      ubar = [ubar;uhat_1];
     end;
     ubar = ubar(2:end);
     
    
     % Compute Prediction Set
     bayes_num = zeros(n_clevels,3);
     MN_num = zeros(n_clevels,3);
     [bayes_num1,MN_num1] ... 
         = lr_pred_set_mc_split(u,ubar,noos_0,noos_1,q,hfrac_vec,clevel_vec,bcd_mat,sigma_save,gdist,w_mat,lfd_mat,cvalue_mat);
     bayes_num(:,:) = bayes_num1(:,1,:);
     MN_num(:,:) = MN_num1(:,1,:);
     bayes_save = bayes_save + bayes_num(:,1);
     MN_save = MN_save+MN_num(:,1);
     if bis == 0;
         n0 = n0+1;
         bayes_save_0 = bayes_save_0 + bayes_num(:,1);
         MN_save_0 = MN_save_0+MN_num(:,1);
         n0_0 = n0_0 + noos_0;
         bayes_save_0_0 = bayes_save_0_0 + bayes_num(:,2);
         MN_save_0_0 = MN_save_0_0+MN_num(:,2);
         n0_1 = n0_1 + noos_1;
         bayes_save_0_1 = bayes_save_0_1 + bayes_num(:,3);
         MN_save_0_1 = MN_save_0_1+MN_num(:,3); 
     end;
     if bis == 1;
         n1 = n1+1;
         bayes_save_1 = bayes_save_1 + bayes_num(:,1);
         MN_save_1 = MN_save_1+MN_num(:,1);
         n1_0 = n1_0 + noos_0;
         bayes_save_1_0 = bayes_save_1_0 + bayes_num(:,2);
         MN_save_1_0 = MN_save_1_0+MN_num(:,2);
         n1_1 = n1_1 + noos_1;
         bayes_save_1_1 = bayes_save_1_1 + bayes_num(:,3);
         MN_save_1_1 = MN_save_1_1+MN_num(:,3); 
     end;     
      
  end;
 
   N = N1_mc*N2_mc;
   N0 = n0*N2_mc;
   N1 = n1*N2_mc;
   frac_0 = N0/N;
   frac_1 = N1/N;
   frac_0_0 = n0_0/N;
   frac_1_0 = n1_0/N;
   frac_0_1 = n0_1/N;
   frac_1_1 = n1_1/N;
 
   bc = bayes_save/N;
   mnc = MN_save/N;
   
   bc_0 = bayes_save_0/N0;
   mnc_0 = MN_save_0/N0;
   
   bc_1 = bayes_save_1/N1;
   mnc_1 = MN_save_1/N1;
   
   bc_0_0 = bayes_save_0_0/n0_0;
   mnc_0_0 = MN_save_0_0/n0_0;
   
   bc_0_1 = bayes_save_0_1/n0_1;
   mnc_0_1 = MN_save_0_1/n0_1;
   
   bc_1_0 = bayes_save_1_0/n1_0;
   mnc_1_0 = MN_save_1_0/n1_0;
   
   bc_1_1 = bayes_save_1_1/n1_1;
   mnc_1_1 = MN_save_1_1/n1_1;
   
   
   fprintf(fid_out,'\n\n\n');
   fprintf(fid_out,[' Date and Time: ' datestr(now) '\n']);
   fprintf(fid_out,' q: %3i \n',q);
   fprintf(fid_out,' N1_mc = %5.0f\n',N1_mc);
   fprintf(fid_out,' N2_mc = %6.0f\n',N2_mc);
   fprintf(fid_out,'  Discrete Break in I(0) mean \n');
   fprintf(fid_out,'   Probability of break = %6.4f \n',p);
   fprintf(fid_out,'   Size of break = %6.3f \n',sz);
   fprintf(fid_out,'   Scale Factor = %6.3f \n',scl);
   fprintf(fid_out,'\n');
   fprintf(fid_out,'   Overall Coverage \n');
   fprintf(fid_out,'     Bayes Sets %6.3f,%6.3f\n',bc(1),bc(2));
   fprintf(fid_out,'     MN Sets    %6.3f,%6.3f\n',mnc(1),mnc(2));
   fprintf(fid_out,'\n');
     
     fprintf(fid_out,'     In-Sample Break = 0 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_0); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_0(1),bc_0(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_0(1),mnc_0(2));
     fprintf(fid_out,'\n');
     
     fprintf(fid_out,'     In-Sample Break = 1 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_1); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_1(1),bc_1(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_1(1),mnc_1(2));
     fprintf(fid_out,'\n');
   
     fprintf(fid_out,'     In-Sample Break = 0, Out of Sample = 0 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_0_0); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_0_0(1),bc_0_0(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_0_0(1),mnc_0_0(2));
     fprintf(fid_out,'\n');
   
     fprintf(fid_out,'     In-Sample Break = 0, Out of Sample = 1 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_0_1); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_0_1(1),bc_0_1(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_0_1(1),mnc_0_1(2));
     fprintf(fid_out,'\n');
     
     fprintf(fid_out,'     In-Sample Break = 1, Out of Sample = 0 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_1_0); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_1_0(1),bc_1_0(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_1_0(1),mnc_1_0(2));
     fprintf(fid_out,'\n');
   
     fprintf(fid_out,'     In-Sample Break = 1, Out of Sample = 1 \n');
     fprintf(fid_out,'     Fraction of total = %6.3f\n',frac_1_1); 
     fprintf(fid_out,'       Bayes Sets %6.3f,%6.3f\n',bc_1_1(1),bc_1_1(2));
     fprintf(fid_out,'       MN Sets    %6.3f,%6.3f\n',mnc_1_1(1),mnc_1_1(2));
     fprintf(fid_out,'\n');
   
 
 end;
 
end;
   

type(str_out);

end;