% Aug 30, 2015, MWW
function [Bayes_pred_sets,MN_pred_sets,like_bcd,post_bcd,sigma_save] ... 
         = lr_pred_set(u,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)        
%   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_pred_sets:  This contains the lower and upper endpoints of the
%                     equal-tail Bayes predictive sets.  The dimension is
%                        Dimension: 2 x (n_clevels) x (n_horizons)
%                     where (1,.,.) corresponds to the lower endpoint and
%                           (2,.,.) corresponds to the upper endpoint.
%   MN_pred_sets:  This contains the lower and upper endpoints of the
%                     equal-tail MN-predictive sets.  The dimension is
%                        Dimension: 2 x (n_clevels) x (n_horizons)
%                     where (1,.,.) corresponds to the lower endpoint and
%                           (2,.,.) corresponds to the upper endpoint. 
%                     Note: missing values are returned if LFD is not
%                     supplied to the function
% like_bcd:        This contains the likelihood values for each value of
%                  bcd. The dimension is 
%                      Dimension: (n_bcd) x 1
% post_bcd:        This is the posterior distribution for the bcd values
%                       Dimension: (n_bcd) x 1
% sigma_save:      This returns the value of the matrix sigma_save.  If
%                  sigma_save is input, its value on the return is
%                  unaltered.  If sigma_save is not input, the value
%                  internally computed by the function is returned. (This
%                  can be useful for later calls to the function.
% ier_sigma_save:  This matrix contains error codes for the computation of
%                  the Sigma matrice for each value of bcd. 0 indicates no
%                  error. 1 indicates a non-pd Sigma, which is corrected by
%                  taking the absolute value of tbe offending eigenvalue. 
%                    Dimension: (n_bcd) x 1
%
%
% -- 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
  
   % 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
 
 %------------------- Step 2: Compute Covariance Matrix for each value of hfrac
 if size(sigma_save,1) == 1;  % Compute sigma_save 
     sigma_save = NaN*zeros(q+n_horizons,q+n_horizons,n_bcd);
     fprintf('Computing %5i Sigma Matrices \n',n_bcd);
     for i = 1:n_bcd;
         b = bcd_mat(i,1);
         c = bcd_mat(i,2);
         d = bcd_mat(i,3);
         tic;
         [sigma] = Sigma_Compute(b,c,d,q,hfrac_vec);
         sigma_save(:,:,i)=sigma;
         toc_save = toc;
         if i == 1;
             fprintf('  Time per matrix %6.3f seconds \n',toc_save);
             fprintf('    Approximate total %6.2f minutes \n',toc_save*n_bcd/60);
         end;
     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 
   Bayes_pred_sets_Ys = NaN*zeros(2,n_clevels,n_horizons);
    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;
     
     % 3.d.ii: Compute Equal-Tail Prediction Sets for u
     Bayes_pred_sets = xbar + sqrt(X'*X)*Bayes_pred_sets_Ys;
     
     
     % Construct MN_prediction_sets -- find endpoints of set 
      MN_pred_sets_Ys = NaN*zeros(2,n_clevels,n_horizons);
      MN_pred_sets = NaN*zeros(2,n_clevels,n_horizons);
      if size(lfd_mat,1) > 1;
           % 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;
           
           tolerance = .001*sqrt(q);
           for ih = 1:n_horizons  % Compute prediction set for each horizon
             for ic = 1:n_clevels;
               w_weight = w_mat(:,ih,ic);
               lamdist = lfd_mat(:,ih,ic);
               cvalue = cvalue_mat(ih,ic);
               Bayes_cl = Bayes_pred_sets_Ys(1,ic,ih);
               Bayes_cu = Bayes_pred_sets_Ys(2,ic,ih);
               Bayes_width = cu - cl;
               % -------------  Lower -------------- 
               yl = Bayes_cl;
               a_in = A_check(yl,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
               if a_in == 0;
                 MN_pred_sets_Ys(1,ic,ih) = Bayes_cl;
               else;               
                 % Find yll so that a(yll) = 0
                 yll = yl;           
                 while a_in == 1;
                    yll = yll-Bayes_width;
                    a_in = A_check(yll,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
                 end;
                 % A(yl) = 1; A(yll) = 0; Update
                 while (yl-yll) > tolerance;   
                   ymid = (yl+yll)/2;
                   a_in = A_check(ymid,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
                   if a_in == 0;
                    yll = ymid;
                   else;
                    yl = ymid;
                   end;
                 end;
                 MN_pred_sets_Ys(1,ic,ih) = yl; 
               end;
               % -------------  Upper -------------- 
               yu = Bayes_cu;
               a_in = A_check(yu,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
               if a_in == 0;
                 MN_pred_sets_Ys(2,ic,ih) = Bayes_cu;
               else;
                 % Find yuu so that a(yuu) = 1
                 yuu = yu;           
                 while a_in == 1;
                    yuu = yuu + Bayes_width;
                    a_in = A_check(yuu,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
                 end;
                 % A(yu) = 1; A(yuu) = 0; Update
                 while (yuu-yu) > tolerance;    
                   ymid = (yuu+yu)/2;
                   a_in = A_check(ymid,Xs,sigma_save,ih,f_xg,w_weight,lamdist,cvalue);
                   if a_in == 0;
                    yuu = ymid;
                   else;
                    yu = ymid;
                   end;
                 end;
                 MN_pred_sets_Ys(2,ic,ih) = yu; 
               end;
            
             end;  %ic    
               
           end; %ih   
        MN_pred_sets = xbar + sqrt(X'*X)*MN_pred_sets_Ys; 
           
      end; %size lfd_mat > 1                   