% Aug 18, 2014, MWW
function [Bayes_pred_sets,MLE_pred_sets,I0_pred_sets,MN_pred_sets] ... 
         = lr_pred_sp500(u,q,hfrac_vec,clevel_vec,bcd_mat,sigma_save,gdist,w_mat,lfd_mat,cvalue_mat);

%This function calculates the long-run predictions

% -- 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
 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;
 
 % Find Column of BCD Mat corresponding to I(0) model
 tmp = (bcd_mat == 0);
 tmp = sum(tmp,2);
 [tmp,i_i0] = max(tmp);
  
 % ------------------- 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
   [tmp,i_mle] = max(f_x);
      
   % 3.c.i: Compute Predictive Density for y 
   Bayes_pred_sets_Ys = NaN*zeros(2,n_clevels,n_horizons);
   MLE_pred_sets_Ys = NaN*zeros(2,n_clevels,n_horizons);
   I0_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
             
             sigma_xx= sigma_save(1:q,1:q,i_mle);
             sigma_xy = sigma_save(1:q,q+ih,i_mle);
             sigma_yy = sigma_save(q+ih,q+ih,i_mle);
             sigma_xx_inv = inv(sigma_xx);
             qf = Xs'*sigma_xx_inv*Xs;
             var_x = qf/q;
             reg_coef = sigma_xx_inv*sigma_xy;
             reg_xs = Xs'*reg_coef;
             sigma_uu = sigma_yy-sigma_xy'*reg_coef;
             for ic = 1:n_clevels;
               pct_lower = (1-clevel_vec(ic))/2;
               pct_upper = 1-pct_lower;
               t_lower = tinv(pct_lower,q);
               t_upper = tinv(pct_upper,q);
               MLE_pred_sets_Ys(1,ic,ih) = reg_xs+sqrt(var_x*sigma_uu)*t_lower; 
               MLE_pred_sets_Ys(2,ic,ih) = reg_xs+sqrt(var_x*sigma_uu)*t_upper;
             end;             
             
             sigma_xx= sigma_save(1:q,1:q,i_i0);
             sigma_xy = sigma_save(1:q,q+ih,i_i0);
             sigma_yy = sigma_save(q+ih,q+ih,i_i0);
             sigma_xx_inv = inv(sigma_xx);
             qf = Xs'*sigma_xx_inv*Xs;
             var_x = qf/q;
             reg_coef = sigma_xx_inv*sigma_xy;
             reg_xs = Xs'*reg_coef;
             sigma_uu = sigma_yy-sigma_xy'*reg_coef;
             for ic = 1:n_clevels;
               pct_lower = (1-clevel_vec(ic))/2;
               pct_upper = 1-pct_lower;
               t_lower = tinv(pct_lower,q);
               t_upper = tinv(pct_upper,q);
               I0_pred_sets_Ys(1,ic,ih) = reg_xs+sqrt(var_x*sigma_uu)*t_lower; 
               I0_pred_sets_Ys(2,ic,ih) = reg_xs+sqrt(var_x*sigma_uu)*t_upper;
             end; 
             
     end;
     
     Bayes_pred_sets = xbar + sqrt(X'*X)*Bayes_pred_sets_Ys;
     MLE_pred_sets = xbar + sqrt(X'*X)*MLE_pred_sets_Ys;
     I0_pred_sets = xbar + sqrt(X'*X)*I0_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;    