% mc_ppt.m
% August 19, 2015, MWW
% Monte Carlo Experiment 1
% ----------- 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
 
  % ---- Output file
  str_out = [outdir 'mc_ppt.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
     scl = 1;
  
     % -- DGP Parameters 
     T_mc = 12*65;  % 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; 
     
   
      % -- Matrices to Save Result
      % -- Overall --
      bayes_save = zeros(n_clevels,1);
      MN_save = zeros(n_clevels,1);
   
      for i_mc = 1:N1_mc
        % Select d
        % Generate Data
        [u,uhat_oos] = dgp_ppt(T_mc,T_oos,N2_mc);
     
        % 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,uhat_oos,0,0,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);
     end;
 
     N = N1_mc*N2_mc;
     bc = bayes_save/N;
     mnc = MN_save/N;
   
   
     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,'   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');  

type(str_out);

end;