% Aug 9, 2104, MWW
function [Bayes_num,MN_num] ... 
         = lr_pred_set_mc_split(u,uhat_mat,n0,n1,q,hfrac_vec,clevel_vec,bcd_mat,sigma_save,gdist,w_mat,lfd_mat,cvalue_mat);
%This function calculates the long-run prediction in Mueller and Watson's
%Measuring Uncertainty about Long-Run Predictions
%The steps listed in the function correspond to the steps listed in
%Appendix of that paper
% ------------ Input Parameters
%   u:  Data vector ... missing values are not allowed
%            Dimension: (Tx1)      
%   uhat_mat: Out-of-sample values of uhat: (n_mc x n_horizons)  
%   q:  Number of low-frequency cosine transformations to use
%            Dimension: scalar
%   hfrac_vec: This vector contains the forecast horizons expressed as
%              fractions of T. Because this a vector, multiple horizons can
%              be computed with one call to the function
%            Dimension: (n_horizons x 1)
%   clevel_vec: This vector contains the "confidence" levels. that is the
%               coverage level of the prediction sets (e.g., 50%, 80%, 90%)
%            Dimension: (n_clevels x 1)
%   bcd_mat: This matrix contains the grid of values of 'bcd' ('theta' in the
%            paper) under consideration.
%             Dimension: (n_bcd x 1)
%   sigma_save:  This 3-dimensional matrix contains the "sigma" matrices for each row of bcd_mat 
%                and including all forecast horizons.
%                  If sigma_save = scalar, then the sigma matrices are
%                  computed by the program.
%                  Otherwise the dimension of sigma_mat is:
%                  Demension:(q+n_horizons)x(q_n_horizons)x(n_bcd);
%   gdist: This is (a) the Bayes prior for 'bcd'. It is ALSO the weighting
%          matrix for the set volume used in the frequentist calcuation. In
%          both cases this is called "Gamma" in the paper.
%              Dimesnion: n_bcd x 1
%   w_mat: This is a matrix of for weighting the length(over bcd) for frequentist intervals. 
%            As discussed in the paper, the w depends on q,"gdist',
%            horizon, and coverage level.  If lfd_mat = scalar, the MN
%            prediction sets are not computed.  Otherwise the w_mat is a
%            3 dimensional matrix with 
%              Dimension: (n_bcd) x (n_horizons) x (n_clevels)
%   lfd_mat: This is a matrix of least favorable distributions (over bcd). 
%            As discussed in the paper, the LFD depends on q,"gdist',
%            horizon, and coverage level.  If lfd_mat = scalar, the MN
%            prediction sets are not computed.  Otherwise the lfd_mat is a
%            3 dimensional matrix with 
%              Dimension: (n_bcd) x (n_horizons) x (n_clevels)
%   cvalue_mat: This is a matrix of critical values for each of the least favorable 
%            distributions. If the lfd_mat is not given, then this is not used (so it 
%            can be an arbitrary scalar. Otherwise, its dimension is
%               Dimension: (n_horizons) x (n_clevels) 
% --------------- Output parameters
%
%   Bayes_num:  number of times Bayes set includes uhat_vec
%                 This is n_clevels x n_horizons
%   MN_num: number of times Bayes set includes uhat_vec
%                 This is n_clevels x n_horizons
% -- Parameter Names
  T = size(u,1);
  n_horizons = size(hfrac_vec,1);  % Number of Forecast Horizons
  n_clevels = size(clevel_vec,1);  % Number of "Confidence" levels
  n_bcd = size(bcd_mat,1);         % number of 'bcd' parameter values
  n_mc = size(uhat_mat,1);         % number of oos mc values of uhat
  
   % Check for missing values, sizes of arrays and so forth
  if isnan(u) > 0;
      error('missing values in data vector');
  end;
  if size(gdist,1) ~= n_bcd;
      error('gdist has wrong dimension');
  end;
 
 %------------------ Step 1: Compute xbar(1:T), Cosine Transform and Standardised Cosine transform
 psi = psi_compute(T,q);
 xbar = mean(u);          % Sample Mean
 X = psi'*(u - xbar);     % Cosine Transformaion
 Xs = X/sqrt(X'*X);       % Normalize Cos Transformation
 
 % Transform uhat_mat to Ys_mat
 Ys_mat = (uhat_mat-xbar)/sqrt(X'*X);   % Normalize 
 
 %------------------- Step 2: Compute Covariance Matrix for each value of hfrac
 ier_sigma_save = zeros(n_bcd,1);
 if size(sigma_save,1) == 1;  % Compute sigma_save 
     sigma_save = NaN*zeros(q+n_horizons,q+n_horizons,n_bcd);
     for i = 1:n_bcd;
         b = bcd_mat(i,1);
         c = bcd_mat(i,2);
         d = bcd_mat(i,3);
         [sigma,ier] = Sigma_Compute(b,c,d,q,hfrac_vec);
         tmp = eig(sigma);
         sigma_save(:,:,i)=sigma;
         ier_sigma_save(i) = ier;
     end;
 end;
  
 % ------------------- Step 3: Bayes Analysis
   % 3.a: gdist is prior 
   % 3.b: Compute likelihood for each value of theta= 
   f_x = NaN*zeros(n_bcd,1);
   xs_qf = NaN*zeros(n_bcd,1);
   for i=1:n_bcd;
       sigma = sigma_save(1:q,1:q,i);
       sigma_inv = inv(sigma);
       f_x(i) = den_invariant(Xs,sigma_inv);
       xs_qf(i) = Xs'*sigma_inv*Xs;
   end;
   like_bcd = f_x;                     % Likelihood values
   f_marg_x = gdist'*f_x;              % Marginal Likelihood
   post_bcd = (gdist.*f_x)/f_marg_x;   % Posterior
      
   % 3.c.i: Compute Predictive Density for y over a grid of values of y
    for ih = 1:n_horizons  % Compute prediction set for each horizon
             reg_xs = NaN*zeros(n_bcd,1);    % Regression Coefficients times Xs
             sigma_uu = NaN*zeros(n_bcd,1);  % Variance of regression error
             for i = 1:n_bcd; 
                sigma_xx= sigma_save(1:q,1:q,i); 
                sigma_xy = sigma_save(1:q,q+ih,i);
                sigma_yy = sigma_save(q+ih,q+ih,i);
                sigma_xx_inv = inv(sigma_xx);
                reg_coef = sigma_xx_inv*sigma_xy;
                reg_xs(i) = Xs'*reg_coef;
                sigma_uu(i) = sigma_yy-sigma_xy'*reg_coef;
             end;
                   
             tolerance = .0001;
             for ic = 1:n_clevels;
               pct_lower = (1-clevel_vec(ic))/2;
               pct_upper = 1-pct_lower;
               [cl,pct_actual] = t_mixture(pct_lower,post_bcd,reg_xs,sigma_uu,xs_qf,q,tolerance);
               [cu,pct_actual] = t_mixture(pct_upper,post_bcd,reg_xs,sigma_uu,xs_qf,q,tolerance); 
               Bayes_pred_sets_Ys(1,ic,ih) = cl;
               Bayes_pred_sets_Ys(2,ic,ih) = cu;
             end
     end;
     
     Bayes_in_mat = NaN*zeros(n_mc,n_clevels,n_horizons);
           for ih = 1:n_horizons  % Compute prediction set for each horizon
               for ic = 1:n_clevels;
                  B_lt = Ys_mat(:,ih) > Bayes_pred_sets_Ys(1,ic,ih);
                  B_gt = Ys_mat(:,ih) < Bayes_pred_sets_Ys(2,ic,ih);
                  Bayes_in_mat(:,ic,ih) = B_lt.*B_gt;
               end;
           end;
     
     
     % Compute joint density of Ys and Xs over a grid of values of Ys; Choose Grid points
     %         over support that is three times as wide as Bayes sets 
      f_yx_save = NaN*zeros(n_mc,n_bcd,n_horizons);   % Joint density for each horizon
      for ih = 1:n_horizons  % Compute prediction set for each horizon
         Ys_grid = Ys_mat(:,ih);   % grid of values for Ys
         Xs_Ys_grid = [repmat(Xs,1,size(Ys_grid,1)); Ys_grid'];   % grid of values for Xs and Ys (Xs held at sample values)
         f_yx = NaN*zeros(size(Ys_grid,1),n_bcd);  % for joint density
         for i = 1:n_bcd; 
                sigma = zeros(q+1,q+1);
                sigma(1:q,1:q) = sigma_save(1:q,1:q,i); 
                sigma(1:q,q+1) = sigma_save(1:q,q+ih,i);
                sigma(q+1,1:q) = sigma_save(q+ih,1:q,i);
                sigma(q+1,q+1) = sigma_save(q+ih,q+ih,i);
                sigma_inv = inv(sigma);
                f_yx(:,i) = den_invariant(Xs_Ys_grid,sigma_inv);
         end;
         f_yx_save(:,:,ih) = f_yx;                          % Save Joint Density
      end;

% ------------------- Step 4: Construct MN-Frequentist Prediction Sets
      % Compute (0,1) for Ys_mat in A prediction set
           % Compute Denominator g(x) = E(sqrt(X'X)|Xs) for each value of theta
           g = NaN*zeros(n_bcd,1);
           c = sqrt(2)*gamma((q+1)/2)/gamma(q/2);
           for i=1:n_bcd;
             sigma = sigma_save(1:q,1:q,i);
             sigma_inv = inv(sigma);
             qf = Xs'*sigma_inv*Xs;
             g(i) = c/sqrt(qf);
           end;
           f_xg = f_x.*g;
           A_in_mat = NaN*zeros(n_mc,n_clevels,n_horizons);
           for ih = 1:n_horizons  % Compute prediction set for each horizon
               f_yx = f_yx_save(:,:,ih);
               for ic = 1:n_clevels;
                  w_weight = w_mat(:,ih,ic);
                  lamdist = lfd_mat(:,ih,ic);
                  cvalue = cvalue_mat(ih,ic);
                  f_yx_lam = f_yx*lamdist;
                  f_xg_w = f_xg'*w_weight;
                  A_in_mat(:,ic,ih) = f_yx_lam >= cvalue*f_xg_w;
               end;
           end;
     
           MN_in_mat = NaN*zeros(n_mc,n_clevels,n_horizons);
           for ih = 1:n_horizons  % Compute prediction set for each horizon
               for ic = 1:n_clevels;
                  B = Bayes_in_mat(:,ic,ih);
                  A = A_in_mat(:,ic,ih);
                  MN_in_mat(:,ic,ih) = (A+B)>0;
               end;
           end;
           Bayes_num = NaN*zeros(n_clevels,n_horizons,3);
           MN_num = NaN*zeros(n_clevels,n_horizons,3);
           for ih = 1:n_horizons  % Compute prediction set for each horizon
               for ic = 1:n_clevels;
                  Bayes_num(ic,ih,1) = sum(Bayes_in_mat(:,ic,ih));
                  MN_num(ic,ih,1) = sum(MN_in_mat(:,ic,ih));
                  Bayes_num(ic,ih,2) = sum(Bayes_in_mat(1:n0,ic,ih));
                  MN_num(ic,ih,2) = sum(MN_in_mat(1:n0,ic,ih));
                  Bayes_num(ic,ih,3) = sum(Bayes_in_mat(n0+1:n0+n1,ic,ih));
                  MN_num(ic,ih,3) = sum(MN_in_mat(n0+1:n0+n1,ic,ih));
               end;
           end;
                       
end
